![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTEBZ 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSTEBZ + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstebz.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstebz.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstebz.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSTEBZ( 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 * DOUBLE PRECISION ABSTOL, VL, VU 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 00032 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> DSTEBZ 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 DOUBLE PRECISION 00090 *> \endverbatim 00091 *> 00092 *> \param[in] VU 00093 *> \verbatim 00094 *> VU is DOUBLE PRECISION 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 DOUBLE PRECISION 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*DLAMCH('S'), not zero. 00128 *> \endverbatim 00129 *> 00130 *> \param[in] D 00131 *> \verbatim 00132 *> D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 00159 *> On exit, the first M elements of W will contain the 00160 *> eigenvalues. (DSTEBZ 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. (DSTEBZ 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 DOUBLE PRECISION 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 DOUBLE PRECISION, 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 DOUBLE PRECISION, 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 DSTEBZ( 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 DOUBLE PRECISION ABSTOL, VL, VU 00275 * .. 00276 * .. Array Arguments .. 00277 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ) 00278 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 00279 * .. 00280 * 00281 * ===================================================================== 00282 * 00283 * .. Parameters .. 00284 DOUBLE PRECISION ZERO, ONE, TWO, HALF 00285 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00286 $ HALF = 1.0D0 / TWO ) 00287 DOUBLE PRECISION FUDGE, RELFAC 00288 PARAMETER ( FUDGE = 2.1D0, RELFAC = 2.0D0 ) 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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00306 EXTERNAL LSAME, ILAENV, DLAMCH 00307 * .. 00308 * .. External Subroutines .. 00309 EXTERNAL DLAEBZ, 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 ) 00350 $ INFO = -5 00351 ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) 00352 $ THEN 00353 INFO = -6 00354 ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) 00355 $ THEN 00356 INFO = -7 00357 END IF 00358 * 00359 IF( INFO.NE.0 ) THEN 00360 CALL XERBLA( 'DSTEBZ', -INFO ) 00361 RETURN 00362 END IF 00363 * 00364 * Initialize error flags 00365 * 00366 INFO = 0 00367 NCNVRG = .FALSE. 00368 TOOFEW = .FALSE. 00369 * 00370 * Quick return if possible 00371 * 00372 M = 0 00373 IF( N.EQ.0 ) 00374 $ RETURN 00375 * 00376 * Simplifications: 00377 * 00378 IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N ) 00379 $ IRANGE = 1 00380 * 00381 * Get machine constants 00382 * NB is the minimum vector length for vector bisection, or 0 00383 * if only scalar is to be done. 00384 * 00385 SAFEMN = DLAMCH( 'S' ) 00386 ULP = DLAMCH( 'P' ) 00387 RTOLI = ULP*RELFAC 00388 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) 00389 IF( NB.LE.1 ) 00390 $ NB = 0 00391 * 00392 * Special Case when N=1 00393 * 00394 IF( N.EQ.1 ) THEN 00395 NSPLIT = 1 00396 ISPLIT( 1 ) = 1 00397 IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN 00398 M = 0 00399 ELSE 00400 W( 1 ) = D( 1 ) 00401 IBLOCK( 1 ) = 1 00402 M = 1 00403 END IF 00404 RETURN 00405 END IF 00406 * 00407 * Compute Splitting Points 00408 * 00409 NSPLIT = 1 00410 WORK( N ) = ZERO 00411 PIVMIN = ONE 00412 * 00413 DO 10 J = 2, N 00414 TMP1 = E( J-1 )**2 00415 IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN 00416 ISPLIT( NSPLIT ) = J - 1 00417 NSPLIT = NSPLIT + 1 00418 WORK( J-1 ) = ZERO 00419 ELSE 00420 WORK( J-1 ) = TMP1 00421 PIVMIN = MAX( PIVMIN, TMP1 ) 00422 END IF 00423 10 CONTINUE 00424 ISPLIT( NSPLIT ) = N 00425 PIVMIN = PIVMIN*SAFEMN 00426 * 00427 * Compute Interval and ATOLI 00428 * 00429 IF( IRANGE.EQ.3 ) THEN 00430 * 00431 * RANGE='I': Compute the interval containing eigenvalues 00432 * IL through IU. 00433 * 00434 * Compute Gershgorin interval for entire (split) matrix 00435 * and use it as the initial interval 00436 * 00437 GU = D( 1 ) 00438 GL = D( 1 ) 00439 TMP1 = ZERO 00440 * 00441 DO 20 J = 1, N - 1 00442 TMP2 = SQRT( WORK( J ) ) 00443 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00444 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00445 TMP1 = TMP2 00446 20 CONTINUE 00447 * 00448 GU = MAX( GU, D( N )+TMP1 ) 00449 GL = MIN( GL, D( N )-TMP1 ) 00450 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00451 GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 00452 GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 00453 * 00454 * Compute Iteration parameters 00455 * 00456 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00457 $ LOG( TWO ) ) + 2 00458 IF( ABSTOL.LE.ZERO ) THEN 00459 ATOLI = ULP*TNORM 00460 ELSE 00461 ATOLI = ABSTOL 00462 END IF 00463 * 00464 WORK( N+1 ) = GL 00465 WORK( N+2 ) = GL 00466 WORK( N+3 ) = GU 00467 WORK( N+4 ) = GU 00468 WORK( N+5 ) = GL 00469 WORK( N+6 ) = GU 00470 IWORK( 1 ) = -1 00471 IWORK( 2 ) = -1 00472 IWORK( 3 ) = N + 1 00473 IWORK( 4 ) = N + 1 00474 IWORK( 5 ) = IL - 1 00475 IWORK( 6 ) = IU 00476 * 00477 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, 00478 $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 00479 $ IWORK, W, IBLOCK, IINFO ) 00480 * 00481 IF( IWORK( 6 ).EQ.IU ) THEN 00482 WL = WORK( N+1 ) 00483 WLU = WORK( N+3 ) 00484 NWL = IWORK( 1 ) 00485 WU = WORK( N+4 ) 00486 WUL = WORK( N+2 ) 00487 NWU = IWORK( 4 ) 00488 ELSE 00489 WL = WORK( N+2 ) 00490 WLU = WORK( N+4 ) 00491 NWL = IWORK( 2 ) 00492 WU = WORK( N+3 ) 00493 WUL = WORK( N+1 ) 00494 NWU = IWORK( 3 ) 00495 END IF 00496 * 00497 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 00498 INFO = 4 00499 RETURN 00500 END IF 00501 ELSE 00502 * 00503 * RANGE='A' or 'V' -- Set ATOLI 00504 * 00505 TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ), 00506 $ ABS( D( N ) )+ABS( E( N-1 ) ) ) 00507 * 00508 DO 30 J = 2, N - 1 00509 TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+ 00510 $ ABS( E( J ) ) ) 00511 30 CONTINUE 00512 * 00513 IF( ABSTOL.LE.ZERO ) THEN 00514 ATOLI = ULP*TNORM 00515 ELSE 00516 ATOLI = ABSTOL 00517 END IF 00518 * 00519 IF( IRANGE.EQ.2 ) THEN 00520 WL = VL 00521 WU = VU 00522 ELSE 00523 WL = ZERO 00524 WU = ZERO 00525 END IF 00526 END IF 00527 * 00528 * Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. 00529 * NWL accumulates the number of eigenvalues .le. WL, 00530 * NWU accumulates the number of eigenvalues .le. WU 00531 * 00532 M = 0 00533 IEND = 0 00534 INFO = 0 00535 NWL = 0 00536 NWU = 0 00537 * 00538 DO 70 JB = 1, NSPLIT 00539 IOFF = IEND 00540 IBEGIN = IOFF + 1 00541 IEND = ISPLIT( JB ) 00542 IN = IEND - IOFF 00543 * 00544 IF( IN.EQ.1 ) THEN 00545 * 00546 * Special Case -- IN=1 00547 * 00548 IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN ) 00549 $ NWL = NWL + 1 00550 IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN ) 00551 $ NWU = NWU + 1 00552 IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE. 00553 $ D( IBEGIN )-PIVMIN ) ) THEN 00554 M = M + 1 00555 W( M ) = D( IBEGIN ) 00556 IBLOCK( M ) = JB 00557 END IF 00558 ELSE 00559 * 00560 * General Case -- IN > 1 00561 * 00562 * Compute Gershgorin Interval 00563 * and use it as the initial interval 00564 * 00565 GU = D( IBEGIN ) 00566 GL = D( IBEGIN ) 00567 TMP1 = ZERO 00568 * 00569 DO 40 J = IBEGIN, IEND - 1 00570 TMP2 = ABS( E( J ) ) 00571 GU = MAX( GU, D( J )+TMP1+TMP2 ) 00572 GL = MIN( GL, D( J )-TMP1-TMP2 ) 00573 TMP1 = TMP2 00574 40 CONTINUE 00575 * 00576 GU = MAX( GU, D( IEND )+TMP1 ) 00577 GL = MIN( GL, D( IEND )-TMP1 ) 00578 BNORM = MAX( ABS( GL ), ABS( GU ) ) 00579 GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN 00580 GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN 00581 * 00582 * Compute ATOLI for the current submatrix 00583 * 00584 IF( ABSTOL.LE.ZERO ) THEN 00585 ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) ) 00586 ELSE 00587 ATOLI = ABSTOL 00588 END IF 00589 * 00590 IF( IRANGE.GT.1 ) THEN 00591 IF( GU.LT.WL ) THEN 00592 NWL = NWL + IN 00593 NWU = NWU + IN 00594 GO TO 70 00595 END IF 00596 GL = MAX( GL, WL ) 00597 GU = MIN( GU, WU ) 00598 IF( GL.GE.GU ) 00599 $ GO TO 70 00600 END IF 00601 * 00602 * Set Up Initial Interval 00603 * 00604 WORK( N+1 ) = GL 00605 WORK( N+IN+1 ) = GU 00606 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00607 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00608 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 00609 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00610 * 00611 NWL = NWL + IWORK( 1 ) 00612 NWU = NWU + IWORK( IN+1 ) 00613 IWOFF = M - IWORK( 1 ) 00614 * 00615 * Compute Eigenvalues 00616 * 00617 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 00618 $ LOG( TWO ) ) + 2 00619 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 00620 $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), 00621 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 00622 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 00623 * 00624 * Copy Eigenvalues Into W and IBLOCK 00625 * Use -JB for block number for unconverged eigenvalues. 00626 * 00627 DO 60 J = 1, IOUT 00628 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 00629 * 00630 * Flag non-convergence. 00631 * 00632 IF( J.GT.IOUT-IINFO ) THEN 00633 NCNVRG = .TRUE. 00634 IB = -JB 00635 ELSE 00636 IB = JB 00637 END IF 00638 DO 50 JE = IWORK( J ) + 1 + IWOFF, 00639 $ IWORK( J+IN ) + IWOFF 00640 W( JE ) = TMP1 00641 IBLOCK( JE ) = IB 00642 50 CONTINUE 00643 60 CONTINUE 00644 * 00645 M = M + IM 00646 END IF 00647 70 CONTINUE 00648 * 00649 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 00650 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 00651 * 00652 IF( IRANGE.EQ.3 ) THEN 00653 IM = 0 00654 IDISCL = IL - 1 - NWL 00655 IDISCU = NWU - IU 00656 * 00657 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00658 DO 80 JE = 1, M 00659 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 00660 IDISCL = IDISCL - 1 00661 ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 00662 IDISCU = IDISCU - 1 00663 ELSE 00664 IM = IM + 1 00665 W( IM ) = W( JE ) 00666 IBLOCK( IM ) = IBLOCK( JE ) 00667 END IF 00668 80 CONTINUE 00669 M = IM 00670 END IF 00671 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 00672 * 00673 * Code to deal with effects of bad arithmetic: 00674 * Some low eigenvalues to be discarded are not in (WL,WLU], 00675 * or high eigenvalues to be discarded are not in (WUL,WU] 00676 * so just kill off the smallest IDISCL/largest IDISCU 00677 * eigenvalues, by simply finding the smallest/largest 00678 * eigenvalue(s). 00679 * 00680 * (If N(w) is monotone non-decreasing, this should never 00681 * happen.) 00682 * 00683 IF( IDISCL.GT.0 ) THEN 00684 WKILL = WU 00685 DO 100 JDISC = 1, IDISCL 00686 IW = 0 00687 DO 90 JE = 1, M 00688 IF( IBLOCK( JE ).NE.0 .AND. 00689 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 00690 IW = JE 00691 WKILL = W( JE ) 00692 END IF 00693 90 CONTINUE 00694 IBLOCK( IW ) = 0 00695 100 CONTINUE 00696 END IF 00697 IF( IDISCU.GT.0 ) THEN 00698 * 00699 WKILL = WL 00700 DO 120 JDISC = 1, IDISCU 00701 IW = 0 00702 DO 110 JE = 1, M 00703 IF( IBLOCK( JE ).NE.0 .AND. 00704 $ ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN 00705 IW = JE 00706 WKILL = W( JE ) 00707 END IF 00708 110 CONTINUE 00709 IBLOCK( IW ) = 0 00710 120 CONTINUE 00711 END IF 00712 IM = 0 00713 DO 130 JE = 1, M 00714 IF( IBLOCK( JE ).NE.0 ) THEN 00715 IM = IM + 1 00716 W( IM ) = W( JE ) 00717 IBLOCK( IM ) = IBLOCK( JE ) 00718 END IF 00719 130 CONTINUE 00720 M = IM 00721 END IF 00722 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 00723 TOOFEW = .TRUE. 00724 END IF 00725 END IF 00726 * 00727 * If ORDER='B', do nothing -- the eigenvalues are already sorted 00728 * by block. 00729 * If ORDER='E', sort the eigenvalues from smallest to largest 00730 * 00731 IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN 00732 DO 150 JE = 1, M - 1 00733 IE = 0 00734 TMP1 = W( JE ) 00735 DO 140 J = JE + 1, M 00736 IF( W( J ).LT.TMP1 ) THEN 00737 IE = J 00738 TMP1 = W( J ) 00739 END IF 00740 140 CONTINUE 00741 * 00742 IF( IE.NE.0 ) THEN 00743 ITMP1 = IBLOCK( IE ) 00744 W( IE ) = W( JE ) 00745 IBLOCK( IE ) = IBLOCK( JE ) 00746 W( JE ) = TMP1 00747 IBLOCK( JE ) = ITMP1 00748 END IF 00749 150 CONTINUE 00750 END IF 00751 * 00752 INFO = 0 00753 IF( NCNVRG ) 00754 $ INFO = INFO + 1 00755 IF( TOOFEW ) 00756 $ INFO = INFO + 2 00757 RETURN 00758 * 00759 * End of DSTEBZ 00760 * 00761 END