template_lapack_stevx.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_STEVX_HEADER
00036 #define TEMPLATE_LAPACK_STEVX_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_stevx(const char *jobz, const char *range, const integer *n, Treal *
00041         d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il, 
00042         const integer *iu, const Treal *abstol, integer *m, Treal *w, 
00043         Treal *z__, const integer *ldz, Treal *work, integer *iwork, 
00044         integer *ifail, integer *info)
00045 {
00046 /*  -- LAPACK driver routine (version 3.0) --   
00047        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00048        Courant Institute, Argonne National Lab, and Rice University   
00049        June 30, 1999   
00050 
00051 
00052     Purpose   
00053     =======   
00054 
00055     DSTEVX computes selected eigenvalues and, optionally, eigenvectors   
00056     of a real symmetric tridiagonal matrix A.  Eigenvalues and   
00057     eigenvectors can be selected by specifying either a range of values   
00058     or a range of indices for the desired eigenvalues.   
00059 
00060     Arguments   
00061     =========   
00062 
00063     JOBZ    (input) CHARACTER*1   
00064             = 'N':  Compute eigenvalues only;   
00065             = 'V':  Compute eigenvalues and eigenvectors.   
00066 
00067     RANGE   (input) CHARACTER*1   
00068             = 'A': all eigenvalues will be found.   
00069             = 'V': all eigenvalues in the half-open interval (VL,VU]   
00070                    will be found.   
00071             = 'I': the IL-th through IU-th eigenvalues will be found.   
00072 
00073     N       (input) INTEGER   
00074             The order of the matrix.  N >= 0.   
00075 
00076     D       (input/output) DOUBLE PRECISION array, dimension (N)   
00077             On entry, the n diagonal elements of the tridiagonal matrix   
00078             A.   
00079             On exit, D may be multiplied by a constant factor chosen   
00080             to avoid over/underflow in computing the eigenvalues.   
00081 
00082     E       (input/output) DOUBLE PRECISION array, dimension (N)   
00083             On entry, the (n-1) subdiagonal elements of the tridiagonal   
00084             matrix A in elements 1 to N-1 of E; E(N) need not be set.   
00085             On exit, E may be multiplied by a constant factor chosen   
00086             to avoid over/underflow in computing the eigenvalues.   
00087 
00088     VL      (input) DOUBLE PRECISION   
00089     VU      (input) DOUBLE PRECISION   
00090             If RANGE='V', the lower and upper bounds of the interval to   
00091             be searched for eigenvalues. VL < VU.   
00092             Not referenced if RANGE = 'A' or 'I'.   
00093 
00094     IL      (input) INTEGER   
00095     IU      (input) INTEGER   
00096             If RANGE='I', the indices (in ascending order) of the   
00097             smallest and largest eigenvalues to be returned.   
00098             1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.   
00099             Not referenced if RANGE = 'A' or 'V'.   
00100 
00101     ABSTOL  (input) DOUBLE PRECISION   
00102             The absolute error tolerance for the eigenvalues.   
00103             An approximate eigenvalue is accepted as converged   
00104             when it is determined to lie in an interval [a,b]   
00105             of width less than or equal to   
00106 
00107                     ABSTOL + EPS *   max( |a|,|b| ) ,   
00108 
00109             where EPS is the machine precision.  If ABSTOL is less   
00110             than or equal to zero, then  EPS*|T|  will be used in   
00111             its place, where |T| is the 1-norm of the tridiagonal   
00112             matrix.   
00113 
00114             Eigenvalues will be computed most accurately when ABSTOL is   
00115             set to twice the underflow threshold 2*DLAMCH('S'), not zero.   
00116             If this routine returns with INFO>0, indicating that some   
00117             eigenvectors did not converge, try setting ABSTOL to   
00118             2*DLAMCH('S').   
00119 
00120             See "Computing Small Singular Values of Bidiagonal Matrices   
00121             with Guaranteed High Relative Accuracy," by Demmel and   
00122             Kahan, LAPACK Working Note #3.   
00123 
00124     M       (output) INTEGER   
00125             The total number of eigenvalues found.  0 <= M <= N.   
00126             If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.   
00127 
00128     W       (output) DOUBLE PRECISION array, dimension (N)   
00129             The first M elements contain the selected eigenvalues in   
00130             ascending order.   
00131 
00132     Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )   
00133             If JOBZ = 'V', then if INFO = 0, the first M columns of Z   
00134             contain the orthonormal eigenvectors of the matrix A   
00135             corresponding to the selected eigenvalues, with the i-th   
00136             column of Z holding the eigenvector associated with W(i).   
00137             If an eigenvector fails to converge (INFO > 0), then that   
00138             column of Z contains the latest approximation to the   
00139             eigenvector, and the index of the eigenvector is returned   
00140             in IFAIL.  If JOBZ = 'N', then Z is not referenced.   
00141             Note: the user must ensure that at least max(1,M) columns are   
00142             supplied in the array Z; if RANGE = 'V', the exact value of M   
00143             is not known in advance and an upper bound must be used.   
00144 
00145     LDZ     (input) INTEGER   
00146             The leading dimension of the array Z.  LDZ >= 1, and if   
00147             JOBZ = 'V', LDZ >= max(1,N).   
00148 
00149     WORK    (workspace) DOUBLE PRECISION array, dimension (5*N)   
00150 
00151     IWORK   (workspace) INTEGER array, dimension (5*N)   
00152 
00153     IFAIL   (output) INTEGER array, dimension (N)   
00154             If JOBZ = 'V', then if INFO = 0, the first M elements of   
00155             IFAIL are zero.  If INFO > 0, then IFAIL contains the   
00156             indices of the eigenvectors that failed to converge.   
00157             If JOBZ = 'N', then IFAIL is not referenced.   
00158 
00159     INFO    (output) INTEGER   
00160             = 0:  successful exit   
00161             < 0:  if INFO = -i, the i-th argument had an illegal value   
00162             > 0:  if INFO = i, then i eigenvectors failed to converge.   
00163                   Their indices are stored in array IFAIL.   
00164 
00165     =====================================================================   
00166 
00167 
00168        Test the input parameters.   
00169 
00170        Parameter adjustments */
00171     /* Table of constant values */
00172      integer c__1 = 1;
00173     
00174     /* System generated locals */
00175     integer z_dim1, z_offset, i__1, i__2;
00176     Treal d__1, d__2;
00177     /* Local variables */
00178      integer imax;
00179      Treal rmin, rmax, tnrm;
00180      integer itmp1, i__, j;
00181      Treal sigma;
00182      char order[1];
00183      logical wantz;
00184      integer jj;
00185      logical alleig, indeig;
00186      integer iscale, indibl;
00187      logical valeig;
00188      Treal safmin;
00189      Treal bignum;
00190      integer indisp;
00191      integer indiwo;
00192      integer indwrk;
00193      integer nsplit;
00194      Treal smlnum, eps, vll, vuu, tmp1;
00195 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
00196 
00197 
00198     --d__;
00199     --e;
00200     --w;
00201     z_dim1 = *ldz;
00202     z_offset = 1 + z_dim1 * 1;
00203     z__ -= z_offset;
00204     --work;
00205     --iwork;
00206     --ifail;
00207 
00208     /* Function Body */
00209     wantz = template_blas_lsame(jobz, "V");
00210     alleig = template_blas_lsame(range, "A");
00211     valeig = template_blas_lsame(range, "V");
00212     indeig = template_blas_lsame(range, "I");
00213 
00214     *info = 0;
00215     if (! (wantz || template_blas_lsame(jobz, "N"))) {
00216         *info = -1;
00217     } else if (! (alleig || valeig || indeig)) {
00218         *info = -2;
00219     } else if (*n < 0) {
00220         *info = -3;
00221     } else {
00222         if (valeig) {
00223             if (*n > 0 && *vu <= *vl) {
00224                 *info = -7;
00225             }
00226         } else if (indeig) {
00227             if (*il < 1 || *il > maxMACRO(1,*n)) {
00228                 *info = -8;
00229             } else if (*iu < minMACRO(*n,*il) || *iu > *n) {
00230                 *info = -9;
00231             }
00232         }
00233     }
00234     if (*info == 0) {
00235         if (*ldz < 1 || wantz && *ldz < *n) {
00236             *info = -14;
00237         }
00238     }
00239 
00240     if (*info != 0) {
00241         i__1 = -(*info);
00242         template_blas_erbla("STEVX ", &i__1);
00243         return 0;
00244     }
00245 
00246 /*     Quick return if possible */
00247 
00248     *m = 0;
00249     if (*n == 0) {
00250         return 0;
00251     }
00252 
00253     if (*n == 1) {
00254         if (alleig || indeig) {
00255             *m = 1;
00256             w[1] = d__[1];
00257         } else {
00258             if (*vl < d__[1] && *vu >= d__[1]) {
00259                 *m = 1;
00260                 w[1] = d__[1];
00261             }
00262         }
00263         if (wantz) {
00264             z___ref(1, 1) = 1.;
00265         }
00266         return 0;
00267     }
00268 
00269 /*     Get machine constants. */
00270 
00271     safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00272     eps = template_lapack_lamch("Precision", (Treal)0);
00273     smlnum = safmin / eps;
00274     bignum = 1. / smlnum;
00275     rmin = template_blas_sqrt(smlnum);
00276 /* Computing MIN */
00277     d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
00278     rmax = minMACRO(d__1,d__2);
00279 
00280 /*     Scale matrix to allowable range, if necessary. */
00281 
00282     iscale = 0;
00283     if (valeig) {
00284         vll = *vl;
00285         vuu = *vu;
00286     } else {
00287         vll = 0.;
00288         vuu = 0.;
00289     }
00290     tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
00291     if (tnrm > 0. && tnrm < rmin) {
00292         iscale = 1;
00293         sigma = rmin / tnrm;
00294     } else if (tnrm > rmax) {
00295         iscale = 1;
00296         sigma = rmax / tnrm;
00297     }
00298     if (iscale == 1) {
00299         template_blas_scal(n, &sigma, &d__[1], &c__1);
00300         i__1 = *n - 1;
00301         template_blas_scal(&i__1, &sigma, &e[1], &c__1);
00302         if (valeig) {
00303             vll = *vl * sigma;
00304             vuu = *vu * sigma;
00305         }
00306     }
00307 
00308 /*     If all eigenvalues are desired and ABSTOL is less than zero, then   
00309        call DSTERF or SSTEQR.  If this fails for some eigenvalue, then   
00310        try DSTEBZ. */
00311 
00312     if ((alleig || indeig && *il == 1 && *iu == *n) && *abstol <= 0.) {
00313         template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1);
00314         i__1 = *n - 1;
00315         template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1);
00316         indwrk = *n + 1;
00317         if (! wantz) {
00318             template_lapack_sterf(n, &w[1], &work[1], info);
00319         } else {
00320             template_lapack_steqr("I", n, &w[1], &work[1], &z__[z_offset], ldz, &work[
00321                     indwrk], info);
00322             if (*info == 0) {
00323                 i__1 = *n;
00324                 for (i__ = 1; i__ <= i__1; ++i__) {
00325                     ifail[i__] = 0;
00326 /* L10: */
00327                 }
00328             }
00329         }
00330         if (*info == 0) {
00331             *m = *n;
00332             goto L20;
00333         }
00334         *info = 0;
00335     }
00336 
00337 /*     Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. */
00338 
00339     if (wantz) {
00340         *(unsigned char *)order = 'B';
00341     } else {
00342         *(unsigned char *)order = 'E';
00343     }
00344     indwrk = 1;
00345     indibl = 1;
00346     indisp = indibl + *n;
00347     indiwo = indisp + *n;
00348     template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, &
00349             nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[indwrk], &
00350             iwork[indiwo], info);
00351 
00352     if (wantz) {
00353         template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], &
00354                 z__[z_offset], ldz, &work[indwrk], &iwork[indiwo], &ifail[1], 
00355                 info);
00356     }
00357 
00358 /*     If matrix was scaled, then rescale eigenvalues appropriately. */
00359 
00360 L20:
00361     if (iscale == 1) {
00362         if (*info == 0) {
00363             imax = *m;
00364         } else {
00365             imax = *info - 1;
00366         }
00367         d__1 = 1. / sigma;
00368         template_blas_scal(&imax, &d__1, &w[1], &c__1);
00369     }
00370 
00371 /*     If eigenvalues are not in order, then sort them, along with   
00372        eigenvectors. */
00373 
00374     if (wantz) {
00375         i__1 = *m - 1;
00376         for (j = 1; j <= i__1; ++j) {
00377             i__ = 0;
00378             tmp1 = w[j];
00379             i__2 = *m;
00380             for (jj = j + 1; jj <= i__2; ++jj) {
00381                 if (w[jj] < tmp1) {
00382                     i__ = jj;
00383                     tmp1 = w[jj];
00384                 }
00385 /* L30: */
00386             }
00387 
00388             if (i__ != 0) {
00389                 itmp1 = iwork[indibl + i__ - 1];
00390                 w[i__] = w[j];
00391                 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
00392                 w[j] = tmp1;
00393                 iwork[indibl + j - 1] = itmp1;
00394                 template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, j), &c__1);
00395                 if (*info != 0) {
00396                     itmp1 = ifail[i__];
00397                     ifail[i__] = ifail[j];
00398                     ifail[j] = itmp1;
00399                 }
00400             }
00401 /* L40: */
00402         }
00403     }
00404 
00405     return 0;
00406 
00407 /*     End of DSTEVX */
00408 
00409 } /* dstevx_ */
00410 
00411 #undef z___ref
00412 
00413 
00414 #endif

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7