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