ergo
|
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