|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 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 *DIR$ NOVECTOR 00305 DO 10 J = 2, N 00306 TMP1 = E( J-1 )**2 00307 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 00308 ISPLIT( NSPLIT ) = J - 1 00309 NSPLIT = NSPLIT + 1 00310 WORK( J-1 ) = ZERO 00311 ELSE 00312 WORK( J-1 ) = TMP1 00313 PIVMIN = MAX( PIVMIN, TMP1 ) 00314 END IF 00315 10 CONTINUE 00316 ISPLIT( NSPLIT ) = N 00317 PIVMIN = PIVMIN*SAFEMN 00318 * 00319 * Compute Interval and ATOLI 00320 * 00321 IF( IRANGE.EQ.3 ) THEN 00322 * 00323 * RANGE='I': Compute the interval containing eigenvalues 00324 * IL through IU. 00325 * 00326 * Compute Gershgorin interval for entire (split) matrix 00327 * and use it as the initial interval 00328 * 00329 GU = D( 1 ) 00330 GL = D( 1 ) 00331 TMP1 = ZERO 00332 * 00333 DO 20 J = 1, N - 1 00334 TMP2 = SQRT( WORK( J ) ) 00335 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00336 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00337 TMP1 = TMP2 00338 20 CONTINUE 00339 * 00340 GU = MAX( GU, D( N )+TMP1 ) 00341 GL = MIN( GL, D( N )-TMP1 ) 00342 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00343 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 00344 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 00345 * 00346 * Compute Iteration parameters 00347 * 00348 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00349 $ LOG( TWO ) ) + 2 00350 IF( ABSTOL.LE.ZERO ) THEN 00351 ATOLI = ULP*TNORM 00352 ELSE 00353 ATOLI = ABSTOL 00354 END IF 00355 * 00356 WORK( N+1 ) = GL 00357 WORK( N+2 ) = GL 00358 WORK( N+3 ) = GU 00359 WORK( N+4 ) = GU 00360 WORK( N+5 ) = GL 00361 WORK( N+6 ) = GU 00362 IWORK( 1 ) = -1 00363 IWORK( 2 ) = -1 00364 IWORK( 3 ) = N + 1 00365 IWORK( 4 ) = N + 1 00366 IWORK( 5 ) = IL - 1 00367 IWORK( 6 ) = IU 00368 * 00369 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 00370 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 00371 $ IWORK, W, IBLOCK, IINFO ) 00372 * 00373 IF( IWORK( 6 ).EQ.IU ) THEN 00374 WL = WORK( N+1 ) 00375 WLU = WORK( N+3 ) 00376 NWL = IWORK( 1 ) 00377 WU = WORK( N+4 ) 00378 WUL = WORK( N+2 ) 00379 NWU = IWORK( 4 ) 00380 ELSE 00381 WL = WORK( N+2 ) 00382 WLU = WORK( N+4 ) 00383 NWL = IWORK( 2 ) 00384 WU = WORK( N+3 ) 00385 WUL = WORK( N+1 ) 00386 NWU = IWORK( 3 ) 00387 END IF 00388 * 00389 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 00390 INFO = 4 00391 RETURN 00392 END IF 00393 ELSE 00394 * 00395 * RANGE='A' or 'V' -- Set ATOLI 00396 * 00397 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 00398 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 00399 * 00400 DO 30 J = 2, N - 1 00401 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 00402 $ ABS( E( J ) ) ) 00403 30 CONTINUE 00404 * 00405 IF( ABSTOL.LE.ZERO ) THEN 00406 ATOLI = ULP*TNORM 00407 ELSE 00408 ATOLI = ABSTOL 00409 END IF 00410 * 00411 IF( IRANGE.EQ.2 ) THEN 00412 WL = VL 00413 WU = VU 00414 ELSE 00415 WL = ZERO 00416 WU = ZERO 00417 END IF 00418 END IF 00419 * 00420 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 00421 * NWL accumulates the number of eigenvalues .le. WL, 00422 * NWU accumulates the number of eigenvalues .le. WU 00423 * 00424 M = 0 00425 IEND = 0 00426 INFO = 0 00427 NWL = 0 00428 NWU = 0 00429 * 00430 DO 70 JB = 1, NSPLIT 00431 IOFF = IEND 00432 IBEGIN = IOFF + 1 00433 IEND = ISPLIT( JB ) 00434 IN = IEND - IOFF 00435 * 00436 IF( IN.EQ.1 ) THEN 00437 * 00438 * Special Case -- IN=1 00439 * 00440 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 00441 $ NWL = NWL + 1 00442 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 00443 $ NWU = NWU + 1 00444 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 00445 $ D( IBEGIN )-PIVMIN ) ) THEN 00446 M = M + 1 00447 W( M ) = D( IBEGIN ) 00448 IBLOCK( M ) = JB 00449 END IF 00450 ELSE 00451 * 00452 * General Case -- IN > 1 00453 * 00454 * Compute Gershgorin Interval 00455 * and use it as the initial interval 00456 * 00457 GU = D( IBEGIN ) 00458 GL = D( IBEGIN ) 00459 TMP1 = ZERO 00460 * 00461 DO 40 J = IBEGIN, IEND - 1 00462 TMP2 = ABS( E( J ) ) 00463 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00464 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00465 TMP1 = TMP2 00466 40 CONTINUE 00467 * 00468 GU = MAX( GU, D( IEND )+TMP1 ) 00469 GL = MIN( GL, D( IEND )-TMP1 ) 00470 BNORM = MAX( ABS( GL ), ABS( GU ) ) 00471 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 00472 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 00473 * 00474 * Compute ATOLI for the current submatrix 00475 * 00476 IF( ABSTOL.LE.ZERO ) THEN 00477 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 00478 ELSE 00479 ATOLI = ABSTOL 00480 END IF 00481 * 00482 IF( IRANGE.GT.1 ) THEN 00483 IF( GU.LT.WL ) THEN 00484 NWL = NWL + IN 00485 NWU = NWU + IN 00486 GO TO 70 00487 END IF 00488 GL = MAX( GL, WL ) 00489 GU = MIN( GU, WU ) 00490 IF( GL.GE.GU ) 00491 $ GO TO 70 00492 END IF 00493 * 00494 * Set Up Initial Interval 00495 * 00496 WORK( N+1 ) = GL 00497 WORK( N+IN+1 ) = GU 00498 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00499 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00500 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 00501 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00502 * 00503 NWL = NWL + IWORK( 1 ) 00504 NWU = NWU + IWORK( IN+1 ) 00505 IWOFF = M - IWORK( 1 ) 00506 * 00507 * Compute Eigenvalues 00508 * 00509 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 00510 $ LOG( TWO ) ) + 2 00511 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00512 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00513 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 00514 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00515 * 00516 * Copy Eigenvalues Into W and IBLOCK 00517 * Use -JB for block number for unconverged eigenvalues. 00518 * 00519 DO 60 J = 1, IOUT 00520 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 00521 * 00522 * Flag non-convergence. 00523 * 00524 IF( J.GT.IOUT-IINFO ) THEN 00525 NCNVRG = .TRUE. 00526 IB = -JB 00527 ELSE 00528 IB = JB 00529 END IF 00530 DO 50 JE = IWORK( J ) + 1 + IWOFF, 00531 $ IWORK( J+IN ) + IWOFF 00532 W( JE ) = TMP1 00533 IBLOCK( JE ) = IB 00534 50 CONTINUE 00535 60 CONTINUE 00536 * 00537 M = M + IM 00538 END IF 00539 70 CONTINUE 00540 * 00541 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00542 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 00543 * 00544 IF( IRANGE.EQ.3 ) THEN 00545 IM = 0 00546 IDISCL = IL - 1 - NWL 00547 IDISCU = NWU - IU 00548 * 00549 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00550 DO 80 JE = 1, M 00551 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 00552 IDISCL = IDISCL - 1 00553 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 00554 IDISCU = IDISCU - 1 00555 ELSE 00556 IM = IM + 1 00557 W( IM ) = W( JE ) 00558 IBLOCK( IM ) = IBLOCK( JE ) 00559 END IF 00560 80 CONTINUE 00561 M = IM 00562 END IF 00563 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00564 * 00565 * Code to deal with effects of bad arithmetic: 00566 * Some low eigenvalues to be discarded are not in (WL,WLU], 00567 * or high eigenvalues to be discarded are not in (WUL,WU] 00568 * so just kill off the smallest IDISCL/largest IDISCU 00569 * eigenvalues, by simply finding the smallest/largest 00570 * eigenvalue(s). 00571 * 00572 * (If N(w) is monotone non-decreasing, this should never 00573 * happen.) 00574 * 00575 IF( IDISCL.GT.0 ) THEN 00576 WKILL = WU 00577 DO 100 JDISC = 1, IDISCL 00578 IW = 0 00579 DO 90 JE = 1, M 00580 IF( IBLOCK( JE ).NE.0 .AND. 00581 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 00582 IW = JE 00583 WKILL = W( JE ) 00584 END IF 00585 90 CONTINUE 00586 IBLOCK( IW ) = 0 00587 100 CONTINUE 00588 END IF 00589 IF( IDISCU.GT.0 ) THEN 00590 * 00591 WKILL = WL 00592 DO 120 JDISC = 1, IDISCU 00593 IW = 0 00594 DO 110 JE = 1, M 00595 IF( IBLOCK( JE ).NE.0 .AND. 00596 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 00597 IW = JE 00598 WKILL = W( JE ) 00599 END IF 00600 110 CONTINUE 00601 IBLOCK( IW ) = 0 00602 120 CONTINUE 00603 END IF 00604 IM = 0 00605 DO 130 JE = 1, M 00606 IF( IBLOCK( JE ).NE.0 ) THEN 00607 IM = IM + 1 00608 W( IM ) = W( JE ) 00609 IBLOCK( IM ) = IBLOCK( JE ) 00610 END IF 00611 130 CONTINUE 00612 M = IM 00613 END IF 00614 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 00615 TOOFEW = .TRUE. 00616 END IF 00617 END IF 00618 * 00619 * If ORDER='B', do nothing -- the eigenvalues are already sorted 00620 * by block. 00621 * If ORDER='E', sort the eigenvalues from smallest to largest 00622 * 00623 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 00624 DO 150 JE = 1, M - 1 00625 IE = 0 00626 TMP1 = W( JE ) 00627 DO 140 J = JE + 1, M 00628 IF( W( J ).LT.TMP1 ) THEN 00629 IE = J 00630 TMP1 = W( J ) 00631 END IF 00632 140 CONTINUE 00633 * 00634 IF( IE.NE.0 ) THEN 00635 ITMP1 = IBLOCK( IE ) 00636 W( IE ) = W( JE ) 00637 IBLOCK( IE ) = IBLOCK( JE ) 00638 W( JE ) = TMP1 00639 IBLOCK( JE ) = ITMP1 00640 END IF 00641 150 CONTINUE 00642 END IF 00643 * 00644 INFO = 0 00645 IF( NCNVRG ) 00646 $ INFO = INFO + 1 00647 IF( TOOFEW ) 00648 $ INFO = INFO + 2 00649 RETURN 00650 * 00651 * End of DSTEBZ 00652 * 00653 END
1.7.3 
