template_lapack_sterf.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027  
00028  /* This file belongs to the template_lapack part of the Ergo source 
00029   * code. The source files in the template_lapack directory are modified
00030   * versions of files originally distributed as CLAPACK, see the
00031   * Copyright/license notice in the file template_lapack/COPYING.
00032   */
00033  
00034 
00035 #ifndef TEMPLATE_LAPACK_STERF_HEADER
00036 #define TEMPLATE_LAPACK_STERF_HEADER
00037 
00038 #include "template_lapack_common.h"
00039 
00040 template<class Treal>
00041 int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, 
00042         integer *info)
00043 {
00044 /*  -- LAPACK routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        June 30, 1999   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DSTERF computes all eigenvalues of a symmetric tridiagonal matrix   
00054     using the Pal-Walker-Kahan variant of the QL or QR algorithm.   
00055 
00056     Arguments   
00057     =========   
00058 
00059     N       (input) INTEGER   
00060             The order of the matrix.  N >= 0.   
00061 
00062     D       (input/output) DOUBLE PRECISION array, dimension (N)   
00063             On entry, the n diagonal elements of the tridiagonal matrix.   
00064             On exit, if INFO = 0, the eigenvalues in ascending order.   
00065 
00066     E       (input/output) DOUBLE PRECISION array, dimension (N-1)   
00067             On entry, the (n-1) subdiagonal elements of the tridiagonal   
00068             matrix.   
00069             On exit, E has been destroyed.   
00070 
00071     INFO    (output) INTEGER   
00072             = 0:  successful exit   
00073             < 0:  if INFO = -i, the i-th argument had an illegal value   
00074             > 0:  the algorithm failed to find all of the eigenvalues in   
00075                   a total of 30*N iterations; if INFO = i, then i   
00076                   elements of E have not converged to zero.   
00077 
00078     =====================================================================   
00079 
00080 
00081        Test the input parameters.   
00082 
00083        Parameter adjustments */
00084     /* Table of constant values */
00085      integer c__0 = 0;
00086      integer c__1 = 1;
00087      Treal c_b32 = 1.;
00088     
00089     /* System generated locals */
00090     integer i__1;
00091     Treal d__1, d__2, d__3;
00092     /* Local variables */
00093      Treal oldc;
00094      integer lend, jtot;
00095      Treal c__;
00096      integer i__, l, m;
00097      Treal p, gamma, r__, s, alpha, sigma, anorm;
00098      integer l1;
00099      Treal bb;
00100      integer iscale;
00101      Treal oldgam, safmin;
00102      Treal safmax;
00103      integer lendsv;
00104      Treal ssfmin;
00105      integer nmaxit;
00106      Treal ssfmax, rt1, rt2, eps, rte;
00107      integer lsv;
00108      Treal eps2;
00109 
00110 
00111     --e;
00112     --d__;
00113 
00114     /* Function Body */
00115     *info = 0;
00116 
00117 /*     Quick return if possible */
00118 
00119     if (*n < 0) {
00120         *info = -1;
00121         i__1 = -(*info);
00122         template_blas_erbla("STERF ", &i__1);
00123         return 0;
00124     }
00125     if (*n <= 1) {
00126         return 0;
00127     }
00128 
00129 /*     Determine the unit roundoff for this environment. */
00130 
00131     eps = template_lapack_lamch("E", (Treal)0);
00132 /* Computing 2nd power */
00133     d__1 = eps;
00134     eps2 = d__1 * d__1;
00135     safmin = template_lapack_lamch("S", (Treal)0);
00136     safmax = 1. / safmin;
00137     ssfmax = template_blas_sqrt(safmax) / 3.;
00138     ssfmin = template_blas_sqrt(safmin) / eps2;
00139 
00140 /*     Compute the eigenvalues of the tridiagonal matrix. */
00141 
00142     nmaxit = *n * 30;
00143     sigma = 0.;
00144     jtot = 0;
00145 
00146 /*     Determine where the matrix splits and choose QL or QR iteration   
00147        for each block, according to whether top or bottom diagonal   
00148        element is smaller. */
00149 
00150     l1 = 1;
00151 
00152 L10:
00153     if (l1 > *n) {
00154         goto L170;
00155     }
00156     if (l1 > 1) {
00157         e[l1 - 1] = 0.;
00158     }
00159     i__1 = *n - 1;
00160     for (m = l1; m <= i__1; ++m) {
00161         if ((d__3 = e[m], absMACRO(d__3)) <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) * 
00162                 template_blas_sqrt((d__2 = d__[m + 1], absMACRO(d__2))) * eps) {
00163             e[m] = 0.;
00164             goto L30;
00165         }
00166 /* L20: */
00167     }
00168     m = *n;
00169 
00170 L30:
00171     l = l1;
00172     lsv = l;
00173     lend = m;
00174     lendsv = lend;
00175     l1 = m + 1;
00176     if (lend == l) {
00177         goto L10;
00178     }
00179 
00180 /*     Scale submatrix in rows and columns L to LEND */
00181 
00182     i__1 = lend - l + 1;
00183     anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
00184     iscale = 0;
00185     if (anorm > ssfmax) {
00186         iscale = 1;
00187         i__1 = lend - l + 1;
00188         template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, 
00189                 info);
00190         i__1 = lend - l;
00191         template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, 
00192                 info);
00193     } else if (anorm < ssfmin) {
00194         iscale = 2;
00195         i__1 = lend - l + 1;
00196         template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, 
00197                 info);
00198         i__1 = lend - l;
00199         template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, 
00200                 info);
00201     }
00202 
00203     i__1 = lend - 1;
00204     for (i__ = l; i__ <= i__1; ++i__) {
00205 /* Computing 2nd power */
00206         d__1 = e[i__];
00207         e[i__] = d__1 * d__1;
00208 /* L40: */
00209     }
00210 
00211 /*     Choose between QL and QR iteration */
00212 
00213     if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
00214         lend = lsv;
00215         l = lendsv;
00216     }
00217 
00218     if (lend >= l) {
00219 
00220 /*        QL Iteration   
00221 
00222           Look for small subdiagonal element. */
00223 
00224 L50:
00225         if (l != lend) {
00226             i__1 = lend - 1;
00227             for (m = l; m <= i__1; ++m) {
00228                 if ((d__2 = e[m], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m 
00229                         + 1], absMACRO(d__1))) {
00230                     goto L70;
00231                 }
00232 /* L60: */
00233             }
00234         }
00235         m = lend;
00236 
00237 L70:
00238         if (m < lend) {
00239             e[m] = 0.;
00240         }
00241         p = d__[l];
00242         if (m == l) {
00243             goto L90;
00244         }
00245 
00246 /*        If remaining matrix is 2 by 2, use DLAE2 to compute its   
00247           eigenvalues. */
00248 
00249         if (m == l + 1) {
00250             rte = template_blas_sqrt(e[l]);
00251             template_lapack_lae2(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
00252             d__[l] = rt1;
00253             d__[l + 1] = rt2;
00254             e[l] = 0.;
00255             l += 2;
00256             if (l <= lend) {
00257                 goto L50;
00258             }
00259             goto L150;
00260         }
00261 
00262         if (jtot == nmaxit) {
00263             goto L150;
00264         }
00265         ++jtot;
00266 
00267 /*        Form shift. */
00268 
00269         rte = template_blas_sqrt(e[l]);
00270         sigma = (d__[l + 1] - p) / (rte * 2.);
00271         r__ = template_lapack_lapy2(&sigma, &c_b32);
00272         sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
00273 
00274         c__ = 1.;
00275         s = 0.;
00276         gamma = d__[m] - sigma;
00277         p = gamma * gamma;
00278 
00279 /*        Inner loop */
00280 
00281         i__1 = l;
00282         for (i__ = m - 1; i__ >= i__1; --i__) {
00283             bb = e[i__];
00284             r__ = p + bb;
00285             if (i__ != m - 1) {
00286                 e[i__ + 1] = s * r__;
00287             }
00288             oldc = c__;
00289             c__ = p / r__;
00290             s = bb / r__;
00291             oldgam = gamma;
00292             alpha = d__[i__];
00293             gamma = c__ * (alpha - sigma) - s * oldgam;
00294             d__[i__ + 1] = oldgam + (alpha - gamma);
00295             if (c__ != 0.) {
00296                 p = gamma * gamma / c__;
00297             } else {
00298                 p = oldc * bb;
00299             }
00300 /* L80: */
00301         }
00302 
00303         e[l] = s * p;
00304         d__[l] = sigma + gamma;
00305         goto L50;
00306 
00307 /*        Eigenvalue found. */
00308 
00309 L90:
00310         d__[l] = p;
00311 
00312         ++l;
00313         if (l <= lend) {
00314             goto L50;
00315         }
00316         goto L150;
00317 
00318     } else {
00319 
00320 /*        QR Iteration   
00321 
00322           Look for small superdiagonal element. */
00323 
00324 L100:
00325         i__1 = lend + 1;
00326         for (m = l; m >= i__1; --m) {
00327             if ((d__2 = e[m - 1], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m 
00328                     - 1], absMACRO(d__1))) {
00329                 goto L120;
00330             }
00331 /* L110: */
00332         }
00333         m = lend;
00334 
00335 L120:
00336         if (m > lend) {
00337             e[m - 1] = 0.;
00338         }
00339         p = d__[l];
00340         if (m == l) {
00341             goto L140;
00342         }
00343 
00344 /*        If remaining matrix is 2 by 2, use DLAE2 to compute its   
00345           eigenvalues. */
00346 
00347         if (m == l - 1) {
00348             rte = template_blas_sqrt(e[l - 1]);
00349             template_lapack_lae2(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
00350             d__[l] = rt1;
00351             d__[l - 1] = rt2;
00352             e[l - 1] = 0.;
00353             l += -2;
00354             if (l >= lend) {
00355                 goto L100;
00356             }
00357             goto L150;
00358         }
00359 
00360         if (jtot == nmaxit) {
00361             goto L150;
00362         }
00363         ++jtot;
00364 
00365 /*        Form shift. */
00366 
00367         rte = template_blas_sqrt(e[l - 1]);
00368         sigma = (d__[l - 1] - p) / (rte * 2.);
00369         r__ = template_lapack_lapy2(&sigma, &c_b32);
00370         sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
00371 
00372         c__ = 1.;
00373         s = 0.;
00374         gamma = d__[m] - sigma;
00375         p = gamma * gamma;
00376 
00377 /*        Inner loop */
00378 
00379         i__1 = l - 1;
00380         for (i__ = m; i__ <= i__1; ++i__) {
00381             bb = e[i__];
00382             r__ = p + bb;
00383             if (i__ != m) {
00384                 e[i__ - 1] = s * r__;
00385             }
00386             oldc = c__;
00387             c__ = p / r__;
00388             s = bb / r__;
00389             oldgam = gamma;
00390             alpha = d__[i__ + 1];
00391             gamma = c__ * (alpha - sigma) - s * oldgam;
00392             d__[i__] = oldgam + (alpha - gamma);
00393             if (c__ != 0.) {
00394                 p = gamma * gamma / c__;
00395             } else {
00396                 p = oldc * bb;
00397             }
00398 /* L130: */
00399         }
00400 
00401         e[l - 1] = s * p;
00402         d__[l] = sigma + gamma;
00403         goto L100;
00404 
00405 /*        Eigenvalue found. */
00406 
00407 L140:
00408         d__[l] = p;
00409 
00410         --l;
00411         if (l >= lend) {
00412             goto L100;
00413         }
00414         goto L150;
00415 
00416     }
00417 
00418 /*     Undo scaling if necessary */
00419 
00420 L150:
00421     if (iscale == 1) {
00422         i__1 = lendsv - lsv + 1;
00423         template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], 
00424                 n, info);
00425     }
00426     if (iscale == 2) {
00427         i__1 = lendsv - lsv + 1;
00428         template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], 
00429                 n, info);
00430     }
00431 
00432 /*     Check for no convergence to an eigenvalue after a total   
00433        of N*MAXIT iterations. */
00434 
00435     if (jtot < nmaxit) {
00436         goto L10;
00437     }
00438     i__1 = *n - 1;
00439     for (i__ = 1; i__ <= i__1; ++i__) {
00440         if (e[i__] != 0.) {
00441             ++(*info);
00442         }
00443 /* L160: */
00444     }
00445     goto L180;
00446 
00447 /*     Sort eigenvalues in increasing order. */
00448 
00449 L170:
00450     template_lapack_lasrt("I", n, &d__[1], info);
00451 
00452 L180:
00453     return 0;
00454 
00455 /*     End of DSTERF */
00456 
00457 } /* dsterf_ */
00458 
00459 #endif

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7