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