00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_STEVR_HEADER 00038 #define TEMPLATE_LAPACK_STEVR_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_stevr(const char *jobz, const char *range, const integer *n, 00042 Treal * d__, Treal *e, const Treal *vl, 00043 const Treal *vu, const integer *il, 00044 const integer *iu, const Treal *abstol, 00045 integer *m, Treal *w, 00046 Treal *z__, const integer *ldz, integer *isuppz, 00047 Treal *work, 00048 integer *lwork, integer *iwork, integer *liwork, 00049 integer *info) 00050 { 00051 /* System generated locals */ 00052 integer z_dim1, z_offset, i__1, i__2; 00053 Treal d__1, d__2; 00054 00055 /* Builtin functions */ 00056 00057 /* Local variables */ 00058 integer i__, j, jj; 00059 Treal eps, vll, vuu, tmp1; 00060 integer imax; 00061 Treal rmin, rmax; 00062 logical test; 00063 Treal tnrm; 00064 integer itmp1; 00065 Treal sigma; 00066 char order[1]; 00067 integer lwmin; 00068 logical wantz; 00069 logical alleig, indeig; 00070 integer iscale, ieeeok, indibl, indifl; 00071 logical valeig; 00072 Treal safmin; 00073 Treal bignum; 00074 integer indisp; 00075 integer indiwo; 00076 integer liwmin; 00077 logical tryrac; 00078 integer nsplit; 00079 Treal smlnum; 00080 logical lquery; 00081 00082 00083 /* -- LAPACK driver routine (version 3.2) -- */ 00084 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00085 /* November 2006 */ 00086 00087 /* .. Scalar Arguments .. */ 00088 /* .. */ 00089 /* .. Array Arguments .. */ 00090 /* .. */ 00091 00092 /* Purpose */ 00093 /* ======= */ 00094 00095 /* DSTEVR computes selected eigenvalues and, optionally, eigenvectors */ 00096 /* of a real symmetric tridiagonal matrix T. Eigenvalues and */ 00097 /* eigenvectors can be selected by specifying either a range of values */ 00098 /* or a range of indices for the desired eigenvalues. */ 00099 00100 /* Whenever possible, DSTEVR calls DSTEMR to compute the */ 00101 /* eigenspectrum using Relatively Robust Representations. DSTEMR */ 00102 /* computes eigenvalues by the dqds algorithm, while orthogonal */ 00103 /* eigenvectors are computed from various "good" L D L^T representations */ 00104 /* (also known as Relatively Robust Representations). Gram-Schmidt */ 00105 /* orthogonalization is avoided as far as possible. More specifically, */ 00106 /* the various steps of the algorithm are as follows. For the i-th */ 00107 /* unreduced block of T, */ 00108 /* (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T */ 00109 /* is a relatively robust representation, */ 00110 /* (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high */ 00111 /* relative accuracy by the dqds algorithm, */ 00112 /* (c) If there is a cluster of close eigenvalues, "choose" sigma_i */ 00113 /* close to the cluster, and go to step (a), */ 00114 /* (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, */ 00115 /* compute the corresponding eigenvector by forming a */ 00116 /* rank-revealing twisted factorization. */ 00117 /* The desired accuracy of the output can be specified by the input */ 00118 /* parameter ABSTOL. */ 00119 00120 /* For more details, see "A new O(n^2) algorithm for the symmetric */ 00121 /* tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, */ 00122 /* Computer Science Division Technical Report No. UCB//CSD-97-971, */ 00123 /* UC Berkeley, May 1997. */ 00124 00125 00126 /* Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested */ 00127 /* on machines which conform to the ieee-754 floating point standard. */ 00128 /* DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and */ 00129 /* when partial spectrum requests are made. */ 00130 00131 /* Normal execution of DSTEMR may create NaNs and infinities and */ 00132 /* hence may abort due to a floating point exception in environments */ 00133 /* which do not handle NaNs and infinities in the ieee standard default */ 00134 /* manner. */ 00135 00136 /* Arguments */ 00137 /* ========= */ 00138 00139 /* JOBZ (input) CHARACTER*1 */ 00140 /* = 'N': Compute eigenvalues only; */ 00141 /* = 'V': Compute eigenvalues and eigenvectors. */ 00142 00143 /* RANGE (input) CHARACTER*1 */ 00144 /* = 'A': all eigenvalues will be found. */ 00145 /* = 'V': all eigenvalues in the half-open interval (VL,VU] */ 00146 /* will be found. */ 00147 /* = 'I': the IL-th through IU-th eigenvalues will be found. */ 00148 /* ********* For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and */ 00149 /* ********* DSTEIN are called */ 00150 00151 /* N (input) INTEGER */ 00152 /* The order of the matrix. N >= 0. */ 00153 00154 /* D (input/output) DOUBLE PRECISION array, dimension (N) */ 00155 /* On entry, the n diagonal elements of the tridiagonal matrix */ 00156 /* A. */ 00157 /* On exit, D may be multiplied by a constant factor chosen */ 00158 /* to avoid over/underflow in computing the eigenvalues. */ 00159 00160 /* E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1)) */ 00161 /* On entry, the (n-1) subdiagonal elements of the tridiagonal */ 00162 /* matrix A in elements 1 to N-1 of E. */ 00163 /* On exit, E may be multiplied by a constant factor chosen */ 00164 /* to avoid over/underflow in computing the eigenvalues. */ 00165 00166 /* VL (input) DOUBLE PRECISION */ 00167 /* VU (input) DOUBLE PRECISION */ 00168 /* If RANGE='V', the lower and upper bounds of the interval to */ 00169 /* be searched for eigenvalues. VL < VU. */ 00170 /* Not referenced if RANGE = 'A' or 'I'. */ 00171 00172 /* IL (input) INTEGER */ 00173 /* IU (input) INTEGER */ 00174 /* If RANGE='I', the indices (in ascending order) of the */ 00175 /* smallest and largest eigenvalues to be returned. */ 00176 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */ 00177 /* Not referenced if RANGE = 'A' or 'V'. */ 00178 00179 /* ABSTOL (input) DOUBLE PRECISION */ 00180 /* The absolute error tolerance for the eigenvalues. */ 00181 /* An approximate eigenvalue is accepted as converged */ 00182 /* when it is determined to lie in an interval [a,b] */ 00183 /* of width less than or equal to */ 00184 00185 /* ABSTOL + EPS * max( |a|,|b| ) , */ 00186 00187 /* where EPS is the machine precision. If ABSTOL is less than */ 00188 /* or equal to zero, then EPS*|T| will be used in its place, */ 00189 /* where |T| is the 1-norm of the tridiagonal matrix obtained */ 00190 /* by reducing A to tridiagonal form. */ 00191 00192 /* See "Computing Small Singular Values of Bidiagonal Matrices */ 00193 /* with Guaranteed High Relative Accuracy," by Demmel and */ 00194 /* Kahan, LAPACK Working Note #3. */ 00195 00196 /* If high relative accuracy is important, set ABSTOL to */ 00197 /* DLAMCH( 'Safe minimum' ). Doing so will guarantee that */ 00198 /* eigenvalues are computed to high relative accuracy when */ 00199 /* possible in future releases. The current code does not */ 00200 /* make any guarantees about high relative accuracy, but */ 00201 /* future releases will. See J. Barlow and J. Demmel, */ 00202 /* "Computing Accurate Eigensystems of Scaled Diagonally */ 00203 /* Dominant Matrices", LAPACK Working Note #7, for a discussion */ 00204 /* of which matrices define their eigenvalues to high relative */ 00205 /* accuracy. */ 00206 00207 /* M (output) INTEGER */ 00208 /* The total number of eigenvalues found. 0 <= M <= N. */ 00209 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */ 00210 00211 /* W (output) DOUBLE PRECISION array, dimension (N) */ 00212 /* The first M elements contain the selected eigenvalues in */ 00213 /* ascending order. */ 00214 00215 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */ 00216 /* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */ 00217 /* contain the orthonormal eigenvectors of the matrix A */ 00218 /* corresponding to the selected eigenvalues, with the i-th */ 00219 /* column of Z holding the eigenvector associated with W(i). */ 00220 /* Note: the user must ensure that at least max(1,M) columns are */ 00221 /* supplied in the array Z; if RANGE = 'V', the exact value of M */ 00222 /* is not known in advance and an upper bound must be used. */ 00223 00224 /* LDZ (input) INTEGER */ 00225 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00226 /* JOBZ = 'V', LDZ >= max(1,N). */ 00227 00228 /* ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) */ 00229 /* The support of the eigenvectors in Z, i.e., the indices */ 00230 /* indicating the nonzero elements in Z. The i-th eigenvector */ 00231 /* is nonzero only in elements ISUPPZ( 2*i-1 ) through */ 00232 /* ISUPPZ( 2*i ). */ 00233 /* ********* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 */ 00234 00235 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ 00236 /* On exit, if INFO = 0, WORK(1) returns the optimal (and */ 00237 /* minimal) LWORK. */ 00238 00239 /* LWORK (input) INTEGER */ 00240 /* The dimension of the array WORK. LWORK >= max(1,20*N). */ 00241 00242 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00243 /* only calculates the optimal sizes of the WORK and IWORK */ 00244 /* arrays, returns these values as the first entries of the WORK */ 00245 /* and IWORK arrays, and no error message related to LWORK or */ 00246 /* LIWORK is issued by XERBLA. */ 00247 00248 /* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */ 00249 /* On exit, if INFO = 0, IWORK(1) returns the optimal (and */ 00250 /* minimal) LIWORK. */ 00251 00252 /* LIWORK (input) INTEGER */ 00253 /* The dimension of the array IWORK. LIWORK >= max(1,10*N). */ 00254 00255 /* If LIWORK = -1, then a workspace query is assumed; the */ 00256 /* routine only calculates the optimal sizes of the WORK and */ 00257 /* IWORK arrays, returns these values as the first entries of */ 00258 /* the WORK and IWORK arrays, and no error message related to */ 00259 /* LWORK or LIWORK is issued by XERBLA. */ 00260 00261 /* INFO (output) INTEGER */ 00262 /* = 0: successful exit */ 00263 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00264 /* > 0: Internal error */ 00265 00266 /* Further Details */ 00267 /* =============== */ 00268 00269 /* Based on contributions by */ 00270 /* Inderjit Dhillon, IBM Almaden, USA */ 00271 /* Osni Marques, LBNL/NERSC, USA */ 00272 /* Ken Stanley, Computer Science Division, University of */ 00273 /* California at Berkeley, USA */ 00274 00275 /* ===================================================================== */ 00276 00277 /* .. Parameters .. */ 00278 /* .. */ 00279 /* .. Local Scalars .. */ 00280 /* .. */ 00281 /* .. External Functions .. */ 00282 /* .. */ 00283 /* .. External Subroutines .. */ 00284 /* .. */ 00285 /* .. Intrinsic Functions .. */ 00286 /* .. */ 00287 /* .. Executable Statements .. */ 00288 00289 00290 /* Test the input parameters. */ 00291 00292 /* Parameter adjustments */ 00293 /* Table of constant values */ 00294 00295 integer c__10 = 10; 00296 integer c__1 = 1; 00297 integer c__2 = 2; 00298 integer c__3 = 3; 00299 integer c__4 = 4; 00300 00301 --d__; 00302 --e; 00303 --w; 00304 z_dim1 = *ldz; 00305 z_offset = 1 + z_dim1; 00306 z__ -= z_offset; 00307 --isuppz; 00308 --work; 00309 --iwork; 00310 00311 /* Function Body */ 00312 ieeeok = template_lapack_ilaenv(&c__10, "DSTEVR", "N", &c__1, &c__2, &c__3, &c__4, (ftnlen)6, (ftnlen)1); 00313 00314 wantz = template_blas_lsame(jobz, "V"); 00315 alleig = template_blas_lsame(range, "A"); 00316 valeig = template_blas_lsame(range, "V"); 00317 indeig = template_blas_lsame(range, "I"); 00318 00319 lquery = *lwork == -1 || *liwork == -1; 00320 /* Computing MAX */ 00321 i__1 = 1, i__2 = *n * 20; 00322 lwmin = maxMACRO(i__1,i__2); 00323 /* Computing MAX */ 00324 i__1 = 1, i__2 = *n * 10; 00325 liwmin = maxMACRO(i__1,i__2); 00326 00327 00328 *info = 0; 00329 if (! (wantz || template_blas_lsame(jobz, "N"))) { 00330 *info = -1; 00331 } else if (! (alleig || valeig || indeig)) { 00332 *info = -2; 00333 } else if (*n < 0) { 00334 *info = -3; 00335 } else { 00336 if (valeig) { 00337 if (*n > 0 && *vu <= *vl) { 00338 *info = -7; 00339 } 00340 } else if (indeig) { 00341 if (*il < 1 || *il > maxMACRO(1,*n)) { 00342 *info = -8; 00343 } else if (*iu < minMACRO(*n,*il) || *iu > *n) { 00344 *info = -9; 00345 } 00346 } 00347 } 00348 if (*info == 0) { 00349 if (*ldz < 1 || ( wantz && *ldz < *n ) ) { 00350 *info = -14; 00351 } 00352 } 00353 00354 if (*info == 0) { 00355 work[1] = (Treal) lwmin; 00356 iwork[1] = liwmin; 00357 00358 if (*lwork < lwmin && ! lquery) { 00359 *info = -17; 00360 } else if (*liwork < liwmin && ! lquery) { 00361 *info = -19; 00362 } 00363 } 00364 00365 if (*info != 0) { 00366 i__1 = -(*info); 00367 template_blas_erbla("STEVR", &i__1); 00368 return 0; 00369 } else if (lquery) { 00370 return 0; 00371 } 00372 00373 /* Quick return if possible */ 00374 00375 *m = 0; 00376 if (*n == 0) { 00377 return 0; 00378 } 00379 00380 if (*n == 1) { 00381 if (alleig || indeig) { 00382 *m = 1; 00383 w[1] = d__[1]; 00384 } else { 00385 if (*vl < d__[1] && *vu >= d__[1]) { 00386 *m = 1; 00387 w[1] = d__[1]; 00388 } 00389 } 00390 if (wantz) { 00391 z__[z_dim1 + 1] = 1.; 00392 } 00393 return 0; 00394 } 00395 00396 /* Get machine constants. */ 00397 00398 safmin = template_lapack_lamch("Safe minimum", (Treal)0); 00399 eps = template_lapack_lamch("Precision", (Treal)0); 00400 smlnum = safmin / eps; 00401 bignum = 1. / smlnum; 00402 rmin = template_blas_sqrt(smlnum); 00403 /* Computing MIN */ 00404 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin)); 00405 rmax = minMACRO(d__1,d__2); 00406 00407 00408 /* Scale matrix to allowable range, if necessary. */ 00409 00410 iscale = 0; 00411 vll = *vl; 00412 vuu = *vu; 00413 00414 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]); 00415 if (tnrm > 0. && tnrm < rmin) { 00416 iscale = 1; 00417 sigma = rmin / tnrm; 00418 } else if (tnrm > rmax) { 00419 iscale = 1; 00420 sigma = rmax / tnrm; 00421 } 00422 if (iscale == 1) { 00423 template_blas_scal(n, &sigma, &d__[1], &c__1); 00424 i__1 = *n - 1; 00425 template_blas_scal(&i__1, &sigma, &e[1], &c__1); 00426 if (valeig) { 00427 vll = *vl * sigma; 00428 vuu = *vu * sigma; 00429 } 00430 } 00431 /* Initialize indices into workspaces. Note: These indices are used only */ 00432 /* if DSTERF or DSTEMR fail. */ 00433 /* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and */ 00434 /* stores the block indices of each of the M<=N eigenvalues. */ 00435 indibl = 1; 00436 /* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and */ 00437 /* stores the starting and finishing indices of each block. */ 00438 indisp = indibl + *n; 00439 /* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors */ 00440 /* that corresponding to eigenvectors that fail to converge in */ 00441 /* DSTEIN. This information is discarded; if any fail, the driver */ 00442 /* returns INFO > 0. */ 00443 indifl = indisp + *n; 00444 /* INDIWO is the offset of the remaining integer workspace. */ 00445 indiwo = indisp + *n; 00446 00447 /* If all eigenvalues are desired, then */ 00448 /* call DSTERF or DSTEMR. If this fails for some eigenvalue, then */ 00449 /* try DSTEBZ. */ 00450 00451 00452 test = FALSE_; 00453 if (indeig) { 00454 if (*il == 1 && *iu == *n) { 00455 test = TRUE_; 00456 } 00457 } 00458 if ((alleig || test) && ieeeok == 1) { 00459 i__1 = *n - 1; 00460 template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1); 00461 if (! wantz) { 00462 template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1); 00463 template_lapack_sterf(n, &w[1], &work[1], info); 00464 } else { 00465 template_blas_copy(n, &d__[1], &c__1, &work[*n + 1], &c__1); 00466 if (*abstol <= *n * 2. * eps) { 00467 tryrac = TRUE_; 00468 } else { 00469 tryrac = FALSE_; 00470 } 00471 i__1 = *lwork - (*n << 1); 00472 template_lapack_stemr(jobz, "A", n, &work[*n + 1], &work[1], vl, vu, il, iu, m, 00473 &w[1], &z__[z_offset], ldz, n, &isuppz[1], &tryrac, &work[ 00474 (*n << 1) + 1], &i__1, &iwork[1], liwork, info); 00475 00476 } 00477 if (*info == 0) { 00478 *m = *n; 00479 goto L10; 00480 } 00481 *info = 0; 00482 } 00483 00484 /* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN. */ 00485 00486 if (wantz) { 00487 *(unsigned char *)order = 'B'; 00488 } else { 00489 *(unsigned char *)order = 'E'; 00490 } 00491 template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, & 00492 nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[1], &iwork[ 00493 indiwo], info); 00494 00495 if (wantz) { 00496 template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], & 00497 z__[z_offset], ldz, &work[1], &iwork[indiwo], &iwork[indifl], 00498 info); 00499 } 00500 00501 /* If matrix was scaled, then rescale eigenvalues appropriately. */ 00502 00503 L10: 00504 if (iscale == 1) { 00505 if (*info == 0) { 00506 imax = *m; 00507 } else { 00508 imax = *info - 1; 00509 } 00510 d__1 = 1. / sigma; 00511 template_blas_scal(&imax, &d__1, &w[1], &c__1); 00512 } 00513 00514 /* If eigenvalues are not in order, then sort them, along with */ 00515 /* eigenvectors. */ 00516 00517 if (wantz) { 00518 i__1 = *m - 1; 00519 for (j = 1; j <= i__1; ++j) { 00520 i__ = 0; 00521 tmp1 = w[j]; 00522 i__2 = *m; 00523 for (jj = j + 1; jj <= i__2; ++jj) { 00524 if (w[jj] < tmp1) { 00525 i__ = jj; 00526 tmp1 = w[jj]; 00527 } 00528 /* L20: */ 00529 } 00530 00531 if (i__ != 0) { 00532 itmp1 = iwork[i__]; 00533 w[i__] = w[j]; 00534 iwork[i__] = iwork[j]; 00535 w[j] = tmp1; 00536 iwork[j] = itmp1; 00537 template_blas_swap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1], 00538 &c__1); 00539 } 00540 /* L30: */ 00541 } 00542 } 00543 00544 /* Causes problems with tests 19 & 20: */ 00545 /* IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002 */ 00546 00547 00548 work[1] = (Treal) lwmin; 00549 iwork[1] = liwmin; 00550 return 0; 00551 00552 /* End of DSTEVR */ 00553 00554 } /* dstevr_ */ 00555 00556 #endif 00557