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_LARRD_HEADER 00036 #define TEMPLATE_LAPACK_LARRD_HEADER 00037 00038 template<class Treal> 00039 int template_lapack_larrd(const char *range, const char *order, const integer *n, Treal 00040 *vl, Treal *vu, integer *il, integer *iu, Treal *gers, 00041 Treal *reltol, Treal *d__, Treal *e, Treal *e2, 00042 Treal *pivmin, integer *nsplit, integer *isplit, integer *m, 00043 Treal *w, Treal *werr, Treal *wl, Treal *wu, 00044 integer *iblock, integer *indexw, Treal *work, integer *iwork, 00045 integer *info) 00046 { 00047 /* System generated locals */ 00048 integer i__1, i__2, i__3; 00049 Treal d__1, d__2; 00050 00051 00052 /* Local variables */ 00053 integer i__, j, ib, ie, je, nb; 00054 Treal gl; 00055 integer im, in; 00056 Treal gu; 00057 integer iw, jee; 00058 Treal eps; 00059 integer nwl; 00060 Treal wlu = 0; 00061 Treal wul = 0; 00062 integer nwu; 00063 Treal tmp1, tmp2; 00064 integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc; 00065 integer iinfo; 00066 Treal atoli; 00067 integer iwoff, itmax; 00068 Treal wkill, rtoli, uflow, tnorm; 00069 integer ibegin; 00070 integer irange, idiscl, idumma[1]; 00071 integer idiscu; 00072 logical ncnvrg, toofew; 00073 00074 00075 /* -- LAPACK auxiliary routine (version 3.2.1) -- */ 00076 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00077 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ 00078 /* -- April 2009 -- */ 00079 00080 /* .. Scalar Arguments .. */ 00081 /* .. */ 00082 /* .. Array Arguments .. */ 00083 /* .. */ 00084 00085 /* Purpose */ 00086 /* ======= */ 00087 00088 /* DLARRD computes the eigenvalues of a symmetric tridiagonal */ 00089 /* matrix T to suitable accuracy. This is an auxiliary code to be */ 00090 /* called from DSTEMR. */ 00091 /* The user may ask for all eigenvalues, all eigenvalues */ 00092 /* in the half-open interval (VL, VU], or the IL-th through IU-th */ 00093 /* eigenvalues. */ 00094 00095 /* To avoid overflow, the matrix must be scaled so that its */ 00096 /* largest element is no greater than overflow**(1/2) * */ 00097 /* underflow**(1/4) in absolute value, and for greatest */ 00098 /* accuracy, it should not be much smaller than that. */ 00099 00100 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */ 00101 /* Matrix", Report CS41, Computer Science Dept., Stanford */ 00102 /* University, July 21, 1966. */ 00103 00104 /* Arguments */ 00105 /* ========= */ 00106 00107 /* RANGE (input) CHARACTER */ 00108 /* = 'A': ("All") all eigenvalues will be found. */ 00109 /* = 'V': ("Value") all eigenvalues in the half-open interval */ 00110 /* (VL, VU] will be found. */ 00111 /* = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */ 00112 /* entire matrix) will be found. */ 00113 00114 /* ORDER (input) CHARACTER */ 00115 /* = 'B': ("By Block") the eigenvalues will be grouped by */ 00116 /* split-off block (see IBLOCK, ISPLIT) and */ 00117 /* ordered from smallest to largest within */ 00118 /* the block. */ 00119 /* = 'E': ("Entire matrix") */ 00120 /* the eigenvalues for the entire matrix */ 00121 /* will be ordered from smallest to */ 00122 /* largest. */ 00123 00124 /* N (input) INTEGER */ 00125 /* The order of the tridiagonal matrix T. N >= 0. */ 00126 00127 /* VL (input) DOUBLE PRECISION */ 00128 /* VU (input) DOUBLE PRECISION */ 00129 /* If RANGE='V', the lower and upper bounds of the interval to */ 00130 /* be searched for eigenvalues. Eigenvalues less than or equal */ 00131 /* to VL, or greater than VU, will not be returned. VL < VU. */ 00132 /* Not referenced if RANGE = 'A' or 'I'. */ 00133 00134 /* IL (input) INTEGER */ 00135 /* IU (input) INTEGER */ 00136 /* If RANGE='I', the indices (in ascending order) of the */ 00137 /* smallest and largest eigenvalues to be returned. */ 00138 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */ 00139 /* Not referenced if RANGE = 'A' or 'V'. */ 00140 00141 /* GERS (input) DOUBLE PRECISION array, dimension (2*N) */ 00142 /* The N Gerschgorin intervals (the i-th Gerschgorin interval */ 00143 /* is (GERS(2*i-1), GERS(2*i)). */ 00144 00145 /* RELTOL (input) DOUBLE PRECISION */ 00146 /* The minimum relative width of an interval. When an interval */ 00147 /* is narrower than RELTOL times the larger (in */ 00148 /* magnitude) endpoint, then it is considered to be */ 00149 /* sufficiently small, i.e., converged. Note: this should */ 00150 /* always be at least radix*machine epsilon. */ 00151 00152 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00153 /* The n diagonal elements of the tridiagonal matrix T. */ 00154 00155 /* E (input) DOUBLE PRECISION array, dimension (N-1) */ 00156 /* The (n-1) off-diagonal elements of the tridiagonal matrix T. */ 00157 00158 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */ 00159 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */ 00160 00161 /* PIVMIN (input) DOUBLE PRECISION */ 00162 /* The minimum pivot allowed in the Sturm sequence for T. */ 00163 00164 /* NSPLIT (input) INTEGER */ 00165 /* The number of diagonal blocks in the matrix T. */ 00166 /* 1 <= NSPLIT <= N. */ 00167 00168 /* ISPLIT (input) INTEGER array, dimension (N) */ 00169 /* The splitting points, at which T breaks up into submatrices. */ 00170 /* The first submatrix consists of rows/columns 1 to ISPLIT(1), */ 00171 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */ 00172 /* etc., and the NSPLIT-th consists of rows/columns */ 00173 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */ 00174 /* (Only the first NSPLIT elements will actually be used, but */ 00175 /* since the user cannot know a priori what value NSPLIT will */ 00176 /* have, N words must be reserved for ISPLIT.) */ 00177 00178 /* M (output) INTEGER */ 00179 /* The actual number of eigenvalues found. 0 <= M <= N. */ 00180 /* (See also the description of INFO=2,3.) */ 00181 00182 /* W (output) DOUBLE PRECISION array, dimension (N) */ 00183 /* On exit, the first M elements of W will contain the */ 00184 /* eigenvalue approximations. DLARRD computes an interval */ 00185 /* I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */ 00186 /* approximation is given as the interval midpoint */ 00187 /* W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */ 00188 /* WERR(j) = abs( a_j - b_j)/2 */ 00189 00190 /* WERR (output) DOUBLE PRECISION array, dimension (N) */ 00191 /* The error bound on the corresponding eigenvalue approximation */ 00192 /* in W. */ 00193 00194 /* WL (output) DOUBLE PRECISION */ 00195 /* WU (output) DOUBLE PRECISION */ 00196 /* The interval (WL, WU] contains all the wanted eigenvalues. */ 00197 /* If RANGE='V', then WL=VL and WU=VU. */ 00198 /* If RANGE='A', then WL and WU are the global Gerschgorin bounds */ 00199 /* on the spectrum. */ 00200 /* If RANGE='I', then WL and WU are computed by DLAEBZ from the */ 00201 /* index range specified. */ 00202 00203 /* IBLOCK (output) INTEGER array, dimension (N) */ 00204 /* At each row/column j where E(j) is zero or small, the */ 00205 /* matrix T is considered to split into a block diagonal */ 00206 /* matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which */ 00207 /* block (from 1 to the number of blocks) the eigenvalue W(i) */ 00208 /* belongs. (DLARRD may use the remaining N-M elements as */ 00209 /* workspace.) */ 00210 00211 /* INDEXW (output) INTEGER array, dimension (N) */ 00212 /* The indices of the eigenvalues within each block (submatrix); */ 00213 /* for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */ 00214 /* i-th eigenvalue W(i) is the j-th eigenvalue in block k. */ 00215 00216 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */ 00217 00218 /* IWORK (workspace) INTEGER array, dimension (3*N) */ 00219 00220 /* INFO (output) INTEGER */ 00221 /* = 0: successful exit */ 00222 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00223 /* > 0: some or all of the eigenvalues failed to converge or */ 00224 /* were not computed: */ 00225 /* =1 or 3: Bisection failed to converge for some */ 00226 /* eigenvalues; these eigenvalues are flagged by a */ 00227 /* negative block number. The effect is that the */ 00228 /* eigenvalues may not be as accurate as the */ 00229 /* absolute and relative tolerances. This is */ 00230 /* generally caused by unexpectedly inaccurate */ 00231 /* arithmetic. */ 00232 /* =2 or 3: RANGE='I' only: Not all of the eigenvalues */ 00233 /* IL:IU were found. */ 00234 /* Effect: M < IU+1-IL */ 00235 /* Cause: non-monotonic arithmetic, causing the */ 00236 /* Sturm sequence to be non-monotonic. */ 00237 /* Cure: recalculate, using RANGE='A', and pick */ 00238 /* out eigenvalues IL:IU. In some cases, */ 00239 /* increasing the PARAMETER "FUDGE" may */ 00240 /* make things work. */ 00241 /* = 4: RANGE='I', and the Gershgorin interval */ 00242 /* initially used was too small. No eigenvalues */ 00243 /* were computed. */ 00244 /* Probable cause: your machine has sloppy */ 00245 /* floating-point arithmetic. */ 00246 /* Cure: Increase the PARAMETER "FUDGE", */ 00247 /* recompile, and try again. */ 00248 00249 /* Internal Parameters */ 00250 /* =================== */ 00251 00252 /* FUDGE DOUBLE PRECISION, default = 2 */ 00253 /* A "fudge factor" to widen the Gershgorin intervals. Ideally, */ 00254 /* a value of 1 should work, but on machines with sloppy */ 00255 /* arithmetic, this needs to be larger. The default for */ 00256 /* publicly released versions should be large enough to handle */ 00257 /* the worst machine around. Note that this has no effect */ 00258 /* on accuracy of the solution. */ 00259 00260 /* Based on contributions by */ 00261 /* W. Kahan, University of California, Berkeley, USA */ 00262 /* Beresford Parlett, University of California, Berkeley, USA */ 00263 /* Jim Demmel, University of California, Berkeley, USA */ 00264 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00265 /* Osni Marques, LBNL/NERSC, USA */ 00266 /* Christof Voemel, University of California, Berkeley, USA */ 00267 00268 /* ===================================================================== */ 00269 00270 /* .. Parameters .. */ 00271 /* .. */ 00272 /* .. Local Scalars .. */ 00273 /* .. */ 00274 /* .. Local Arrays .. */ 00275 /* .. */ 00276 /* .. External Functions .. */ 00277 /* .. */ 00278 /* .. External Subroutines .. */ 00279 /* .. */ 00280 /* .. Intrinsic Functions .. */ 00281 /* .. */ 00282 /* .. Executable Statements .. */ 00283 00284 /* Table of constant values */ 00285 00286 integer c__1 = 1; 00287 integer c_n1 = -1; 00288 integer c__3 = 3; 00289 integer c__2 = 2; 00290 integer c__0 = 0; 00291 00292 /* Parameter adjustments */ 00293 --iwork; 00294 --work; 00295 --indexw; 00296 --iblock; 00297 --werr; 00298 --w; 00299 --isplit; 00300 --e2; 00301 --e; 00302 --d__; 00303 --gers; 00304 00305 /* Function Body */ 00306 *info = 0; 00307 00308 /* Decode RANGE */ 00309 00310 if (template_blas_lsame(range, "A")) { 00311 irange = 1; 00312 } else if (template_blas_lsame(range, "V")) { 00313 irange = 2; 00314 } else if (template_blas_lsame(range, "I")) { 00315 irange = 3; 00316 } else { 00317 irange = 0; 00318 } 00319 00320 /* Check for Errors */ 00321 00322 if (irange <= 0) { 00323 *info = -1; 00324 } else if (! (template_blas_lsame(order, "B") || template_blas_lsame(order, 00325 "E"))) { 00326 *info = -2; 00327 } else if (*n < 0) { 00328 *info = -3; 00329 } else if (irange == 2) { 00330 if (*vl >= *vu) { 00331 *info = -5; 00332 } 00333 } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) { 00334 *info = -6; 00335 } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) { 00336 *info = -7; 00337 } 00338 00339 if (*info != 0) { 00340 return 0; 00341 } 00342 /* Initialize error flags */ 00343 *info = 0; 00344 ncnvrg = FALSE_; 00345 toofew = FALSE_; 00346 /* Quick return if possible */ 00347 *m = 0; 00348 if (*n == 0) { 00349 return 0; 00350 } 00351 /* Simplification: */ 00352 if (irange == 3 && *il == 1 && *iu == *n) { 00353 irange = 1; 00354 } 00355 /* Get machine constants */ 00356 eps = template_lapack_lamch("P",(Treal)0); 00357 uflow = template_lapack_lamch("U",(Treal)0); 00358 /* Special Case when N=1 */ 00359 /* Treat case of 1x1 matrix for quick return */ 00360 if (*n == 1) { 00361 if ( irange == 1 || ( irange == 2 && d__[1] > *vl && d__[1] <= *vu ) || 00362 ( irange == 3 && *il == 1 && *iu == 1 ) ) { 00363 *m = 1; 00364 w[1] = d__[1]; 00365 /* The computation error of the eigenvalue is zero */ 00366 werr[1] = 0.; 00367 iblock[1] = 1; 00368 indexw[1] = 1; 00369 } 00370 return 0; 00371 } 00372 /* NB is the minimum vector length for vector bisection, or 0 */ 00373 /* if only scalar is to be done. */ 00374 nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); 00375 if (nb <= 1) { 00376 nb = 0; 00377 } 00378 /* Find global spectral radius */ 00379 gl = d__[1]; 00380 gu = d__[1]; 00381 i__1 = *n; 00382 for (i__ = 1; i__ <= i__1; ++i__) { 00383 /* Computing MIN */ 00384 d__1 = gl, d__2 = gers[(i__ << 1) - 1]; 00385 gl = minMACRO(d__1,d__2); 00386 /* Computing MAX */ 00387 d__1 = gu, d__2 = gers[i__ * 2]; 00388 gu = maxMACRO(d__1,d__2); 00389 /* L5: */ 00390 } 00391 /* Compute global Gerschgorin bounds and spectral diameter */ 00392 /* Computing MAX */ 00393 d__1 = absMACRO(gl), d__2 = absMACRO(gu); 00394 tnorm = maxMACRO(d__1,d__2); 00395 gl = gl - tnorm * 2. * eps * *n - *pivmin * 4.; 00396 gu = gu + tnorm * 2. * eps * *n + *pivmin * 4.; 00397 /* [JAN/28/2009] remove the line below since SPDIAM variable not use */ 00398 /* SPDIAM = GU - GL */ 00399 /* Input arguments for DLAEBZ: */ 00400 /* The relative tolerance. An interval (a,b] lies within */ 00401 /* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), */ 00402 rtoli = *reltol; 00403 /* Set the absolute tolerance for interval convergence to zero to force */ 00404 /* interval convergence based on relative size of the interval. */ 00405 /* This is dangerous because intervals might not converge when RELTOL is */ 00406 /* small. But at least a very small number should be selected so that for */ 00407 /* strongly graded matrices, the code can get relatively accurate */ 00408 /* eigenvalues. */ 00409 atoli = uflow * 4. + *pivmin * 4.; 00410 if (irange == 3) { 00411 /* RANGE='I': Compute an interval containing eigenvalues */ 00412 /* IL through IU. The initial interval [GL,GU] from the global */ 00413 /* Gerschgorin bounds GL and GU is refined by DLAEBZ. */ 00414 itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 00415 2; 00416 work[*n + 1] = gl; 00417 work[*n + 2] = gl; 00418 work[*n + 3] = gu; 00419 work[*n + 4] = gu; 00420 work[*n + 5] = gl; 00421 work[*n + 6] = gu; 00422 iwork[1] = -1; 00423 iwork[2] = -1; 00424 iwork[3] = *n + 1; 00425 iwork[4] = *n + 1; 00426 iwork[5] = *il - 1; 00427 iwork[6] = *iu; 00428 00429 template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, & 00430 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5] 00431 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo); 00432 if (iinfo != 0) { 00433 *info = iinfo; 00434 return 0; 00435 } 00436 /* On exit, output intervals may not be ordered by ascending negcount */ 00437 if (iwork[6] == *iu) { 00438 *wl = work[*n + 1]; 00439 wlu = work[*n + 3]; 00440 nwl = iwork[1]; 00441 *wu = work[*n + 4]; 00442 wul = work[*n + 2]; 00443 nwu = iwork[4]; 00444 } else { 00445 *wl = work[*n + 2]; 00446 wlu = work[*n + 4]; 00447 nwl = iwork[2]; 00448 *wu = work[*n + 3]; 00449 wul = work[*n + 1]; 00450 nwu = iwork[3]; 00451 } 00452 /* On exit, the interval [WL, WLU] contains a value with negcount NWL, */ 00453 /* and [WUL, WU] contains a value with negcount NWU. */ 00454 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) { 00455 *info = 4; 00456 return 0; 00457 } 00458 } else if (irange == 2) { 00459 *wl = *vl; 00460 *wu = *vu; 00461 } else if (irange == 1) { 00462 *wl = gl; 00463 *wu = gu; 00464 } 00465 /* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */ 00466 /* NWL accumulates the number of eigenvalues .le. WL, */ 00467 /* NWU accumulates the number of eigenvalues .le. WU */ 00468 *m = 0; 00469 iend = 0; 00470 *info = 0; 00471 nwl = 0; 00472 nwu = 0; 00473 00474 i__1 = *nsplit; 00475 for (jblk = 1; jblk <= i__1; ++jblk) { 00476 ioff = iend; 00477 ibegin = ioff + 1; 00478 iend = isplit[jblk]; 00479 in = iend - ioff; 00480 00481 if (in == 1) { 00482 /* 1x1 block */ 00483 if (*wl >= d__[ibegin] - *pivmin) { 00484 ++nwl; 00485 } 00486 if (*wu >= d__[ibegin] - *pivmin) { 00487 ++nwu; 00488 } 00489 if (irange == 1 || ( *wl < d__[ibegin] - *pivmin && *wu >= d__[ 00490 ibegin] - *pivmin ) ) { 00491 ++(*m); 00492 w[*m] = d__[ibegin]; 00493 werr[*m] = 0.; 00494 /* The gap for a single block doesn't matter for the later */ 00495 /* algorithm and is assigned an arbitrary large value */ 00496 iblock[*m] = jblk; 00497 indexw[*m] = 1; 00498 } 00499 /* Disabled 2x2 case because of a failure on the following matrix */ 00500 /* RANGE = 'I', IL = IU = 4 */ 00501 /* Original Tridiagonal, d = [ */ 00502 /* -0.150102010615740E+00 */ 00503 /* -0.849897989384260E+00 */ 00504 /* -0.128208148052635E-15 */ 00505 /* 0.128257718286320E-15 */ 00506 /* ]; */ 00507 /* e = [ */ 00508 /* -0.357171383266986E+00 */ 00509 /* -0.180411241501588E-15 */ 00510 /* -0.175152352710251E-15 */ 00511 /* ]; */ 00512 00513 /* ELSE IF( IN.EQ.2 ) THEN */ 00514 /* * 2x2 block */ 00515 /* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */ 00516 /* TMP1 = HALF*(D(IBEGIN)+D(IEND)) */ 00517 /* L1 = TMP1 - DISC */ 00518 /* IF( WL.GE. L1-PIVMIN ) */ 00519 /* $ NWL = NWL + 1 */ 00520 /* IF( WU.GE. L1-PIVMIN ) */ 00521 /* $ NWU = NWU + 1 */ 00522 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */ 00523 /* $ L1-PIVMIN ) ) THEN */ 00524 /* M = M + 1 */ 00525 /* W( M ) = L1 */ 00526 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */ 00527 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */ 00528 /* IBLOCK( M ) = JBLK */ 00529 /* INDEXW( M ) = 1 */ 00530 /* ENDIF */ 00531 /* L2 = TMP1 + DISC */ 00532 /* IF( WL.GE. L2-PIVMIN ) */ 00533 /* $ NWL = NWL + 1 */ 00534 /* IF( WU.GE. L2-PIVMIN ) */ 00535 /* $ NWU = NWU + 1 */ 00536 /* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */ 00537 /* $ L2-PIVMIN ) ) THEN */ 00538 /* M = M + 1 */ 00539 /* W( M ) = L2 */ 00540 /* * The uncertainty of eigenvalues of a 2x2 matrix is very small */ 00541 /* WERR( M ) = EPS * ABS( W( M ) ) * TWO */ 00542 /* IBLOCK( M ) = JBLK */ 00543 /* INDEXW( M ) = 2 */ 00544 /* ENDIF */ 00545 } else { 00546 /* General Case - block of size IN >= 2 */ 00547 /* Compute local Gerschgorin interval and use it as the initial */ 00548 /* interval for DLAEBZ */ 00549 gu = d__[ibegin]; 00550 gl = d__[ibegin]; 00551 tmp1 = 0.; 00552 i__2 = iend; 00553 for (j = ibegin; j <= i__2; ++j) { 00554 /* Computing MIN */ 00555 d__1 = gl, d__2 = gers[(j << 1) - 1]; 00556 gl = minMACRO(d__1,d__2); 00557 /* Computing MAX */ 00558 d__1 = gu, d__2 = gers[j * 2]; 00559 gu = maxMACRO(d__1,d__2); 00560 /* L40: */ 00561 } 00562 /* [JAN/28/2009] */ 00563 /* change SPDIAM by TNORM in lines 2 and 3 thereafter */ 00564 /* line 1: remove computation of SPDIAM (not useful anymore) */ 00565 /* SPDIAM = GU - GL */ 00566 /* GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN */ 00567 /* GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN */ 00568 gl = gl - tnorm * 2. * eps * in - *pivmin * 2.; 00569 gu = gu + tnorm * 2. * eps * in + *pivmin * 2.; 00570 00571 if (irange > 1) { 00572 if (gu < *wl) { 00573 /* the local block contains none of the wanted eigenvalues */ 00574 nwl += in; 00575 nwu += in; 00576 goto L70; 00577 } 00578 /* refine search interval if possible, only range (WL,WU] matters */ 00579 gl = maxMACRO(gl,*wl); 00580 gu = minMACRO(gu,*wu); 00581 if (gl >= gu) { 00582 goto L70; 00583 } 00584 } 00585 /* Find negcount of initial interval boundaries GL and GU */ 00586 work[*n + 1] = gl; 00587 work[*n + in + 1] = gu; 00588 template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, 00589 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, & 00590 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], & 00591 w[*m + 1], &iblock[*m + 1], &iinfo); 00592 if (iinfo != 0) { 00593 *info = iinfo; 00594 return 0; 00595 } 00596 00597 nwl += iwork[1]; 00598 nwu += iwork[in + 1]; 00599 iwoff = *m - iwork[1]; 00600 /* Compute Eigenvalues */ 00601 itmax = (integer) ((template_blas_log(gu - gl + *pivmin) - template_blas_log(*pivmin)) / template_blas_log( 00602 2.)) + 2; 00603 template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, 00604 pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, & 00605 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1], 00606 &w[*m + 1], &iblock[*m + 1], &iinfo); 00607 if (iinfo != 0) { 00608 *info = iinfo; 00609 return 0; 00610 } 00611 00612 /* Copy eigenvalues into W and IBLOCK */ 00613 /* Use -JBLK for block number for unconverged eigenvalues. */ 00614 /* Loop over the number of output intervals from DLAEBZ */ 00615 i__2 = iout; 00616 for (j = 1; j <= i__2; ++j) { 00617 /* eigenvalue approximation is middle point of interval */ 00618 tmp1 = (work[j + *n] + work[j + in + *n]) * .5; 00619 /* semi length of error interval */ 00620 tmp2 = (d__1 = work[j + *n] - work[j + in + *n], absMACRO(d__1)) * 00621 .5; 00622 if (j > iout - iinfo) { 00623 /* Flag non-convergence. */ 00624 ncnvrg = TRUE_; 00625 ib = -jblk; 00626 } else { 00627 ib = jblk; 00628 } 00629 i__3 = iwork[j + in] + iwoff; 00630 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) { 00631 w[je] = tmp1; 00632 werr[je] = tmp2; 00633 indexw[je] = je - iwoff; 00634 iblock[je] = ib; 00635 /* L50: */ 00636 } 00637 /* L60: */ 00638 } 00639 00640 *m += im; 00641 } 00642 L70: 00643 ; 00644 } 00645 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */ 00646 /* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */ 00647 if (irange == 3) { 00648 idiscl = *il - 1 - nwl; 00649 idiscu = nwu - *iu; 00650 00651 if (idiscl > 0) { 00652 im = 0; 00653 i__1 = *m; 00654 for (je = 1; je <= i__1; ++je) { 00655 /* Remove some of the smallest eigenvalues from the left so that */ 00656 /* at the end IDISCL =0. Move all eigenvalues up to the left. */ 00657 if (w[je] <= wlu && idiscl > 0) { 00658 --idiscl; 00659 } else { 00660 ++im; 00661 w[im] = w[je]; 00662 werr[im] = werr[je]; 00663 indexw[im] = indexw[je]; 00664 iblock[im] = iblock[je]; 00665 } 00666 /* L80: */ 00667 } 00668 *m = im; 00669 } 00670 if (idiscu > 0) { 00671 /* Remove some of the largest eigenvalues from the right so that */ 00672 /* at the end IDISCU =0. Move all eigenvalues up to the left. */ 00673 im = *m + 1; 00674 for (je = *m; je >= 1; --je) { 00675 if (w[je] >= wul && idiscu > 0) { 00676 --idiscu; 00677 } else { 00678 --im; 00679 w[im] = w[je]; 00680 werr[im] = werr[je]; 00681 indexw[im] = indexw[je]; 00682 iblock[im] = iblock[je]; 00683 } 00684 /* L81: */ 00685 } 00686 jee = 0; 00687 i__1 = *m; 00688 for (je = im; je <= i__1; ++je) { 00689 ++jee; 00690 w[jee] = w[je]; 00691 werr[jee] = werr[je]; 00692 indexw[jee] = indexw[je]; 00693 iblock[jee] = iblock[je]; 00694 /* L82: */ 00695 } 00696 *m = *m - im + 1; 00697 } 00698 if (idiscl > 0 || idiscu > 0) { 00699 /* Code to deal with effects of bad arithmetic. (If N(w) is */ 00700 /* monotone non-decreasing, this should never happen.) */ 00701 /* Some low eigenvalues to be discarded are not in (WL,WLU], */ 00702 /* or high eigenvalues to be discarded are not in (WUL,WU] */ 00703 /* so just kill off the smallest IDISCL/largest IDISCU */ 00704 /* eigenvalues, by marking the corresponding IBLOCK = 0 */ 00705 if (idiscl > 0) { 00706 wkill = *wu; 00707 i__1 = idiscl; 00708 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00709 iw = 0; 00710 i__2 = *m; 00711 for (je = 1; je <= i__2; ++je) { 00712 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) { 00713 iw = je; 00714 wkill = w[je]; 00715 } 00716 /* L90: */ 00717 } 00718 iblock[iw] = 0; 00719 /* L100: */ 00720 } 00721 } 00722 if (idiscu > 0) { 00723 wkill = *wl; 00724 i__1 = idiscu; 00725 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00726 iw = 0; 00727 i__2 = *m; 00728 for (je = 1; je <= i__2; ++je) { 00729 if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) { 00730 iw = je; 00731 wkill = w[je]; 00732 } 00733 /* L110: */ 00734 } 00735 iblock[iw] = 0; 00736 /* L120: */ 00737 } 00738 } 00739 /* Now erase all eigenvalues with IBLOCK set to zero */ 00740 im = 0; 00741 i__1 = *m; 00742 for (je = 1; je <= i__1; ++je) { 00743 if (iblock[je] != 0) { 00744 ++im; 00745 w[im] = w[je]; 00746 werr[im] = werr[je]; 00747 indexw[im] = indexw[je]; 00748 iblock[im] = iblock[je]; 00749 } 00750 /* L130: */ 00751 } 00752 *m = im; 00753 } 00754 if (idiscl < 0 || idiscu < 0) { 00755 toofew = TRUE_; 00756 } 00757 } 00758 00759 if ( ( irange == 1 && *m != *n ) || ( irange == 3 && *m != *iu - *il + 1 ) ) { 00760 toofew = TRUE_; 00761 } 00762 /* If ORDER='B', do nothing the eigenvalues are already sorted by */ 00763 /* block. */ 00764 /* If ORDER='E', sort the eigenvalues from smallest to largest */ 00765 if (template_blas_lsame(order, "E") && *nsplit > 1) { 00766 i__1 = *m - 1; 00767 for (je = 1; je <= i__1; ++je) { 00768 ie = 0; 00769 tmp1 = w[je]; 00770 i__2 = *m; 00771 for (j = je + 1; j <= i__2; ++j) { 00772 if (w[j] < tmp1) { 00773 ie = j; 00774 tmp1 = w[j]; 00775 } 00776 /* L140: */ 00777 } 00778 if (ie != 0) { 00779 tmp2 = werr[ie]; 00780 itmp1 = iblock[ie]; 00781 itmp2 = indexw[ie]; 00782 w[ie] = w[je]; 00783 werr[ie] = werr[je]; 00784 iblock[ie] = iblock[je]; 00785 indexw[ie] = indexw[je]; 00786 w[je] = tmp1; 00787 werr[je] = tmp2; 00788 iblock[je] = itmp1; 00789 indexw[je] = itmp2; 00790 } 00791 /* L150: */ 00792 } 00793 } 00794 00795 *info = 0; 00796 if (ncnvrg) { 00797 ++(*info); 00798 } 00799 if (toofew) { 00800 *info += 2; 00801 } 00802 return 0; 00803 00804 /* End of DLARRD */ 00805 00806 } /* dlarrd_ */ 00807 00808 #endif