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