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_STEBZ_HEADER 00036 #define TEMPLATE_LAPACK_STEBZ_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal 00041 *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, 00042 const Treal *d__, const Treal *e, integer *m, integer *nsplit, 00043 Treal *w, integer *iblock, integer *isplit, Treal *work, 00044 integer *iwork, integer *info) 00045 { 00046 /* -- LAPACK routine (version 3.0) -- 00047 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00048 Courant Institute, Argonne National Lab, and Rice University 00049 June 30, 1999 00050 00051 00052 Purpose 00053 ======= 00054 00055 DSTEBZ computes the eigenvalues of a symmetric tridiagonal 00056 matrix T. The user may ask for all eigenvalues, all eigenvalues 00057 in the half-open interval (VL, VU], or the IL-th through IU-th 00058 eigenvalues. 00059 00060 To avoid overflow, the matrix must be scaled so that its 00061 largest element is no greater than overflow**(1/2) * 00062 underflow**(1/4) in absolute value, and for greatest 00063 accuracy, it should not be much smaller than that. 00064 00065 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00066 Matrix", Report CS41, Computer Science Dept., Stanford 00067 University, July 21, 1966. 00068 00069 Arguments 00070 ========= 00071 00072 RANGE (input) CHARACTER 00073 = 'A': ("All") all eigenvalues will be found. 00074 = 'V': ("Value") all eigenvalues in the half-open interval 00075 (VL, VU] will be found. 00076 = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 00077 entire matrix) will be found. 00078 00079 ORDER (input) CHARACTER 00080 = 'B': ("By Block") the eigenvalues will be grouped by 00081 split-off block (see IBLOCK, ISPLIT) and 00082 ordered from smallest to largest within 00083 the block. 00084 = 'E': ("Entire matrix") 00085 the eigenvalues for the entire matrix 00086 will be ordered from smallest to 00087 largest. 00088 00089 N (input) INTEGER 00090 The order of the tridiagonal matrix T. N >= 0. 00091 00092 VL (input) DOUBLE PRECISION 00093 VU (input) DOUBLE PRECISION 00094 If RANGE='V', the lower and upper bounds of the interval to 00095 be searched for eigenvalues. Eigenvalues less than or equal 00096 to VL, or greater than VU, will not be returned. VL < VU. 00097 Not referenced if RANGE = 'A' or 'I'. 00098 00099 IL (input) INTEGER 00100 IU (input) INTEGER 00101 If RANGE='I', the indices (in ascending order) of the 00102 smallest and largest eigenvalues to be returned. 00103 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00104 Not referenced if RANGE = 'A' or 'V'. 00105 00106 ABSTOL (input) DOUBLE PRECISION 00107 The absolute tolerance for the eigenvalues. An eigenvalue 00108 (or cluster) is considered to be located if it has been 00109 determined to lie in an interval whose width is ABSTOL or 00110 less. If ABSTOL is less than or equal to zero, then ULP*|T| 00111 will be used, where |T| means the 1-norm of T. 00112 00113 Eigenvalues will be computed most accurately when ABSTOL is 00114 set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00115 00116 D (input) DOUBLE PRECISION array, dimension (N) 00117 The n diagonal elements of the tridiagonal matrix T. 00118 00119 E (input) DOUBLE PRECISION array, dimension (N-1) 00120 The (n-1) off-diagonal elements of the tridiagonal matrix T. 00121 00122 M (output) INTEGER 00123 The actual number of eigenvalues found. 0 <= M <= N. 00124 (See also the description of INFO=2,3.) 00125 00126 NSPLIT (output) INTEGER 00127 The number of diagonal blocks in the matrix T. 00128 1 <= NSPLIT <= N. 00129 00130 W (output) DOUBLE PRECISION array, dimension (N) 00131 On exit, the first M elements of W will contain the 00132 eigenvalues. (DSTEBZ may use the remaining N-M elements as 00133 workspace.) 00134 00135 IBLOCK (output) INTEGER array, dimension (N) 00136 At each row/column j where E(j) is zero or small, the 00137 matrix T is considered to split into a block diagonal 00138 matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which 00139 block (from 1 to the number of blocks) the eigenvalue W(i) 00140 belongs. (DSTEBZ may use the remaining N-M elements as 00141 workspace.) 00142 00143 ISPLIT (output) INTEGER array, dimension (N) 00144 The splitting points, at which T breaks up into submatrices. 00145 The first submatrix consists of rows/columns 1 to ISPLIT(1), 00146 the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 00147 etc., and the NSPLIT-th consists of rows/columns 00148 ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 00149 (Only the first NSPLIT elements will actually be used, but 00150 since the user cannot know a priori what value NSPLIT will 00151 have, N words must be reserved for ISPLIT.) 00152 00153 WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00154 00155 IWORK (workspace) INTEGER array, dimension (3*N) 00156 00157 INFO (output) INTEGER 00158 = 0: successful exit 00159 < 0: if INFO = -i, the i-th argument had an illegal value 00160 > 0: some or all of the eigenvalues failed to converge or 00161 were not computed: 00162 =1 or 3: Bisection failed to converge for some 00163 eigenvalues; these eigenvalues are flagged by a 00164 negative block number. The effect is that the 00165 eigenvalues may not be as accurate as the 00166 absolute and relative tolerances. This is 00167 generally caused by unexpectedly inaccurate 00168 arithmetic. 00169 =2 or 3: RANGE='I' only: Not all of the eigenvalues 00170 IL:IU were found. 00171 Effect: M < IU+1-IL 00172 Cause: non-monotonic arithmetic, causing the 00173 Sturm sequence to be non-monotonic. 00174 Cure: recalculate, using RANGE='A', and pick 00175 out eigenvalues IL:IU. In some cases, 00176 increasing the PARAMETER "FUDGE" may 00177 make things work. 00178 = 4: RANGE='I', and the Gershgorin interval 00179 initially used was too small. No eigenvalues 00180 were computed. 00181 Probable cause: your machine has sloppy 00182 floating-point arithmetic. 00183 Cure: Increase the PARAMETER "FUDGE", 00184 recompile, and try again. 00185 00186 Internal Parameters 00187 =================== 00188 00189 RELFAC DOUBLE PRECISION, default = 2.0e0 00190 The relative tolerance. An interval (a,b] lies within 00191 "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), 00192 where "ulp" is the machine precision (distance from 1 to 00193 the next larger floating point number.) 00194 00195 FUDGE DOUBLE PRECISION, default = 2 00196 A "fudge factor" to widen the Gershgorin intervals. Ideally, 00197 a value of 1 should work, but on machines with sloppy 00198 arithmetic, this needs to be larger. The default for 00199 publicly released versions should be large enough to handle 00200 the worst machine around. Note that this has no effect 00201 on accuracy of the solution. 00202 00203 ===================================================================== 00204 00205 00206 Parameter adjustments */ 00207 /* Table of constant values */ 00208 integer c__1 = 1; 00209 integer c_n1 = -1; 00210 integer c__3 = 3; 00211 integer c__2 = 2; 00212 integer c__0 = 0; 00213 00214 /* System generated locals */ 00215 integer i__1, i__2, i__3; 00216 Treal d__1, d__2, d__3, d__4, d__5; 00217 /* Local variables */ 00218 integer iend, ioff, iout, itmp1, j, jdisc; 00219 integer iinfo; 00220 Treal atoli; 00221 integer iwoff; 00222 Treal bnorm; 00223 integer itmax; 00224 Treal wkill, rtoli, tnorm; 00225 integer ib, jb, ie, je, nb; 00226 Treal gl; 00227 integer im, in; 00228 integer ibegin; 00229 Treal gu; 00230 integer iw; 00231 Treal wl; 00232 integer irange, idiscl; 00233 Treal safemn, wu; 00234 integer idumma[1]; 00235 integer idiscu, iorder; 00236 logical ncnvrg; 00237 Treal pivmin; 00238 logical toofew; 00239 integer nwl; 00240 Treal ulp, wlu, wul; 00241 integer nwu; 00242 Treal tmp1, tmp2; 00243 00244 00245 --iwork; 00246 --work; 00247 --isplit; 00248 --iblock; 00249 --w; 00250 --e; 00251 --d__; 00252 00253 /* Initialization added by Elias to get rid of compiler warnings. */ 00254 wlu = wul = 0; 00255 /* Function Body */ 00256 *info = 0; 00257 00258 /* Decode RANGE */ 00259 00260 if (template_blas_lsame(range, "A")) { 00261 irange = 1; 00262 } else if (template_blas_lsame(range, "V")) { 00263 irange = 2; 00264 } else if (template_blas_lsame(range, "I")) { 00265 irange = 3; 00266 } else { 00267 irange = 0; 00268 } 00269 00270 /* Decode ORDER */ 00271 00272 if (template_blas_lsame(order, "B")) { 00273 iorder = 2; 00274 } else if (template_blas_lsame(order, "E")) { 00275 iorder = 1; 00276 } else { 00277 iorder = 0; 00278 } 00279 00280 /* Check for Errors */ 00281 00282 if (irange <= 0) { 00283 *info = -1; 00284 } else if (iorder <= 0) { 00285 *info = -2; 00286 } else if (*n < 0) { 00287 *info = -3; 00288 } else if (irange == 2) { 00289 if (*vl >= *vu) { 00290 *info = -5; 00291 } 00292 } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) { 00293 *info = -6; 00294 } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) { 00295 *info = -7; 00296 } 00297 00298 if (*info != 0) { 00299 i__1 = -(*info); 00300 template_blas_erbla("STEBZ ", &i__1); 00301 return 0; 00302 } 00303 00304 /* Initialize error flags */ 00305 00306 *info = 0; 00307 ncnvrg = FALSE_; 00308 toofew = FALSE_; 00309 00310 /* Quick return if possible */ 00311 00312 *m = 0; 00313 if (*n == 0) { 00314 return 0; 00315 } 00316 00317 /* Simplifications: */ 00318 00319 if (irange == 3 && *il == 1 && *iu == *n) { 00320 irange = 1; 00321 } 00322 00323 /* Get machine constants 00324 NB is the minimum vector length for vector bisection, or 0 00325 if only scalar is to be done. */ 00326 00327 safemn = template_lapack_lamch("S", (Treal)0); 00328 ulp = template_lapack_lamch("P", (Treal)0); 00329 rtoli = ulp * 2.; 00330 nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( 00331 ftnlen)1); 00332 if (nb <= 1) { 00333 nb = 0; 00334 } 00335 00336 /* Special Case when N=1 */ 00337 00338 if (*n == 1) { 00339 *nsplit = 1; 00340 isplit[1] = 1; 00341 if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) { 00342 *m = 0; 00343 } else { 00344 w[1] = d__[1]; 00345 iblock[1] = 1; 00346 *m = 1; 00347 } 00348 return 0; 00349 } 00350 00351 /* Compute Splitting Points */ 00352 00353 *nsplit = 1; 00354 work[*n] = 0.; 00355 pivmin = 1.; 00356 00357 /* DIR$ NOVECTOR */ 00358 i__1 = *n; 00359 for (j = 2; j <= i__1; ++j) { 00360 /* Computing 2nd power */ 00361 d__1 = e[j - 1]; 00362 tmp1 = d__1 * d__1; 00363 /* Computing 2nd power */ 00364 d__2 = ulp; 00365 if ((d__1 = d__[j] * d__[j - 1], absMACRO(d__1)) * (d__2 * d__2) + safemn 00366 > tmp1) { 00367 isplit[*nsplit] = j - 1; 00368 ++(*nsplit); 00369 work[j - 1] = 0.; 00370 } else { 00371 work[j - 1] = tmp1; 00372 pivmin = maxMACRO(pivmin,tmp1); 00373 } 00374 /* L10: */ 00375 } 00376 isplit[*nsplit] = *n; 00377 pivmin *= safemn; 00378 00379 /* Compute Interval and ATOLI */ 00380 00381 if (irange == 3) { 00382 00383 /* RANGE='I': Compute the interval containing eigenvalues 00384 IL through IU. 00385 00386 Compute Gershgorin interval for entire (split) matrix 00387 and use it as the initial interval */ 00388 00389 gu = d__[1]; 00390 gl = d__[1]; 00391 tmp1 = 0.; 00392 00393 i__1 = *n - 1; 00394 for (j = 1; j <= i__1; ++j) { 00395 tmp2 = template_blas_sqrt(work[j]); 00396 /* Computing MAX */ 00397 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2; 00398 gu = maxMACRO(d__1,d__2); 00399 /* Computing MIN */ 00400 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2; 00401 gl = minMACRO(d__1,d__2); 00402 tmp1 = tmp2; 00403 /* L20: */ 00404 } 00405 00406 /* Computing MAX */ 00407 d__1 = gu, d__2 = d__[*n] + tmp1; 00408 gu = maxMACRO(d__1,d__2); 00409 /* Computing MIN */ 00410 d__1 = gl, d__2 = d__[*n] - tmp1; 00411 gl = minMACRO(d__1,d__2); 00412 /* Computing MAX */ 00413 d__1 = absMACRO(gl), d__2 = absMACRO(gu); 00414 tnorm = maxMACRO(d__1,d__2); 00415 gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.; 00416 gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.; 00417 00418 /* Compute Iteration parameters */ 00419 00420 itmax = (integer) ((template_blas_log(tnorm + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)) + 2; 00421 if (*abstol <= 0.) { 00422 atoli = ulp * tnorm; 00423 } else { 00424 atoli = *abstol; 00425 } 00426 00427 work[*n + 1] = gl; 00428 work[*n + 2] = gl; 00429 work[*n + 3] = gu; 00430 work[*n + 4] = gu; 00431 work[*n + 5] = gl; 00432 work[*n + 6] = gu; 00433 iwork[1] = -1; 00434 iwork[2] = -1; 00435 iwork[3] = *n + 1; 00436 iwork[4] = *n + 1; 00437 iwork[5] = *il - 1; 00438 iwork[6] = *iu; 00439 00440 template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin, 00441 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n 00442 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo); 00443 00444 if (iwork[6] == *iu) { 00445 wl = work[*n + 1]; 00446 wlu = work[*n + 3]; 00447 nwl = iwork[1]; 00448 wu = work[*n + 4]; 00449 wul = work[*n + 2]; 00450 nwu = iwork[4]; 00451 } else { 00452 wl = work[*n + 2]; 00453 wlu = work[*n + 4]; 00454 nwl = iwork[2]; 00455 wu = work[*n + 3]; 00456 wul = work[*n + 1]; 00457 nwu = iwork[3]; 00458 } 00459 00460 if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) { 00461 *info = 4; 00462 return 0; 00463 } 00464 } else { 00465 00466 /* RANGE='A' or 'V' -- Set ATOLI 00467 00468 Computing MAX */ 00469 d__3 = absMACRO(d__[1]) + absMACRO(e[1]), d__4 = (d__1 = d__[*n], absMACRO(d__1)) + ( 00470 d__2 = e[*n - 1], absMACRO(d__2)); 00471 tnorm = maxMACRO(d__3,d__4); 00472 00473 i__1 = *n - 1; 00474 for (j = 2; j <= i__1; ++j) { 00475 /* Computing MAX */ 00476 d__4 = tnorm, d__5 = (d__1 = d__[j], absMACRO(d__1)) + (d__2 = e[j - 1] 00477 , absMACRO(d__2)) + (d__3 = e[j], absMACRO(d__3)); 00478 tnorm = maxMACRO(d__4,d__5); 00479 /* L30: */ 00480 } 00481 00482 if (*abstol <= 0.) { 00483 atoli = ulp * tnorm; 00484 } else { 00485 atoli = *abstol; 00486 } 00487 00488 if (irange == 2) { 00489 wl = *vl; 00490 wu = *vu; 00491 } else { 00492 wl = 0.; 00493 wu = 0.; 00494 } 00495 } 00496 00497 /* Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 00498 NWL accumulates the number of eigenvalues .le. WL, 00499 NWU accumulates the number of eigenvalues .le. WU */ 00500 00501 *m = 0; 00502 iend = 0; 00503 *info = 0; 00504 nwl = 0; 00505 nwu = 0; 00506 00507 i__1 = *nsplit; 00508 for (jb = 1; jb <= i__1; ++jb) { 00509 ioff = iend; 00510 ibegin = ioff + 1; 00511 iend = isplit[jb]; 00512 in = iend - ioff; 00513 00514 if (in == 1) { 00515 00516 /* Special Case -- IN=1 */ 00517 00518 if (irange == 1 || wl >= d__[ibegin] - pivmin) { 00519 ++nwl; 00520 } 00521 if (irange == 1 || wu >= d__[ibegin] - pivmin) { 00522 ++nwu; 00523 } 00524 if (irange == 1 || ( wl < d__[ibegin] - pivmin && wu >= d__[ibegin] 00525 - pivmin ) ) { 00526 ++(*m); 00527 w[*m] = d__[ibegin]; 00528 iblock[*m] = jb; 00529 } 00530 } else { 00531 00532 /* General Case -- IN > 1 00533 00534 Compute Gershgorin Interval 00535 and use it as the initial interval */ 00536 00537 gu = d__[ibegin]; 00538 gl = d__[ibegin]; 00539 tmp1 = 0.; 00540 00541 i__2 = iend - 1; 00542 for (j = ibegin; j <= i__2; ++j) { 00543 tmp2 = (d__1 = e[j], absMACRO(d__1)); 00544 /* Computing MAX */ 00545 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2; 00546 gu = maxMACRO(d__1,d__2); 00547 /* Computing MIN */ 00548 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2; 00549 gl = minMACRO(d__1,d__2); 00550 tmp1 = tmp2; 00551 /* L40: */ 00552 } 00553 00554 /* Computing MAX */ 00555 d__1 = gu, d__2 = d__[iend] + tmp1; 00556 gu = maxMACRO(d__1,d__2); 00557 /* Computing MIN */ 00558 d__1 = gl, d__2 = d__[iend] - tmp1; 00559 gl = minMACRO(d__1,d__2); 00560 /* Computing MAX */ 00561 d__1 = absMACRO(gl), d__2 = absMACRO(gu); 00562 bnorm = maxMACRO(d__1,d__2); 00563 gl = gl - bnorm * 2. * ulp * in - pivmin * 2.; 00564 gu = gu + bnorm * 2. * ulp * in + pivmin * 2.; 00565 00566 /* Compute ATOLI for the current submatrix */ 00567 00568 if (*abstol <= 0.) { 00569 /* Computing MAX */ 00570 d__1 = absMACRO(gl), d__2 = absMACRO(gu); 00571 atoli = ulp * maxMACRO(d__1,d__2); 00572 } else { 00573 atoli = *abstol; 00574 } 00575 00576 if (irange > 1) { 00577 if (gu < wl) { 00578 nwl += in; 00579 nwu += in; 00580 goto L70; 00581 } 00582 gl = maxMACRO(gl,wl); 00583 gu = minMACRO(gu,wu); 00584 if (gl >= gu) { 00585 goto L70; 00586 } 00587 } 00588 00589 /* Set Up Initial Interval */ 00590 00591 work[*n + 1] = gl; 00592 work[*n + in + 1] = gu; 00593 template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, & 00594 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, & 00595 work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], & 00596 w[*m + 1], &iblock[*m + 1], &iinfo); 00597 00598 nwl += iwork[1]; 00599 nwu += iwork[in + 1]; 00600 iwoff = *m - iwork[1]; 00601 00602 /* Compute Eigenvalues */ 00603 00604 itmax = (integer) ((template_blas_log(gu - gl + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.) 00605 ) + 2; 00606 template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, & 00607 pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, & 00608 work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1], 00609 &w[*m + 1], &iblock[*m + 1], &iinfo); 00610 00611 /* Copy Eigenvalues Into W and IBLOCK 00612 Use -JB for block number for unconverged eigenvalues. */ 00613 00614 i__2 = iout; 00615 for (j = 1; j <= i__2; ++j) { 00616 tmp1 = (work[j + *n] + work[j + in + *n]) * .5; 00617 00618 /* Flag non-convergence. */ 00619 00620 if (j > iout - iinfo) { 00621 ncnvrg = TRUE_; 00622 ib = -jb; 00623 } else { 00624 ib = jb; 00625 } 00626 i__3 = iwork[j + in] + iwoff; 00627 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) { 00628 w[je] = tmp1; 00629 iblock[je] = ib; 00630 /* L50: */ 00631 } 00632 /* L60: */ 00633 } 00634 00635 *m += im; 00636 } 00637 L70: 00638 ; 00639 } 00640 00641 /* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00642 If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */ 00643 00644 if (irange == 3) { 00645 im = 0; 00646 idiscl = *il - 1 - nwl; 00647 idiscu = nwu - *iu; 00648 00649 if (idiscl > 0 || idiscu > 0) { 00650 i__1 = *m; 00651 for (je = 1; je <= i__1; ++je) { 00652 if (w[je] <= wlu && idiscl > 0) { 00653 --idiscl; 00654 } else if (w[je] >= wul && idiscu > 0) { 00655 --idiscu; 00656 } else { 00657 ++im; 00658 w[im] = w[je]; 00659 iblock[im] = iblock[je]; 00660 } 00661 /* L80: */ 00662 } 00663 *m = im; 00664 } 00665 if (idiscl > 0 || idiscu > 0) { 00666 00667 /* Code to deal with effects of bad arithmetic: 00668 Some low eigenvalues to be discarded are not in (WL,WLU], 00669 or high eigenvalues to be discarded are not in (WUL,WU] 00670 so just kill off the smallest IDISCL/largest IDISCU 00671 eigenvalues, by simply finding the smallest/largest 00672 eigenvalue(s). 00673 00674 (If N(w) is monotone non-decreasing, this should never 00675 happen.) */ 00676 00677 if (idiscl > 0) { 00678 wkill = wu; 00679 i__1 = idiscl; 00680 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00681 iw = 0; 00682 i__2 = *m; 00683 for (je = 1; je <= i__2; ++je) { 00684 if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) { 00685 iw = je; 00686 wkill = w[je]; 00687 } 00688 /* L90: */ 00689 } 00690 iblock[iw] = 0; 00691 /* L100: */ 00692 } 00693 } 00694 if (idiscu > 0) { 00695 00696 wkill = wl; 00697 i__1 = idiscu; 00698 for (jdisc = 1; jdisc <= i__1; ++jdisc) { 00699 iw = 0; 00700 i__2 = *m; 00701 for (je = 1; je <= i__2; ++je) { 00702 if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) { 00703 iw = je; 00704 wkill = w[je]; 00705 } 00706 /* L110: */ 00707 } 00708 iblock[iw] = 0; 00709 /* L120: */ 00710 } 00711 } 00712 im = 0; 00713 i__1 = *m; 00714 for (je = 1; je <= i__1; ++je) { 00715 if (iblock[je] != 0) { 00716 ++im; 00717 w[im] = w[je]; 00718 iblock[im] = iblock[je]; 00719 } 00720 /* L130: */ 00721 } 00722 *m = im; 00723 } 00724 if (idiscl < 0 || idiscu < 0) { 00725 toofew = TRUE_; 00726 } 00727 } 00728 00729 /* If ORDER='B', do nothing -- the eigenvalues are already sorted 00730 by block. 00731 If ORDER='E', sort the eigenvalues from smallest to largest */ 00732 00733 if (iorder == 1 && *nsplit > 1) { 00734 i__1 = *m - 1; 00735 for (je = 1; je <= i__1; ++je) { 00736 ie = 0; 00737 tmp1 = w[je]; 00738 i__2 = *m; 00739 for (j = je + 1; j <= i__2; ++j) { 00740 if (w[j] < tmp1) { 00741 ie = j; 00742 tmp1 = w[j]; 00743 } 00744 /* L140: */ 00745 } 00746 00747 if (ie != 0) { 00748 itmp1 = iblock[ie]; 00749 w[ie] = w[je]; 00750 iblock[ie] = iblock[je]; 00751 w[je] = tmp1; 00752 iblock[je] = itmp1; 00753 } 00754 /* L150: */ 00755 } 00756 } 00757 00758 *info = 0; 00759 if (ncnvrg) { 00760 ++(*info); 00761 } 00762 if (toofew) { 00763 *info += 2; 00764 } 00765 return 0; 00766 00767 /* End of DSTEBZ */ 00768 00769 } /* dstebz_ */ 00770 00771 #endif