 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, 00002 $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, 00003 $ INFO ) 00004 * 00005 * -- LAPACK routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 8-18-00: Increase FUDGE factor for T3E (eca) 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER ORDER, RANGE 00013 INTEGER IL, INFO, IU, M, N, NSPLIT 00014 DOUBLE PRECISION ABSTOL, VL, VU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 00018 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DSTEBZ computes the eigenvalues of a symmetric tridiagonal 00025 * matrix T. The user may ask for all eigenvalues, all eigenvalues 00026 * in the half-open interval (VL, VU], or the IL-th through IU-th 00027 * eigenvalues. 00028 * 00029 * To avoid overflow, the matrix must be scaled so that its 00030 * largest element is no greater than overflow**(1/2) * 00031 * underflow**(1/4) in absolute value, and for greatest 00032 * accuracy, it should not be much smaller than that. 00033 * 00034 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00035 * Matrix", Report CS41, Computer Science Dept., Stanford 00036 * University, July 21, 1966. 00037 * 00038 * Arguments 00039 * ========= 00040 * 00041 * RANGE (input) CHARACTER*1 00042 * = 'A': ("All") all eigenvalues will be found. 00043 * = 'V': ("Value") all eigenvalues in the half-open interval 00044 * (VL, VU] will be found. 00045 * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 00046 * entire matrix) will be found. 00047 * 00048 * ORDER (input) CHARACTER*1 00049 * = 'B': ("By Block") the eigenvalues will be grouped by 00050 * split-off block (see IBLOCK, ISPLIT) and 00051 * ordered from smallest to largest within 00052 * the block. 00053 * = 'E': ("Entire matrix") 00054 * the eigenvalues for the entire matrix 00055 * will be ordered from smallest to 00056 * largest. 00057 * 00058 * N (input) INTEGER 00059 * The order of the tridiagonal matrix T. N >= 0. 00060 * 00061 * VL (input) DOUBLE PRECISION 00062 * VU (input) DOUBLE PRECISION 00063 * If RANGE='V', the lower and upper bounds of the interval to 00064 * be searched for eigenvalues. Eigenvalues less than or equal 00065 * to VL, or greater than VU, will not be returned. VL < VU. 00066 * Not referenced if RANGE = 'A' or 'I'. 00067 * 00068 * IL (input) INTEGER 00069 * IU (input) INTEGER 00070 * If RANGE='I', the indices (in ascending order) of the 00071 * smallest and largest eigenvalues to be returned. 00072 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00073 * Not referenced if RANGE = 'A' or 'V'. 00074 * 00075 * ABSTOL (input) DOUBLE PRECISION 00076 * The absolute tolerance for the eigenvalues. An eigenvalue 00077 * (or cluster) is considered to be located if it has been 00078 * determined to lie in an interval whose width is ABSTOL or 00079 * less. If ABSTOL is less than or equal to zero, then ULP*|T| 00080 * will be used, where |T| means the 1-norm of T. 00081 * 00082 * Eigenvalues will be computed most accurately when ABSTOL is 00083 * set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00084 * 00085 * D (input) DOUBLE PRECISION array, dimension (N) 00086 * The n diagonal elements of the tridiagonal matrix T. 00087 * 00088 * E (input) DOUBLE PRECISION array, dimension (N-1) 00089 * The (n-1) off-diagonal elements of the tridiagonal matrix T. 00090 * 00091 * M (output) INTEGER 00092 * The actual number of eigenvalues found. 0 <= M <= N. 00093 * (See also the description of INFO=2,3.) 00094 * 00095 * NSPLIT (output) INTEGER 00096 * The number of diagonal blocks in the matrix T. 00097 * 1 <= NSPLIT <= N. 00098 * 00099 * W (output) DOUBLE PRECISION array, dimension (N) 00100 * On exit, the first M elements of W will contain the 00101 * eigenvalues. (DSTEBZ may use the remaining N-M elements as 00102 * workspace.) 00103 * 00104 * IBLOCK (output) INTEGER array, dimension (N) 00105 * At each row/column j where E(j) is zero or small, the 00106 * matrix T is considered to split into a block diagonal 00107 * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which 00108 * block (from 1 to the number of blocks) the eigenvalue W(i) 00109 * belongs. (DSTEBZ may use the remaining N-M elements as 00110 * workspace.) 00111 * 00112 * ISPLIT (output) INTEGER array, dimension (N) 00113 * The splitting points, at which T breaks up into submatrices. 00114 * The first submatrix consists of rows/columns 1 to ISPLIT(1), 00115 * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 00116 * etc., and the NSPLIT-th consists of rows/columns 00117 * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 00118 * (Only the first NSPLIT elements will actually be used, but 00119 * since the user cannot know a priori what value NSPLIT will 00120 * have, N words must be reserved for ISPLIT.) 00121 * 00122 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00123 * 00124 * IWORK (workspace) INTEGER array, dimension (3*N) 00125 * 00126 * INFO (output) INTEGER 00127 * = 0: successful exit 00128 * < 0: if INFO = -i, the i-th argument had an illegal value 00129 * > 0: some or all of the eigenvalues failed to converge or 00130 * were not computed: 00131 * =1 or 3: Bisection failed to converge for some 00132 * eigenvalues; these eigenvalues are flagged by a 00133 * negative block number. The effect is that the 00134 * eigenvalues may not be as accurate as the 00135 * absolute and relative tolerances. This is 00136 * generally caused by unexpectedly inaccurate 00137 * arithmetic. 00138 * =2 or 3: RANGE='I' only: Not all of the eigenvalues 00139 * IL:IU were found. 00140 * Effect: M < IU+1-IL 00141 * Cause: non-monotonic arithmetic, causing the 00142 * Sturm sequence to be non-monotonic. 00143 * Cure: recalculate, using RANGE='A', and pick 00144 * out eigenvalues IL:IU. In some cases, 00145 * increasing the PARAMETER "FUDGE" may 00146 * make things work. 00147 * = 4: RANGE='I', and the Gershgorin interval 00148 * initially used was too small. No eigenvalues 00149 * were computed. 00150 * Probable cause: your machine has sloppy 00151 * floating-point arithmetic. 00152 * Cure: Increase the PARAMETER "FUDGE", 00153 * recompile, and try again. 00154 * 00155 * Internal Parameters 00156 * =================== 00157 * 00158 * RELFAC DOUBLE PRECISION, default = 2.0e0 00159 * The relative tolerance. An interval (a,b] lies within 00160 * "relative tolerance" if b-a < RELFAC*ulp*max(|a|,|b|), 00161 * where "ulp" is the machine precision (distance from 1 to 00162 * the next larger floating point number.) 00163 * 00164 * FUDGE DOUBLE PRECISION, default = 2 00165 * A "fudge factor" to widen the Gershgorin intervals. Ideally, 00166 * a value of 1 should work, but on machines with sloppy 00167 * arithmetic, this needs to be larger. The default for 00168 * publicly released versions should be large enough to handle 00169 * the worst machine around. Note that this has no effect 00170 * on accuracy of the solution. 00171 * 00172 * ===================================================================== 00173 * 00174 * .. Parameters .. 00175 DOUBLE PRECISION ZERO, ONE, TWO, HALF 00176 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00177 $ HALF = 1.0D0 / TWO ) 00178 DOUBLE PRECISION FUDGE, RELFAC 00179 PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 ) 00180 * .. 00181 * .. Local Scalars .. 00182 LOGICAL NCNVRG, TOOFEW 00183 INTEGER IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, 00184 $ IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX, 00185 $ ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL, 00186 $ NWU 00187 DOUBLE PRECISION ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN, 00188 $ TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL 00189 * .. 00190 * .. Local Arrays .. 00191 INTEGER IDUMMA( 1 ) 00192 * .. 00193 * .. External Functions .. 00194 LOGICAL LSAME 00195 INTEGER ILAENV 00196 DOUBLE PRECISION DLAMCH 00197 EXTERNAL LSAME, ILAENV, DLAMCH 00198 * .. 00199 * .. External Subroutines .. 00200 EXTERNAL DLAEBZ, XERBLA 00201 * .. 00202 * .. Intrinsic Functions .. 00203 INTRINSIC ABS, INT, LOG, MAX, MIN, SQRT 00204 * .. 00205 * .. Executable Statements .. 00206 * 00207 INFO = 0 00208 * 00209 * Decode RANGE 00210 * 00211 IF( LSAME( RANGE, 'A' ) ) THEN 00212 IRANGE = 1 00213 ELSE IF( LSAME( RANGE, 'V' ) ) THEN 00214 IRANGE = 2 00215 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 00216 IRANGE = 3 00217 ELSE 00218 IRANGE = 0 00219 END IF 00220 * 00221 * Decode ORDER 00222 * 00223 IF( LSAME( ORDER, 'B' ) ) THEN 00224 IORDER = 2 00225 ELSE IF( LSAME( ORDER, 'E' ) ) THEN 00226 IORDER = 1 00227 ELSE 00228 IORDER = 0 00229 END IF 00230 * 00231 * Check for Errors 00232 * 00233 IF( IRANGE.LE.0 ) THEN 00234 INFO = -1 00235 ELSE IF( IORDER.LE.0 ) THEN 00236 INFO = -2 00237 ELSE IF( N.LT.0 ) THEN 00238 INFO = -3 00239 ELSE IF( IRANGE.EQ.2 ) THEN 00240 IF( VL.GE.VU ) 00241 $ INFO = -5 00242 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 00243 $ THEN 00244 INFO = -6 00245 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 00246 $ THEN 00247 INFO = -7 00248 END IF 00249 * 00250 IF( INFO.NE.0 ) THEN 00251 CALL XERBLA( 'DSTEBZ', -INFO ) 00252 RETURN 00253 END IF 00254 * 00255 * Initialize error flags 00256 * 00257 INFO = 0 00258 NCNVRG = .FALSE. 00259 TOOFEW = .FALSE. 00260 * 00261 * Quick return if possible 00262 * 00263 M = 0 00264 IF( N.EQ.0 ) 00265 $ RETURN 00266 * 00267 * Simplifications: 00268 * 00269 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) 00270 $ IRANGE = 1 00271 * 00272 * Get machine constants 00273 * NB is the minimum vector length for vector bisection, or 0 00274 * if only scalar is to be done. 00275 * 00276 SAFEMN = DLAMCH( 'S' ) 00277 ULP = DLAMCH( 'P' ) 00278 RTOLI = ULP*RELFAC 00279 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) 00280 IF( NB.LE.1 ) 00281 $ NB = 0 00282 * 00283 * Special Case when N=1 00284 * 00285 IF( N.EQ.1 ) THEN 00286 NSPLIT = 1 00287 ISPLIT( 1 ) = 1 00288 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN 00289 M = 0 00290 ELSE 00291 W( 1 ) = D( 1 ) 00292 IBLOCK( 1 ) = 1 00293 M = 1 00294 END IF 00295 RETURN 00296 END IF 00297 * 00298 * Compute Splitting Points 00299 * 00300 NSPLIT = 1 00301 WORK( N ) = ZERO 00302 PIVMIN = ONE 00303 * 00304 DO 10 J = 2, N 00305 TMP1 = E( J-1 )**2 00306 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 00307 ISPLIT( NSPLIT ) = J - 1 00308 NSPLIT = NSPLIT + 1 00309 WORK( J-1 ) = ZERO 00310 ELSE 00311 WORK( J-1 ) = TMP1 00312 PIVMIN = MAX( PIVMIN, TMP1 ) 00313 END IF 00314 10 CONTINUE 00315 ISPLIT( NSPLIT ) = N 00316 PIVMIN = PIVMIN*SAFEMN 00317 * 00318 * Compute Interval and ATOLI 00319 * 00320 IF( IRANGE.EQ.3 ) THEN 00321 * 00322 * RANGE='I': Compute the interval containing eigenvalues 00323 * IL through IU. 00324 * 00325 * Compute Gershgorin interval for entire (split) matrix 00326 * and use it as the initial interval 00327 * 00328 GU = D( 1 ) 00329 GL = D( 1 ) 00330 TMP1 = ZERO 00331 * 00332 DO 20 J = 1, N - 1 00333 TMP2 = SQRT( WORK( J ) ) 00334 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00335 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00336 TMP1 = TMP2 00337 20 CONTINUE 00338 * 00339 GU = MAX( GU, D( N )+TMP1 ) 00340 GL = MIN( GL, D( N )-TMP1 ) 00341 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00342 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 00343 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 00344 * 00345 * Compute Iteration parameters 00346 * 00347 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00348 $ LOG( TWO ) ) + 2 00349 IF( ABSTOL.LE.ZERO ) THEN 00350 ATOLI = ULP*TNORM 00351 ELSE 00352 ATOLI = ABSTOL 00353 END IF 00354 * 00355 WORK( N+1 ) = GL 00356 WORK( N+2 ) = GL 00357 WORK( N+3 ) = GU 00358 WORK( N+4 ) = GU 00359 WORK( N+5 ) = GL 00360 WORK( N+6 ) = GU 00361 IWORK( 1 ) = -1 00362 IWORK( 2 ) = -1 00363 IWORK( 3 ) = N + 1 00364 IWORK( 4 ) = N + 1 00365 IWORK( 5 ) = IL - 1 00366 IWORK( 6 ) = IU 00367 * 00368 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 00369 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 00370 $ IWORK, W, IBLOCK, IINFO ) 00371 * 00372 IF( IWORK( 6 ).EQ.IU ) THEN 00373 WL = WORK( N+1 ) 00374 WLU = WORK( N+3 ) 00375 NWL = IWORK( 1 ) 00376 WU = WORK( N+4 ) 00377 WUL = WORK( N+2 ) 00378 NWU = IWORK( 4 ) 00379 ELSE 00380 WL = WORK( N+2 ) 00381 WLU = WORK( N+4 ) 00382 NWL = IWORK( 2 ) 00383 WU = WORK( N+3 ) 00384 WUL = WORK( N+1 ) 00385 NWU = IWORK( 3 ) 00386 END IF 00387 * 00388 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 00389 INFO = 4 00390 RETURN 00391 END IF 00392 ELSE 00393 * 00394 * RANGE='A' or 'V' -- Set ATOLI 00395 * 00396 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 00397 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 00398 * 00399 DO 30 J = 2, N - 1 00400 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 00401 $ ABS( E( J ) ) ) 00402 30 CONTINUE 00403 * 00404 IF( ABSTOL.LE.ZERO ) THEN 00405 ATOLI = ULP*TNORM 00406 ELSE 00407 ATOLI = ABSTOL 00408 END IF 00409 * 00410 IF( IRANGE.EQ.2 ) THEN 00411 WL = VL 00412 WU = VU 00413 ELSE 00414 WL = ZERO 00415 WU = ZERO 00416 END IF 00417 END IF 00418 * 00419 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 00420 * NWL accumulates the number of eigenvalues .le. WL, 00421 * NWU accumulates the number of eigenvalues .le. WU 00422 * 00423 M = 0 00424 IEND = 0 00425 INFO = 0 00426 NWL = 0 00427 NWU = 0 00428 * 00429 DO 70 JB = 1, NSPLIT 00430 IOFF = IEND 00431 IBEGIN = IOFF + 1 00432 IEND = ISPLIT( JB ) 00433 IN = IEND - IOFF 00434 * 00435 IF( IN.EQ.1 ) THEN 00436 * 00437 * Special Case -- IN=1 00438 * 00439 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 00440 $ NWL = NWL + 1 00441 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 00442 $ NWU = NWU + 1 00443 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 00444 $ D( IBEGIN )-PIVMIN ) ) THEN 00445 M = M + 1 00446 W( M ) = D( IBEGIN ) 00447 IBLOCK( M ) = JB 00448 END IF 00449 ELSE 00450 * 00451 * General Case -- IN > 1 00452 * 00453 * Compute Gershgorin Interval 00454 * and use it as the initial interval 00455 * 00456 GU = D( IBEGIN ) 00457 GL = D( IBEGIN ) 00458 TMP1 = ZERO 00459 * 00460 DO 40 J = IBEGIN, IEND - 1 00461 TMP2 = ABS( E( J ) ) 00462 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00463 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00464 TMP1 = TMP2 00465 40 CONTINUE 00466 * 00467 GU = MAX( GU, D( IEND )+TMP1 ) 00468 GL = MIN( GL, D( IEND )-TMP1 ) 00469 BNORM = MAX( ABS( GL ), ABS( GU ) ) 00470 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 00471 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 00472 * 00473 * Compute ATOLI for the current submatrix 00474 * 00475 IF( ABSTOL.LE.ZERO ) THEN 00476 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 00477 ELSE 00478 ATOLI = ABSTOL 00479 END IF 00480 * 00481 IF( IRANGE.GT.1 ) THEN 00482 IF( GU.LT.WL ) THEN 00483 NWL = NWL + IN 00484 NWU = NWU + IN 00485 GO TO 70 00486 END IF 00487 GL = MAX( GL, WL ) 00488 GU = MIN( GU, WU ) 00489 IF( GL.GE.GU ) 00490 $ GO TO 70 00491 END IF 00492 * 00493 * Set Up Initial Interval 00494 * 00495 WORK( N+1 ) = GL 00496 WORK( N+IN+1 ) = GU 00497 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00498 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00499 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 00500 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00501 * 00502 NWL = NWL + IWORK( 1 ) 00503 NWU = NWU + IWORK( IN+1 ) 00504 IWOFF = M - IWORK( 1 ) 00505 * 00506 * Compute Eigenvalues 00507 * 00508 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 00509 $ LOG( TWO ) ) + 2 00510 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00511 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00512 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 00513 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00514 * 00515 * Copy Eigenvalues Into W and IBLOCK 00516 * Use -JB for block number for unconverged eigenvalues. 00517 * 00518 DO 60 J = 1, IOUT 00519 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 00520 * 00521 * Flag non-convergence. 00522 * 00523 IF( J.GT.IOUT-IINFO ) THEN 00524 NCNVRG = .TRUE. 00525 IB = -JB 00526 ELSE 00527 IB = JB 00528 END IF 00529 DO 50 JE = IWORK( J ) + 1 + IWOFF, 00530 $ IWORK( J+IN ) + IWOFF 00531 W( JE ) = TMP1 00532 IBLOCK( JE ) = IB 00533 50 CONTINUE 00534 60 CONTINUE 00535 * 00536 M = M + IM 00537 END IF 00538 70 CONTINUE 00539 * 00540 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00541 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 00542 * 00543 IF( IRANGE.EQ.3 ) THEN 00544 IM = 0 00545 IDISCL = IL - 1 - NWL 00546 IDISCU = NWU - IU 00547 * 00548 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00549 DO 80 JE = 1, M 00550 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 00551 IDISCL = IDISCL - 1 00552 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 00553 IDISCU = IDISCU - 1 00554 ELSE 00555 IM = IM + 1 00556 W( IM ) = W( JE ) 00557 IBLOCK( IM ) = IBLOCK( JE ) 00558 END IF 00559 80 CONTINUE 00560 M = IM 00561 END IF 00562 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00563 * 00564 * Code to deal with effects of bad arithmetic: 00565 * Some low eigenvalues to be discarded are not in (WL,WLU], 00566 * or high eigenvalues to be discarded are not in (WUL,WU] 00567 * so just kill off the smallest IDISCL/largest IDISCU 00568 * eigenvalues, by simply finding the smallest/largest 00569 * eigenvalue(s). 00570 * 00571 * (If N(w) is monotone non-decreasing, this should never 00572 * happen.) 00573 * 00574 IF( IDISCL.GT.0 ) THEN 00575 WKILL = WU 00576 DO 100 JDISC = 1, IDISCL 00577 IW = 0 00578 DO 90 JE = 1, M 00579 IF( IBLOCK( JE ).NE.0 .AND. 00580 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 00581 IW = JE 00582 WKILL = W( JE ) 00583 END IF 00584 90 CONTINUE 00585 IBLOCK( IW ) = 0 00586 100 CONTINUE 00587 END IF 00588 IF( IDISCU.GT.0 ) THEN 00589 * 00590 WKILL = WL 00591 DO 120 JDISC = 1, IDISCU 00592 IW = 0 00593 DO 110 JE = 1, M 00594 IF( IBLOCK( JE ).NE.0 .AND. 00595 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 00596 IW = JE 00597 WKILL = W( JE ) 00598 END IF 00599 110 CONTINUE 00600 IBLOCK( IW ) = 0 00601 120 CONTINUE 00602 END IF 00603 IM = 0 00604 DO 130 JE = 1, M 00605 IF( IBLOCK( JE ).NE.0 ) THEN 00606 IM = IM + 1 00607 W( IM ) = W( JE ) 00608 IBLOCK( IM ) = IBLOCK( JE ) 00609 END IF 00610 130 CONTINUE 00611 M = IM 00612 END IF 00613 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 00614 TOOFEW = .TRUE. 00615 END IF 00616 END IF 00617 * 00618 * If ORDER='B', do nothing -- the eigenvalues are already sorted 00619 * by block. 00620 * If ORDER='E', sort the eigenvalues from smallest to largest 00621 * 00622 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 00623 DO 150 JE = 1, M - 1 00624 IE = 0 00625 TMP1 = W( JE ) 00626 DO 140 J = JE + 1, M 00627 IF( W( J ).LT.TMP1 ) THEN 00628 IE = J 00629 TMP1 = W( J ) 00630 END IF 00631 140 CONTINUE 00632 * 00633 IF( IE.NE.0 ) THEN 00634 ITMP1 = IBLOCK( IE ) 00635 W( IE ) = W( JE ) 00636 IBLOCK( IE ) = IBLOCK( JE ) 00637 W( JE ) = TMP1 00638 IBLOCK( JE ) = ITMP1 00639 END IF 00640 150 CONTINUE 00641 END IF 00642 * 00643 INFO = 0 00644 IF( NCNVRG ) 00645 $ INFO = INFO + 1 00646 IF( TOOFEW ) 00647 $ INFO = INFO + 2 00648 RETURN 00649 * 00650 * End of DSTEBZ 00651 * 00652 END
1.7.3 
