 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SSTEBZ( 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 REAL ABSTOL, VL, VU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 00018 REAL D( * ), E( * ), W( * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * SSTEBZ 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) REAL 00062 * VU (input) REAL 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) REAL 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*SLAMCH('S'), not zero. 00084 * 00085 * D (input) REAL array, dimension (N) 00086 * The n diagonal elements of the tridiagonal matrix T. 00087 * 00088 * E (input) REAL 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) REAL array, dimension (N) 00100 * On exit, the first M elements of W will contain the 00101 * eigenvalues. (SSTEBZ 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. (SSTEBZ 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) REAL 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 REAL, 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 REAL, 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 REAL ZERO, ONE, TWO, HALF 00176 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, 00177 $ HALF = 1.0E0 / TWO ) 00178 REAL FUDGE, RELFAC 00179 PARAMETER ( FUDGE = 2.1E0, RELFAC = 2.0E0 ) 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 REAL 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 REAL SLAMCH 00197 EXTERNAL LSAME, ILAENV, SLAMCH 00198 * .. 00199 * .. External Subroutines .. 00200 EXTERNAL SLAEBZ, 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 ) INFO = -5 00241 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 00242 $ THEN 00243 INFO = -6 00244 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 00245 $ THEN 00246 INFO = -7 00247 END IF 00248 * 00249 IF( INFO.NE.0 ) THEN 00250 CALL XERBLA( 'SSTEBZ', -INFO ) 00251 RETURN 00252 END IF 00253 * 00254 * Initialize error flags 00255 * 00256 INFO = 0 00257 NCNVRG = .FALSE. 00258 TOOFEW = .FALSE. 00259 * 00260 * Quick return if possible 00261 * 00262 M = 0 00263 IF( N.EQ.0 ) 00264 $ RETURN 00265 * 00266 * Simplifications: 00267 * 00268 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) 00269 $ IRANGE = 1 00270 * 00271 * Get machine constants 00272 * NB is the minimum vector length for vector bisection, or 0 00273 * if only scalar is to be done. 00274 * 00275 SAFEMN = SLAMCH( 'S' ) 00276 ULP = SLAMCH( 'P' ) 00277 RTOLI = ULP*RELFAC 00278 NB = ILAENV( 1, 'SSTEBZ', ' ', N, -1, -1, -1 ) 00279 IF( NB.LE.1 ) 00280 $ NB = 0 00281 * 00282 * Special Case when N=1 00283 * 00284 IF( N.EQ.1 ) THEN 00285 NSPLIT = 1 00286 ISPLIT( 1 ) = 1 00287 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN 00288 M = 0 00289 ELSE 00290 W( 1 ) = D( 1 ) 00291 IBLOCK( 1 ) = 1 00292 M = 1 00293 END IF 00294 RETURN 00295 END IF 00296 * 00297 * Compute Splitting Points 00298 * 00299 NSPLIT = 1 00300 WORK( N ) = ZERO 00301 PIVMIN = ONE 00302 * 00303 DO 10 J = 2, N 00304 TMP1 = E( J-1 )**2 00305 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 00306 ISPLIT( NSPLIT ) = J - 1 00307 NSPLIT = NSPLIT + 1 00308 WORK( J-1 ) = ZERO 00309 ELSE 00310 WORK( J-1 ) = TMP1 00311 PIVMIN = MAX( PIVMIN, TMP1 ) 00312 END IF 00313 10 CONTINUE 00314 ISPLIT( NSPLIT ) = N 00315 PIVMIN = PIVMIN*SAFEMN 00316 * 00317 * Compute Interval and ATOLI 00318 * 00319 IF( IRANGE.EQ.3 ) THEN 00320 * 00321 * RANGE='I': Compute the interval containing eigenvalues 00322 * IL through IU. 00323 * 00324 * Compute Gershgorin interval for entire (split) matrix 00325 * and use it as the initial interval 00326 * 00327 GU = D( 1 ) 00328 GL = D( 1 ) 00329 TMP1 = ZERO 00330 * 00331 DO 20 J = 1, N - 1 00332 TMP2 = SQRT( WORK( J ) ) 00333 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00334 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00335 TMP1 = TMP2 00336 20 CONTINUE 00337 * 00338 GU = MAX( GU, D( N )+TMP1 ) 00339 GL = MIN( GL, D( N )-TMP1 ) 00340 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00341 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 00342 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 00343 * 00344 * Compute Iteration parameters 00345 * 00346 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00347 $ LOG( TWO ) ) + 2 00348 IF( ABSTOL.LE.ZERO ) THEN 00349 ATOLI = ULP*TNORM 00350 ELSE 00351 ATOLI = ABSTOL 00352 END IF 00353 * 00354 WORK( N+1 ) = GL 00355 WORK( N+2 ) = GL 00356 WORK( N+3 ) = GU 00357 WORK( N+4 ) = GU 00358 WORK( N+5 ) = GL 00359 WORK( N+6 ) = GU 00360 IWORK( 1 ) = -1 00361 IWORK( 2 ) = -1 00362 IWORK( 3 ) = N + 1 00363 IWORK( 4 ) = N + 1 00364 IWORK( 5 ) = IL - 1 00365 IWORK( 6 ) = IU 00366 * 00367 CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 00368 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 00369 $ IWORK, W, IBLOCK, IINFO ) 00370 * 00371 IF( IWORK( 6 ).EQ.IU ) THEN 00372 WL = WORK( N+1 ) 00373 WLU = WORK( N+3 ) 00374 NWL = IWORK( 1 ) 00375 WU = WORK( N+4 ) 00376 WUL = WORK( N+2 ) 00377 NWU = IWORK( 4 ) 00378 ELSE 00379 WL = WORK( N+2 ) 00380 WLU = WORK( N+4 ) 00381 NWL = IWORK( 2 ) 00382 WU = WORK( N+3 ) 00383 WUL = WORK( N+1 ) 00384 NWU = IWORK( 3 ) 00385 END IF 00386 * 00387 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 00388 INFO = 4 00389 RETURN 00390 END IF 00391 ELSE 00392 * 00393 * RANGE='A' or 'V' -- Set ATOLI 00394 * 00395 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 00396 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 00397 * 00398 DO 30 J = 2, N - 1 00399 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 00400 $ ABS( E( J ) ) ) 00401 30 CONTINUE 00402 * 00403 IF( ABSTOL.LE.ZERO ) THEN 00404 ATOLI = ULP*TNORM 00405 ELSE 00406 ATOLI = ABSTOL 00407 END IF 00408 * 00409 IF( IRANGE.EQ.2 ) THEN 00410 WL = VL 00411 WU = VU 00412 ELSE 00413 WL = ZERO 00414 WU = ZERO 00415 END IF 00416 END IF 00417 * 00418 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 00419 * NWL accumulates the number of eigenvalues .le. WL, 00420 * NWU accumulates the number of eigenvalues .le. WU 00421 * 00422 M = 0 00423 IEND = 0 00424 INFO = 0 00425 NWL = 0 00426 NWU = 0 00427 * 00428 DO 70 JB = 1, NSPLIT 00429 IOFF = IEND 00430 IBEGIN = IOFF + 1 00431 IEND = ISPLIT( JB ) 00432 IN = IEND - IOFF 00433 * 00434 IF( IN.EQ.1 ) THEN 00435 * 00436 * Special Case -- IN=1 00437 * 00438 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 00439 $ NWL = NWL + 1 00440 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 00441 $ NWU = NWU + 1 00442 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 00443 $ D( IBEGIN )-PIVMIN ) ) THEN 00444 M = M + 1 00445 W( M ) = D( IBEGIN ) 00446 IBLOCK( M ) = JB 00447 END IF 00448 ELSE 00449 * 00450 * General Case -- IN > 1 00451 * 00452 * Compute Gershgorin Interval 00453 * and use it as the initial interval 00454 * 00455 GU = D( IBEGIN ) 00456 GL = D( IBEGIN ) 00457 TMP1 = ZERO 00458 * 00459 DO 40 J = IBEGIN, IEND - 1 00460 TMP2 = ABS( E( J ) ) 00461 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00462 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00463 TMP1 = TMP2 00464 40 CONTINUE 00465 * 00466 GU = MAX( GU, D( IEND )+TMP1 ) 00467 GL = MIN( GL, D( IEND )-TMP1 ) 00468 BNORM = MAX( ABS( GL ), ABS( GU ) ) 00469 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 00470 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 00471 * 00472 * Compute ATOLI for the current submatrix 00473 * 00474 IF( ABSTOL.LE.ZERO ) THEN 00475 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 00476 ELSE 00477 ATOLI = ABSTOL 00478 END IF 00479 * 00480 IF( IRANGE.GT.1 ) THEN 00481 IF( GU.LT.WL ) THEN 00482 NWL = NWL + IN 00483 NWU = NWU + IN 00484 GO TO 70 00485 END IF 00486 GL = MAX( GL, WL ) 00487 GU = MIN( GU, WU ) 00488 IF( GL.GE.GU ) 00489 $ GO TO 70 00490 END IF 00491 * 00492 * Set Up Initial Interval 00493 * 00494 WORK( N+1 ) = GL 00495 WORK( N+IN+1 ) = GU 00496 CALL SLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00497 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00498 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 00499 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00500 * 00501 NWL = NWL + IWORK( 1 ) 00502 NWU = NWU + IWORK( IN+1 ) 00503 IWOFF = M - IWORK( 1 ) 00504 * 00505 * Compute Eigenvalues 00506 * 00507 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 00508 $ LOG( TWO ) ) + 2 00509 CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00510 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00511 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 00512 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00513 * 00514 * Copy Eigenvalues Into W and IBLOCK 00515 * Use -JB for block number for unconverged eigenvalues. 00516 * 00517 DO 60 J = 1, IOUT 00518 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 00519 * 00520 * Flag non-convergence. 00521 * 00522 IF( J.GT.IOUT-IINFO ) THEN 00523 NCNVRG = .TRUE. 00524 IB = -JB 00525 ELSE 00526 IB = JB 00527 END IF 00528 DO 50 JE = IWORK( J ) + 1 + IWOFF, 00529 $ IWORK( J+IN ) + IWOFF 00530 W( JE ) = TMP1 00531 IBLOCK( JE ) = IB 00532 50 CONTINUE 00533 60 CONTINUE 00534 * 00535 M = M + IM 00536 END IF 00537 70 CONTINUE 00538 * 00539 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00540 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 00541 * 00542 IF( IRANGE.EQ.3 ) THEN 00543 IM = 0 00544 IDISCL = IL - 1 - NWL 00545 IDISCU = NWU - IU 00546 * 00547 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00548 DO 80 JE = 1, M 00549 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 00550 IDISCL = IDISCL - 1 00551 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 00552 IDISCU = IDISCU - 1 00553 ELSE 00554 IM = IM + 1 00555 W( IM ) = W( JE ) 00556 IBLOCK( IM ) = IBLOCK( JE ) 00557 END IF 00558 80 CONTINUE 00559 M = IM 00560 END IF 00561 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00562 * 00563 * Code to deal with effects of bad arithmetic: 00564 * Some low eigenvalues to be discarded are not in (WL,WLU], 00565 * or high eigenvalues to be discarded are not in (WUL,WU] 00566 * so just kill off the smallest IDISCL/largest IDISCU 00567 * eigenvalues, by simply finding the smallest/largest 00568 * eigenvalue(s). 00569 * 00570 * (If N(w) is monotone non-decreasing, this should never 00571 * happen.) 00572 * 00573 IF( IDISCL.GT.0 ) THEN 00574 WKILL = WU 00575 DO 100 JDISC = 1, IDISCL 00576 IW = 0 00577 DO 90 JE = 1, M 00578 IF( IBLOCK( JE ).NE.0 .AND. 00579 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 00580 IW = JE 00581 WKILL = W( JE ) 00582 END IF 00583 90 CONTINUE 00584 IBLOCK( IW ) = 0 00585 100 CONTINUE 00586 END IF 00587 IF( IDISCU.GT.0 ) THEN 00588 * 00589 WKILL = WL 00590 DO 120 JDISC = 1, IDISCU 00591 IW = 0 00592 DO 110 JE = 1, M 00593 IF( IBLOCK( JE ).NE.0 .AND. 00594 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 00595 IW = JE 00596 WKILL = W( JE ) 00597 END IF 00598 110 CONTINUE 00599 IBLOCK( IW ) = 0 00600 120 CONTINUE 00601 END IF 00602 IM = 0 00603 DO 130 JE = 1, M 00604 IF( IBLOCK( JE ).NE.0 ) THEN 00605 IM = IM + 1 00606 W( IM ) = W( JE ) 00607 IBLOCK( IM ) = IBLOCK( JE ) 00608 END IF 00609 130 CONTINUE 00610 M = IM 00611 END IF 00612 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 00613 TOOFEW = .TRUE. 00614 END IF 00615 END IF 00616 * 00617 * If ORDER='B', do nothing -- the eigenvalues are already sorted 00618 * by block. 00619 * If ORDER='E', sort the eigenvalues from smallest to largest 00620 * 00621 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 00622 DO 150 JE = 1, M - 1 00623 IE = 0 00624 TMP1 = W( JE ) 00625 DO 140 J = JE + 1, M 00626 IF( W( J ).LT.TMP1 ) THEN 00627 IE = J 00628 TMP1 = W( J ) 00629 END IF 00630 140 CONTINUE 00631 * 00632 IF( IE.NE.0 ) THEN 00633 ITMP1 = IBLOCK( IE ) 00634 W( IE ) = W( JE ) 00635 IBLOCK( IE ) = IBLOCK( JE ) 00636 W( JE ) = TMP1 00637 IBLOCK( JE ) = ITMP1 00638 END IF 00639 150 CONTINUE 00640 END IF 00641 * 00642 INFO = 0 00643 IF( NCNVRG ) 00644 $ INFO = INFO + 1 00645 IF( TOOFEW ) 00646 $ INFO = INFO + 2 00647 RETURN 00648 * 00649 * End of SSTEBZ 00650 * 00651 END
1.7.3 
