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