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_STEMR_HEADER 00038 #define TEMPLATE_LAPACK_STEMR_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_stemr(const char *jobz, const char *range, const integer *n, Treal * 00042 d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il, 00043 const integer *iu, integer *m, Treal *w, Treal *z__, const integer *ldz, 00044 const integer *nzc, integer *isuppz, logical *tryrac, Treal *work, 00045 integer *lwork, integer *iwork, integer *liwork, integer *info) 00046 { 00047 /* System generated locals */ 00048 integer z_dim1, z_offset, i__1, i__2; 00049 Treal d__1, d__2; 00050 00051 /* Builtin functions */ 00052 00053 /* Local variables */ 00054 integer i__, j; 00055 Treal r1, r2; 00056 integer jj; 00057 Treal cs = 0; 00058 integer in; 00059 Treal sn = 0, wl, wu; 00060 integer iil, iiu; 00061 Treal eps, tmp; 00062 integer indd, iend, jblk, wend; 00063 Treal rmin, rmax; 00064 integer itmp; 00065 Treal tnrm; 00066 integer inde2, itmp2; 00067 Treal rtol1, rtol2; 00068 Treal scale; 00069 integer indgp; 00070 integer iinfo, iindw, ilast; 00071 integer lwmin; 00072 logical wantz; 00073 logical alleig; 00074 integer ibegin; 00075 logical indeig; 00076 integer iindbl; 00077 logical valeig; 00078 integer wbegin; 00079 Treal safmin; 00080 Treal bignum; 00081 integer inderr, iindwk, indgrs, offset; 00082 Treal thresh; 00083 integer iinspl, ifirst, indwrk, liwmin, nzcmin; 00084 Treal pivmin; 00085 integer nsplit; 00086 Treal smlnum; 00087 logical lquery, zquery; 00088 00089 00090 /* -- LAPACK computational routine (version 3.2) -- */ 00091 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00092 /* November 2006 */ 00093 00094 /* .. Scalar Arguments .. */ 00095 /* .. */ 00096 /* .. Array Arguments .. */ 00097 /* .. */ 00098 00099 /* Purpose */ 00100 /* ======= */ 00101 00102 /* DSTEMR computes selected eigenvalues and, optionally, eigenvectors */ 00103 /* of a real symmetric tridiagonal matrix T. Any such unreduced matrix has */ 00104 /* a well defined set of pairwise different real eigenvalues, the corresponding */ 00105 /* real eigenvectors are pairwise orthogonal. */ 00106 00107 /* The spectrum may be computed either completely or partially by specifying */ 00108 /* either an interval (VL,VU] or a range of indices IL:IU for the desired */ 00109 /* eigenvalues. */ 00110 00111 /* Depending on the number of desired eigenvalues, these are computed either */ 00112 /* by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are */ 00113 /* computed by the use of various suitable L D L^T factorizations near clusters */ 00114 /* of close eigenvalues (referred to as RRRs, Relatively Robust */ 00115 /* Representations). An informal sketch of the algorithm follows. */ 00116 00117 /* For each unreduced block (submatrix) of T, */ 00118 /* (a) Compute T - sigma I = L D L^T, so that L and D */ 00119 /* define all the wanted eigenvalues to high relative accuracy. */ 00120 /* This means that small relative changes in the entries of D and L */ 00121 /* cause only small relative changes in the eigenvalues and */ 00122 /* eigenvectors. The standard (unfactored) representation of the */ 00123 /* tridiagonal matrix T does not have this property in general. */ 00124 /* (b) Compute the eigenvalues to suitable accuracy. */ 00125 /* If the eigenvectors are desired, the algorithm attains full */ 00126 /* accuracy of the computed eigenvalues only right before */ 00127 /* the corresponding vectors have to be computed, see steps c) and d). */ 00128 /* (c) For each cluster of close eigenvalues, select a new */ 00129 /* shift close to the cluster, find a new factorization, and refine */ 00130 /* the shifted eigenvalues to suitable accuracy. */ 00131 /* (d) For each eigenvalue with a large enough relative separation compute */ 00132 /* the corresponding eigenvector by forming a rank revealing twisted */ 00133 /* factorization. Go back to (c) for any clusters that remain. */ 00134 00135 /* For more details, see: */ 00136 /* - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations */ 00137 /* to compute orthogonal eigenvectors of symmetric tridiagonal matrices," */ 00138 /* Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. */ 00139 /* - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and */ 00140 /* Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, */ 00141 /* 2004. Also LAPACK Working Note 154. */ 00142 /* - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric */ 00143 /* tridiagonal eigenvalue/eigenvector problem", */ 00144 /* Computer Science Division Technical Report No. UCB/CSD-97-971, */ 00145 /* UC Berkeley, May 1997. */ 00146 00147 /* Notes: */ 00148 /* 1.DSTEMR works only on machines which follow IEEE-754 */ 00149 /* floating-point standard in their handling of infinities and NaNs. */ 00150 /* This permits the use of efficient inner loops avoiding a check for */ 00151 /* zero divisors. */ 00152 00153 /* Arguments */ 00154 /* ========= */ 00155 00156 /* JOBZ (input) CHARACTER*1 */ 00157 /* = 'N': Compute eigenvalues only; */ 00158 /* = 'V': Compute eigenvalues and eigenvectors. */ 00159 00160 /* RANGE (input) CHARACTER*1 */ 00161 /* = 'A': all eigenvalues will be found. */ 00162 /* = 'V': all eigenvalues in the half-open interval (VL,VU] */ 00163 /* will be found. */ 00164 /* = 'I': the IL-th through IU-th eigenvalues will be found. */ 00165 00166 /* N (input) INTEGER */ 00167 /* The order of the matrix. N >= 0. */ 00168 00169 /* D (input/output) DOUBLE PRECISION array, dimension (N) */ 00170 /* On entry, the N diagonal elements of the tridiagonal matrix */ 00171 /* T. On exit, D is overwritten. */ 00172 00173 /* E (input/output) DOUBLE PRECISION array, dimension (N) */ 00174 /* On entry, the (N-1) subdiagonal elements of the tridiagonal */ 00175 /* matrix T in elements 1 to N-1 of E. E(N) need not be set on */ 00176 /* input, but is used internally as workspace. */ 00177 /* On exit, E is overwritten. */ 00178 00179 /* VL (input) DOUBLE PRECISION */ 00180 /* VU (input) DOUBLE PRECISION */ 00181 /* If RANGE='V', the lower and upper bounds of the interval to */ 00182 /* be searched for eigenvalues. VL < VU. */ 00183 /* Not referenced if RANGE = 'A' or 'I'. */ 00184 00185 /* IL (input) INTEGER */ 00186 /* IU (input) INTEGER */ 00187 /* If RANGE='I', the indices (in ascending order) of the */ 00188 /* smallest and largest eigenvalues to be returned. */ 00189 /* 1 <= IL <= IU <= N, if N > 0. */ 00190 /* Not referenced if RANGE = 'A' or 'V'. */ 00191 00192 /* M (output) INTEGER */ 00193 /* The total number of eigenvalues found. 0 <= M <= N. */ 00194 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */ 00195 00196 /* W (output) DOUBLE PRECISION array, dimension (N) */ 00197 /* The first M elements contain the selected eigenvalues in */ 00198 /* ascending order. */ 00199 00200 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */ 00201 /* If JOBZ = 'V', and if INFO = 0, then the first M columns of Z */ 00202 /* contain the orthonormal eigenvectors of the matrix T */ 00203 /* corresponding to the selected eigenvalues, with the i-th */ 00204 /* column of Z holding the eigenvector associated with W(i). */ 00205 /* If JOBZ = 'N', then Z is not referenced. */ 00206 /* Note: the user must ensure that at least max(1,M) columns are */ 00207 /* supplied in the array Z; if RANGE = 'V', the exact value of M */ 00208 /* is not known in advance and can be computed with a workspace */ 00209 /* query by setting NZC = -1, see below. */ 00210 00211 /* LDZ (input) INTEGER */ 00212 /* The leading dimension of the array Z. LDZ >= 1, and if */ 00213 /* JOBZ = 'V', then LDZ >= max(1,N). */ 00214 00215 /* NZC (input) INTEGER */ 00216 /* The number of eigenvectors to be held in the array Z. */ 00217 /* If RANGE = 'A', then NZC >= max(1,N). */ 00218 /* If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. */ 00219 /* If RANGE = 'I', then NZC >= IU-IL+1. */ 00220 /* If NZC = -1, then a workspace query is assumed; the */ 00221 /* routine calculates the number of columns of the array Z that */ 00222 /* are needed to hold the eigenvectors. */ 00223 /* This value is returned as the first entry of the Z array, and */ 00224 /* no error message related to NZC is issued by XERBLA. */ 00225 00226 /* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */ 00227 /* The support of the eigenvectors in Z, i.e., the indices */ 00228 /* indicating the nonzero elements in Z. The i-th computed eigenvector */ 00229 /* is nonzero only in elements ISUPPZ( 2*i-1 ) through */ 00230 /* ISUPPZ( 2*i ). This is relevant in the case when the matrix */ 00231 /* is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. */ 00232 00233 /* TRYRAC (input/output) LOGICAL */ 00234 /* If TRYRAC.EQ..TRUE., indicates that the code should check whether */ 00235 /* the tridiagonal matrix defines its eigenvalues to high relative */ 00236 /* accuracy. If so, the code uses relative-accuracy preserving */ 00237 /* algorithms that might be (a bit) slower depending on the matrix. */ 00238 /* If the matrix does not define its eigenvalues to high relative */ 00239 /* accuracy, the code can uses possibly faster algorithms. */ 00240 /* If TRYRAC.EQ..FALSE., the code is not required to guarantee */ 00241 /* relatively accurate eigenvalues and can use the fastest possible */ 00242 /* techniques. */ 00243 /* On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix */ 00244 /* does not define its eigenvalues to high relative accuracy. */ 00245 00246 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */ 00247 /* On exit, if INFO = 0, WORK(1) returns the optimal */ 00248 /* (and minimal) LWORK. */ 00249 00250 /* LWORK (input) INTEGER */ 00251 /* The dimension of the array WORK. LWORK >= max(1,18*N) */ 00252 /* if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. */ 00253 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00254 /* only calculates the optimal size of the WORK array, returns */ 00255 /* this value as the first entry of the WORK array, and no error */ 00256 /* message related to LWORK is issued by XERBLA. */ 00257 00258 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */ 00259 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */ 00260 00261 /* LIWORK (input) INTEGER */ 00262 /* The dimension of the array IWORK. LIWORK >= max(1,10*N) */ 00263 /* if the eigenvectors are desired, and LIWORK >= max(1,8*N) */ 00264 /* if only the eigenvalues are to be computed. */ 00265 /* If LIWORK = -1, then a workspace query is assumed; the */ 00266 /* routine only calculates the optimal size of the IWORK array, */ 00267 /* returns this value as the first entry of the IWORK array, and */ 00268 /* no error message related to LIWORK is issued by XERBLA. */ 00269 00270 /* INFO (output) INTEGER */ 00271 /* On exit, INFO */ 00272 /* = 0: successful exit */ 00273 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00274 /* > 0: if INFO = 1X, internal error in DLARRE, */ 00275 /* if INFO = 2X, internal error in DLARRV. */ 00276 /* Here, the digit X = ABS( IINFO ) < 10, where IINFO is */ 00277 /* the nonzero error code returned by DLARRE or */ 00278 /* DLARRV, respectively. */ 00279 00280 00281 /* Further Details */ 00282 /* =============== */ 00283 00284 /* Based on contributions by */ 00285 /* Beresford Parlett, University of California, Berkeley, USA */ 00286 /* Jim Demmel, University of California, Berkeley, USA */ 00287 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00288 /* Osni Marques, LBNL/NERSC, USA */ 00289 /* Christof Voemel, University of California, Berkeley, USA */ 00290 00291 /* ===================================================================== */ 00292 00293 /* .. Parameters .. */ 00294 /* .. */ 00295 /* .. Local Scalars .. */ 00296 /* .. */ 00297 /* .. */ 00298 /* .. External Functions .. */ 00299 /* .. */ 00300 /* .. External Subroutines .. */ 00301 /* .. */ 00302 /* .. Intrinsic Functions .. */ 00303 /* .. */ 00304 /* .. Executable Statements .. */ 00305 00306 /* Test the input parameters. */ 00307 00308 /* Parameter adjustments */ 00309 /* Table of constant values */ 00310 integer c__1 = 1; 00311 Treal c_b18 = .001; 00312 00313 --d__; 00314 --e; 00315 --w; 00316 z_dim1 = *ldz; 00317 z_offset = 1 + z_dim1; 00318 z__ -= z_offset; 00319 --isuppz; 00320 --work; 00321 --iwork; 00322 00323 /* Function Body */ 00324 wantz = template_blas_lsame(jobz, "V"); 00325 alleig = template_blas_lsame(range, "A"); 00326 valeig = template_blas_lsame(range, "V"); 00327 indeig = template_blas_lsame(range, "I"); 00328 00329 lquery = *lwork == -1 || *liwork == -1; 00330 zquery = *nzc == -1; 00331 /* DSTEMR needs WORK of size 6*N, IWORK of size 3*N. */ 00332 /* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N. */ 00333 /* Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N. */ 00334 if (wantz) { 00335 lwmin = *n * 18; 00336 liwmin = *n * 10; 00337 } else { 00338 /* need less workspace if only the eigenvalues are wanted */ 00339 lwmin = *n * 12; 00340 liwmin = *n << 3; 00341 } 00342 wl = 0.; 00343 wu = 0.; 00344 iil = 0; 00345 iiu = 0; 00346 if (valeig) { 00347 /* We do not reference VL, VU in the cases RANGE = 'I','A' */ 00348 /* The interval (WL, WU] contains all the wanted eigenvalues. */ 00349 /* It is either given by the user or computed in DLARRE. */ 00350 wl = *vl; 00351 wu = *vu; 00352 } else if (indeig) { 00353 /* We do not reference IL, IU in the cases RANGE = 'V','A' */ 00354 iil = *il; 00355 iiu = *iu; 00356 } 00357 00358 *info = 0; 00359 if (! (wantz || template_blas_lsame(jobz, "N"))) { 00360 *info = -1; 00361 } else if (! (alleig || valeig || indeig)) { 00362 *info = -2; 00363 } else if (*n < 0) { 00364 *info = -3; 00365 } else if (valeig && *n > 0 && wu <= wl) { 00366 *info = -7; 00367 } else if (indeig && (iil < 1 || iil > *n)) { 00368 *info = -8; 00369 } else if (indeig && (iiu < iil || iiu > *n)) { 00370 *info = -9; 00371 } else if (*ldz < 1 || ( wantz && *ldz < *n ) ) { 00372 *info = -13; 00373 } else if (*lwork < lwmin && ! lquery) { 00374 *info = -17; 00375 } else if (*liwork < liwmin && ! lquery) { 00376 *info = -19; 00377 } 00378 00379 /* Get machine constants. */ 00380 00381 safmin = template_lapack_lamch("Safe minimum", (Treal)0); 00382 eps = template_lapack_lamch("Precision", (Treal)0); 00383 smlnum = safmin / eps; 00384 bignum = 1. / smlnum; 00385 rmin = template_blas_sqrt(smlnum); 00386 /* Computing MIN */ 00387 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin)); 00388 rmax = minMACRO(d__1,d__2); 00389 00390 if (*info == 0) { 00391 work[1] = (Treal) lwmin; 00392 iwork[1] = liwmin; 00393 00394 if (wantz && alleig) { 00395 nzcmin = *n; 00396 } else if (wantz && valeig) { 00397 template_lapack_larrc("T", n, vl, vu, &d__[1], &e[1], &safmin, &nzcmin, &itmp, & 00398 itmp2, info); 00399 } else if (wantz && indeig) { 00400 nzcmin = iiu - iil + 1; 00401 } else { 00402 /* WANTZ .EQ. FALSE. */ 00403 nzcmin = 0; 00404 } 00405 if (zquery && *info == 0) { 00406 z__[z_dim1 + 1] = (Treal) nzcmin; 00407 } else if (*nzc < nzcmin && ! zquery) { 00408 *info = -14; 00409 } 00410 } 00411 if (*info != 0) { 00412 00413 i__1 = -(*info); 00414 template_blas_erbla("DSTEMR", &i__1); 00415 00416 return 0; 00417 } else if (lquery || zquery) { 00418 return 0; 00419 } 00420 00421 /* Handle N = 0, 1, and 2 cases immediately */ 00422 00423 *m = 0; 00424 if (*n == 0) { 00425 return 0; 00426 } 00427 00428 if (*n == 1) { 00429 if (alleig || indeig) { 00430 *m = 1; 00431 w[1] = d__[1]; 00432 } else { 00433 if (wl < d__[1] && wu >= d__[1]) { 00434 *m = 1; 00435 w[1] = d__[1]; 00436 } 00437 } 00438 if (wantz && ! zquery) { 00439 z__[z_dim1 + 1] = 1.; 00440 isuppz[1] = 1; 00441 isuppz[2] = 1; 00442 } 00443 return 0; 00444 } 00445 00446 if (*n == 2) { 00447 if (! wantz) { 00448 template_lapack_lae2(&d__[1], &e[1], &d__[2], &r1, &r2); 00449 } else if (wantz && ! zquery) { 00450 template_lapack_laev2(&d__[1], &e[1], &d__[2], &r1, &r2, &cs, &sn); 00451 } 00452 if (alleig || ( valeig && r2 > wl && r2 <= wu ) || ( indeig && iil == 1 ) ) { 00453 ++(*m); 00454 w[*m] = r2; 00455 if (wantz && ! zquery) { 00456 z__[*m * z_dim1 + 1] = -sn; 00457 z__[*m * z_dim1 + 2] = cs; 00458 /* Note: At most one of SN and CS can be zero. */ 00459 if (sn != 0.) { 00460 if (cs != 0.) { 00461 isuppz[(*m << 1) - 1] = 1; 00462 isuppz[(*m << 1) - 1] = 2; 00463 } else { 00464 isuppz[(*m << 1) - 1] = 1; 00465 isuppz[(*m << 1) - 1] = 1; 00466 } 00467 } else { 00468 isuppz[(*m << 1) - 1] = 2; 00469 isuppz[*m * 2] = 2; 00470 } 00471 } 00472 } 00473 if (alleig || ( valeig && r1 > wl && r1 <= wu ) || ( indeig && iiu == 2 ) ) { 00474 ++(*m); 00475 w[*m] = r1; 00476 if (wantz && ! zquery) { 00477 z__[*m * z_dim1 + 1] = cs; 00478 z__[*m * z_dim1 + 2] = sn; 00479 /* Note: At most one of SN and CS can be zero. */ 00480 if (sn != 0.) { 00481 if (cs != 0.) { 00482 isuppz[(*m << 1) - 1] = 1; 00483 isuppz[(*m << 1) - 1] = 2; 00484 } else { 00485 isuppz[(*m << 1) - 1] = 1; 00486 isuppz[(*m << 1) - 1] = 1; 00487 } 00488 } else { 00489 isuppz[(*m << 1) - 1] = 2; 00490 isuppz[*m * 2] = 2; 00491 } 00492 } 00493 } 00494 return 0; 00495 } 00496 /* Continue with general N */ 00497 indgrs = 1; 00498 inderr = (*n << 1) + 1; 00499 indgp = *n * 3 + 1; 00500 indd = (*n << 2) + 1; 00501 inde2 = *n * 5 + 1; 00502 indwrk = *n * 6 + 1; 00503 00504 iinspl = 1; 00505 iindbl = *n + 1; 00506 iindw = (*n << 1) + 1; 00507 iindwk = *n * 3 + 1; 00508 00509 /* Scale matrix to allowable range, if necessary. */ 00510 /* The allowable range is related to the PIVMIN parameter; see the */ 00511 /* comments in DLARRD. The preference for scaling small values */ 00512 /* up is heuristic; we expect users' matrices not to be close to the */ 00513 /* RMAX threshold. */ 00514 00515 scale = 1.; 00516 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]); 00517 if (tnrm > 0. && tnrm < rmin) { 00518 scale = rmin / tnrm; 00519 } else if (tnrm > rmax) { 00520 scale = rmax / tnrm; 00521 } 00522 if (scale != 1.) { 00523 template_blas_scal(n, &scale, &d__[1], &c__1); 00524 i__1 = *n - 1; 00525 template_blas_scal(&i__1, &scale, &e[1], &c__1); 00526 tnrm *= scale; 00527 if (valeig) { 00528 /* If eigenvalues in interval have to be found, */ 00529 /* scale (WL, WU] accordingly */ 00530 wl *= scale; 00531 wu *= scale; 00532 } 00533 } 00534 00535 /* Compute the desired eigenvalues of the tridiagonal after splitting */ 00536 /* into smaller subblocks if the corresponding off-diagonal elements */ 00537 /* are small */ 00538 /* THRESH is the splitting parameter for DLARRE */ 00539 /* A negative THRESH forces the old splitting criterion based on the */ 00540 /* size of the off-diagonal. A positive THRESH switches to splitting */ 00541 /* which preserves relative accuracy. */ 00542 00543 if (*tryrac) { 00544 /* Test whether the matrix warrants the more expensive relative approach. */ 00545 template_lapack_larrr(n, &d__[1], &e[1], &iinfo); 00546 } else { 00547 /* The user does not care about relative accurately eigenvalues */ 00548 iinfo = -1; 00549 } 00550 /* Set the splitting criterion */ 00551 if (iinfo == 0) { 00552 thresh = eps; 00553 } else { 00554 thresh = -eps; 00555 /* relative accuracy is desired but T does not guarantee it */ 00556 *tryrac = FALSE_; 00557 } 00558 00559 if (*tryrac) { 00560 /* Copy original diagonal, needed to guarantee relative accuracy */ 00561 template_blas_copy(n, &d__[1], &c__1, &work[indd], &c__1); 00562 } 00563 /* Store the squares of the offdiagonal values of T */ 00564 i__1 = *n - 1; 00565 for (j = 1; j <= i__1; ++j) { 00566 /* Computing 2nd power */ 00567 d__1 = e[j]; 00568 work[inde2 + j - 1] = d__1 * d__1; 00569 /* L5: */ 00570 } 00571 /* Set the tolerance parameters for bisection */ 00572 if (! wantz) { 00573 /* DLARRE computes the eigenvalues to full precision. */ 00574 rtol1 = eps * 4.; 00575 rtol2 = eps * 4.; 00576 } else { 00577 /* DLARRE computes the eigenvalues to less than full precision. */ 00578 /* DLARRV will refine the eigenvalue approximations, and we can */ 00579 /* need less accurate initial bisection in DLARRE. */ 00580 /* Note: these settings do only affect the subset case and DLARRE */ 00581 rtol1 = template_blas_sqrt(eps); 00582 /* Computing MAX */ 00583 d__1 = template_blas_sqrt(eps) * .005, d__2 = eps * 4.; 00584 rtol2 = maxMACRO(d__1,d__2); 00585 } 00586 template_lapack_larre(range, n, &wl, &wu, &iil, &iiu, &d__[1], &e[1], &work[inde2], & 00587 rtol1, &rtol2, &thresh, &nsplit, &iwork[iinspl], m, &w[1], &work[ 00588 inderr], &work[indgp], &iwork[iindbl], &iwork[iindw], &work[ 00589 indgrs], &pivmin, &work[indwrk], &iwork[iindwk], &iinfo); 00590 if (iinfo != 0) { 00591 *info = absMACRO(iinfo) + 10; 00592 return 0; 00593 } 00594 /* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired */ 00595 /* part of the spectrum. All desired eigenvalues are contained in */ 00596 /* (WL,WU] */ 00597 if (wantz) { 00598 00599 /* Compute the desired eigenvectors corresponding to the computed */ 00600 /* eigenvalues */ 00601 00602 template_lapack_larrv(n, &wl, &wu, &d__[1], &e[1], &pivmin, &iwork[iinspl], m, & 00603 c__1, m, &c_b18, &rtol1, &rtol2, &w[1], &work[inderr], &work[ 00604 indgp], &iwork[iindbl], &iwork[iindw], &work[indgrs], &z__[ 00605 z_offset], ldz, &isuppz[1], &work[indwrk], &iwork[iindwk], & 00606 iinfo); 00607 if (iinfo != 0) { 00608 *info = absMACRO(iinfo) + 20; 00609 return 0; 00610 } 00611 } else { 00612 /* DLARRE computes eigenvalues of the (shifted) root representation */ 00613 /* DLARRV returns the eigenvalues of the unshifted matrix. */ 00614 /* However, if the eigenvectors are not desired by the user, we need */ 00615 /* to apply the corresponding shifts from DLARRE to obtain the */ 00616 /* eigenvalues of the original matrix. */ 00617 i__1 = *m; 00618 for (j = 1; j <= i__1; ++j) { 00619 itmp = iwork[iindbl + j - 1]; 00620 w[j] += e[iwork[iinspl + itmp - 1]]; 00621 /* L20: */ 00622 } 00623 } 00624 00625 if (*tryrac) { 00626 /* Refine computed eigenvalues so that they are relatively accurate */ 00627 /* with respect to the original matrix T. */ 00628 ibegin = 1; 00629 wbegin = 1; 00630 i__1 = iwork[iindbl + *m - 1]; 00631 for (jblk = 1; jblk <= i__1; ++jblk) { 00632 iend = iwork[iinspl + jblk - 1]; 00633 in = iend - ibegin + 1; 00634 wend = wbegin - 1; 00635 /* check if any eigenvalues have to be refined in this block */ 00636 L36: 00637 if (wend < *m) { 00638 if (iwork[iindbl + wend] == jblk) { 00639 ++wend; 00640 goto L36; 00641 } 00642 } 00643 if (wend < wbegin) { 00644 ibegin = iend + 1; 00645 goto L39; 00646 } 00647 offset = iwork[iindw + wbegin - 1] - 1; 00648 ifirst = iwork[iindw + wbegin - 1]; 00649 ilast = iwork[iindw + wend - 1]; 00650 rtol2 = eps * 4.; 00651 template_lapack_larrj(&in, &work[indd + ibegin - 1], &work[inde2 + ibegin - 1], 00652 &ifirst, &ilast, &rtol2, &offset, &w[wbegin], &work[ 00653 inderr + wbegin - 1], &work[indwrk], &iwork[iindwk], & 00654 pivmin, &tnrm, &iinfo); 00655 ibegin = iend + 1; 00656 wbegin = wend + 1; 00657 L39: 00658 ; 00659 } 00660 } 00661 00662 /* If matrix was scaled, then rescale eigenvalues appropriately. */ 00663 00664 if (scale != 1.) { 00665 d__1 = 1. / scale; 00666 template_blas_scal(m, &d__1, &w[1], &c__1); 00667 } 00668 00669 /* If eigenvalues are not in increasing order, then sort them, */ 00670 /* possibly along with eigenvectors. */ 00671 00672 if (nsplit > 1) { 00673 if (! wantz) { 00674 template_lapack_lasrt("I", m, &w[1], &iinfo); 00675 if (iinfo != 0) { 00676 *info = 3; 00677 return 0; 00678 } 00679 } else { 00680 i__1 = *m - 1; 00681 for (j = 1; j <= i__1; ++j) { 00682 i__ = 0; 00683 tmp = w[j]; 00684 i__2 = *m; 00685 for (jj = j + 1; jj <= i__2; ++jj) { 00686 if (w[jj] < tmp) { 00687 i__ = jj; 00688 tmp = w[jj]; 00689 } 00690 /* L50: */ 00691 } 00692 if (i__ != 0) { 00693 w[i__] = w[j]; 00694 w[j] = tmp; 00695 if (wantz) { 00696 template_blas_swap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * 00697 z_dim1 + 1], &c__1); 00698 itmp = isuppz[(i__ << 1) - 1]; 00699 isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1]; 00700 isuppz[(j << 1) - 1] = itmp; 00701 itmp = isuppz[i__ * 2]; 00702 isuppz[i__ * 2] = isuppz[j * 2]; 00703 isuppz[j * 2] = itmp; 00704 } 00705 } 00706 /* L60: */ 00707 } 00708 } 00709 } 00710 00711 00712 work[1] = (Treal) lwmin; 00713 iwork[1] = liwmin; 00714 return 0; 00715 00716 /* End of DSTEMR */ 00717 00718 } /* dstemr_ */ 00719 00720 00721 #endif