1*> \brief \b DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARRD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS, 22* RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 23* M, W, WERR, WL, WU, IBLOCK, INDEXW, 24* WORK, IWORK, INFO ) 25* 26* .. Scalar Arguments .. 27* CHARACTER ORDER, RANGE 28* INTEGER IL, INFO, IU, M, N, NSPLIT 29* DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU 30* .. 31* .. Array Arguments .. 32* INTEGER IBLOCK( * ), INDEXW( * ), 33* $ ISPLIT( * ), IWORK( * ) 34* DOUBLE PRECISION D( * ), E( * ), E2( * ), 35* $ GERS( * ), W( * ), WERR( * ), WORK( * ) 36* .. 37* 38* 39*> \par Purpose: 40* ============= 41*> 42*> \verbatim 43*> 44*> DLARRD computes the eigenvalues of a symmetric tridiagonal 45*> matrix T to suitable accuracy. This is an auxiliary code to be 46*> called from DSTEMR. 47*> The user may ask for all eigenvalues, all eigenvalues 48*> in the half-open interval (VL, VU], or the IL-th through IU-th 49*> eigenvalues. 50*> 51*> To avoid overflow, the matrix must be scaled so that its 52*> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest 53*> accuracy, it should not be much smaller than that. 54*> 55*> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 56*> Matrix", Report CS41, Computer Science Dept., Stanford 57*> University, July 21, 1966. 58*> \endverbatim 59* 60* Arguments: 61* ========== 62* 63*> \param[in] RANGE 64*> \verbatim 65*> RANGE is CHARACTER*1 66*> = 'A': ("All") all eigenvalues will be found. 67*> = 'V': ("Value") all eigenvalues in the half-open interval 68*> (VL, VU] will be found. 69*> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 70*> entire matrix) will be found. 71*> \endverbatim 72*> 73*> \param[in] ORDER 74*> \verbatim 75*> ORDER is CHARACTER*1 76*> = 'B': ("By Block") the eigenvalues will be grouped by 77*> split-off block (see IBLOCK, ISPLIT) and 78*> ordered from smallest to largest within 79*> the block. 80*> = 'E': ("Entire matrix") 81*> the eigenvalues for the entire matrix 82*> will be ordered from smallest to 83*> largest. 84*> \endverbatim 85*> 86*> \param[in] N 87*> \verbatim 88*> N is INTEGER 89*> The order of the tridiagonal matrix T. N >= 0. 90*> \endverbatim 91*> 92*> \param[in] VL 93*> \verbatim 94*> VL is DOUBLE PRECISION 95*> If RANGE='V', the lower bound of the interval to 96*> be searched for eigenvalues. Eigenvalues less than or equal 97*> to VL, or greater than VU, will not be returned. VL < VU. 98*> Not referenced if RANGE = 'A' or 'I'. 99*> \endverbatim 100*> 101*> \param[in] VU 102*> \verbatim 103*> VU is DOUBLE PRECISION 104*> If RANGE='V', the upper bound of the interval to 105*> be searched for eigenvalues. Eigenvalues less than or equal 106*> to VL, or greater than VU, will not be returned. VL < VU. 107*> Not referenced if RANGE = 'A' or 'I'. 108*> \endverbatim 109*> 110*> \param[in] IL 111*> \verbatim 112*> IL is INTEGER 113*> If RANGE='I', the index of the 114*> smallest eigenvalue to be returned. 115*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 116*> Not referenced if RANGE = 'A' or 'V'. 117*> \endverbatim 118*> 119*> \param[in] IU 120*> \verbatim 121*> IU is INTEGER 122*> If RANGE='I', the index of the 123*> largest eigenvalue to be returned. 124*> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 125*> Not referenced if RANGE = 'A' or 'V'. 126*> \endverbatim 127*> 128*> \param[in] GERS 129*> \verbatim 130*> GERS is DOUBLE PRECISION array, dimension (2*N) 131*> The N Gerschgorin intervals (the i-th Gerschgorin interval 132*> is (GERS(2*i-1), GERS(2*i)). 133*> \endverbatim 134*> 135*> \param[in] RELTOL 136*> \verbatim 137*> RELTOL is DOUBLE PRECISION 138*> The minimum relative width of an interval. When an interval 139*> is narrower than RELTOL times the larger (in 140*> magnitude) endpoint, then it is considered to be 141*> sufficiently small, i.e., converged. Note: this should 142*> always be at least radix*machine epsilon. 143*> \endverbatim 144*> 145*> \param[in] D 146*> \verbatim 147*> D is DOUBLE PRECISION array, dimension (N) 148*> The n diagonal elements of the tridiagonal matrix T. 149*> \endverbatim 150*> 151*> \param[in] E 152*> \verbatim 153*> E is DOUBLE PRECISION array, dimension (N-1) 154*> The (n-1) off-diagonal elements of the tridiagonal matrix T. 155*> \endverbatim 156*> 157*> \param[in] E2 158*> \verbatim 159*> E2 is DOUBLE PRECISION array, dimension (N-1) 160*> The (n-1) squared off-diagonal elements of the tridiagonal matrix T. 161*> \endverbatim 162*> 163*> \param[in] PIVMIN 164*> \verbatim 165*> PIVMIN is DOUBLE PRECISION 166*> The minimum pivot allowed in the Sturm sequence for T. 167*> \endverbatim 168*> 169*> \param[in] NSPLIT 170*> \verbatim 171*> NSPLIT is INTEGER 172*> The number of diagonal blocks in the matrix T. 173*> 1 <= NSPLIT <= N. 174*> \endverbatim 175*> 176*> \param[in] ISPLIT 177*> \verbatim 178*> ISPLIT is INTEGER array, dimension (N) 179*> The splitting points, at which T breaks up into submatrices. 180*> The first submatrix consists of rows/columns 1 to ISPLIT(1), 181*> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 182*> etc., and the NSPLIT-th consists of rows/columns 183*> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 184*> (Only the first NSPLIT elements will actually be used, but 185*> since the user cannot know a priori what value NSPLIT will 186*> have, N words must be reserved for ISPLIT.) 187*> \endverbatim 188*> 189*> \param[out] M 190*> \verbatim 191*> M is INTEGER 192*> The actual number of eigenvalues found. 0 <= M <= N. 193*> (See also the description of INFO=2,3.) 194*> \endverbatim 195*> 196*> \param[out] W 197*> \verbatim 198*> W is DOUBLE PRECISION array, dimension (N) 199*> On exit, the first M elements of W will contain the 200*> eigenvalue approximations. DLARRD computes an interval 201*> I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue 202*> approximation is given as the interval midpoint 203*> W(j)= ( a_j + b_j)/2. The corresponding error is bounded by 204*> WERR(j) = abs( a_j - b_j)/2 205*> \endverbatim 206*> 207*> \param[out] WERR 208*> \verbatim 209*> WERR is DOUBLE PRECISION array, dimension (N) 210*> The error bound on the corresponding eigenvalue approximation 211*> in W. 212*> \endverbatim 213*> 214*> \param[out] WL 215*> \verbatim 216*> WL is DOUBLE PRECISION 217*> \endverbatim 218*> 219*> \param[out] WU 220*> \verbatim 221*> WU is DOUBLE PRECISION 222*> The interval (WL, WU] contains all the wanted eigenvalues. 223*> If RANGE='V', then WL=VL and WU=VU. 224*> If RANGE='A', then WL and WU are the global Gerschgorin bounds 225*> on the spectrum. 226*> If RANGE='I', then WL and WU are computed by DLAEBZ from the 227*> index range specified. 228*> \endverbatim 229*> 230*> \param[out] IBLOCK 231*> \verbatim 232*> IBLOCK is INTEGER array, dimension (N) 233*> At each row/column j where E(j) is zero or small, the 234*> matrix T is considered to split into a block diagonal 235*> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which 236*> block (from 1 to the number of blocks) the eigenvalue W(i) 237*> belongs. (DLARRD may use the remaining N-M elements as 238*> workspace.) 239*> \endverbatim 240*> 241*> \param[out] INDEXW 242*> \verbatim 243*> INDEXW is INTEGER array, dimension (N) 244*> The indices of the eigenvalues within each block (submatrix); 245*> for example, INDEXW(i)= j and IBLOCK(i)=k imply that the 246*> i-th eigenvalue W(i) is the j-th eigenvalue in block k. 247*> \endverbatim 248*> 249*> \param[out] WORK 250*> \verbatim 251*> WORK is DOUBLE PRECISION array, dimension (4*N) 252*> \endverbatim 253*> 254*> \param[out] IWORK 255*> \verbatim 256*> IWORK is INTEGER array, dimension (3*N) 257*> \endverbatim 258*> 259*> \param[out] INFO 260*> \verbatim 261*> INFO is INTEGER 262*> = 0: successful exit 263*> < 0: if INFO = -i, the i-th argument had an illegal value 264*> > 0: some or all of the eigenvalues failed to converge or 265*> were not computed: 266*> =1 or 3: Bisection failed to converge for some 267*> eigenvalues; these eigenvalues are flagged by a 268*> negative block number. The effect is that the 269*> eigenvalues may not be as accurate as the 270*> absolute and relative tolerances. This is 271*> generally caused by unexpectedly inaccurate 272*> arithmetic. 273*> =2 or 3: RANGE='I' only: Not all of the eigenvalues 274*> IL:IU were found. 275*> Effect: M < IU+1-IL 276*> Cause: non-monotonic arithmetic, causing the 277*> Sturm sequence to be non-monotonic. 278*> Cure: recalculate, using RANGE='A', and pick 279*> out eigenvalues IL:IU. In some cases, 280*> increasing the PARAMETER "FUDGE" may 281*> make things work. 282*> = 4: RANGE='I', and the Gershgorin interval 283*> initially used was too small. No eigenvalues 284*> were computed. 285*> Probable cause: your machine has sloppy 286*> floating-point arithmetic. 287*> Cure: Increase the PARAMETER "FUDGE", 288*> recompile, and try again. 289*> \endverbatim 290* 291*> \par Internal Parameters: 292* ========================= 293*> 294*> \verbatim 295*> FUDGE DOUBLE PRECISION, default = 2 296*> A "fudge factor" to widen the Gershgorin intervals. Ideally, 297*> a value of 1 should work, but on machines with sloppy 298*> arithmetic, this needs to be larger. The default for 299*> publicly released versions should be large enough to handle 300*> the worst machine around. Note that this has no effect 301*> on accuracy of the solution. 302*> \endverbatim 303*> 304*> \par Contributors: 305* ================== 306*> 307*> W. Kahan, University of California, Berkeley, USA \n 308*> Beresford Parlett, University of California, Berkeley, USA \n 309*> Jim Demmel, University of California, Berkeley, USA \n 310*> Inderjit Dhillon, University of Texas, Austin, USA \n 311*> Osni Marques, LBNL/NERSC, USA \n 312*> Christof Voemel, University of California, Berkeley, USA \n 313* 314* Authors: 315* ======== 316* 317*> \author Univ. of Tennessee 318*> \author Univ. of California Berkeley 319*> \author Univ. of Colorado Denver 320*> \author NAG Ltd. 321* 322*> \ingroup OTHERauxiliary 323* 324* ===================================================================== 325 SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS, 326 $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 327 $ M, W, WERR, WL, WU, IBLOCK, INDEXW, 328 $ WORK, IWORK, INFO ) 329* 330* -- LAPACK auxiliary routine -- 331* -- LAPACK is a software package provided by Univ. of Tennessee, -- 332* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 333* 334* .. Scalar Arguments .. 335 CHARACTER ORDER, RANGE 336 INTEGER IL, INFO, IU, M, N, NSPLIT 337 DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU 338* .. 339* .. Array Arguments .. 340 INTEGER IBLOCK( * ), INDEXW( * ), 341 $ ISPLIT( * ), IWORK( * ) 342 DOUBLE PRECISION D( * ), E( * ), E2( * ), 343 $ GERS( * ), W( * ), WERR( * ), WORK( * ) 344* .. 345* 346* ===================================================================== 347* 348* .. Parameters .. 349 DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE 350 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 351 $ TWO = 2.0D0, HALF = ONE/TWO, 352 $ FUDGE = TWO ) 353 INTEGER ALLRNG, VALRNG, INDRNG 354 PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 ) 355* .. 356* .. Local Scalars .. 357 LOGICAL NCNVRG, TOOFEW 358 INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, 359 $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1, 360 $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB, 361 $ NWL, NWU 362 DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2, 363 $ TNORM, UFLOW, WKILL, WLU, WUL 364 365* .. 366* .. Local Arrays .. 367 INTEGER IDUMMA( 1 ) 368* .. 369* .. External Functions .. 370 LOGICAL LSAME 371 INTEGER ILAENV 372 DOUBLE PRECISION DLAMCH 373 EXTERNAL LSAME, ILAENV, DLAMCH 374* .. 375* .. External Subroutines .. 376 EXTERNAL DLAEBZ 377* .. 378* .. Intrinsic Functions .. 379 INTRINSIC ABS, INT, LOG, MAX, MIN 380* .. 381* .. Executable Statements .. 382* 383 INFO = 0 384* 385* Quick return if possible 386* 387 IF( N.LE.0 ) THEN 388 RETURN 389 END IF 390* 391* Decode RANGE 392* 393 IF( LSAME( RANGE, 'A' ) ) THEN 394 IRANGE = ALLRNG 395 ELSE IF( LSAME( RANGE, 'V' ) ) THEN 396 IRANGE = VALRNG 397 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 398 IRANGE = INDRNG 399 ELSE 400 IRANGE = 0 401 END IF 402* 403* Check for Errors 404* 405 IF( IRANGE.LE.0 ) THEN 406 INFO = -1 407 ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN 408 INFO = -2 409 ELSE IF( N.LT.0 ) THEN 410 INFO = -3 411 ELSE IF( IRANGE.EQ.VALRNG ) THEN 412 IF( VL.GE.VU ) 413 $ INFO = -5 414 ELSE IF( IRANGE.EQ.INDRNG .AND. 415 $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN 416 INFO = -6 417 ELSE IF( IRANGE.EQ.INDRNG .AND. 418 $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 419 INFO = -7 420 END IF 421* 422 IF( INFO.NE.0 ) THEN 423 RETURN 424 END IF 425 426* Initialize error flags 427 INFO = 0 428 NCNVRG = .FALSE. 429 TOOFEW = .FALSE. 430 431* Quick return if possible 432 M = 0 433 IF( N.EQ.0 ) RETURN 434 435* Simplification: 436 IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1 437 438* Get machine constants 439 EPS = DLAMCH( 'P' ) 440 UFLOW = DLAMCH( 'U' ) 441 442 443* Special Case when N=1 444* Treat case of 1x1 matrix for quick return 445 IF( N.EQ.1 ) THEN 446 IF( (IRANGE.EQ.ALLRNG).OR. 447 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. 448 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN 449 M = 1 450 W(1) = D(1) 451* The computation error of the eigenvalue is zero 452 WERR(1) = ZERO 453 IBLOCK( 1 ) = 1 454 INDEXW( 1 ) = 1 455 ENDIF 456 RETURN 457 END IF 458 459* NB is the minimum vector length for vector bisection, or 0 460* if only scalar is to be done. 461 NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) 462 IF( NB.LE.1 ) NB = 0 463 464* Find global spectral radius 465 GL = D(1) 466 GU = D(1) 467 DO 5 I = 1,N 468 GL = MIN( GL, GERS( 2*I - 1)) 469 GU = MAX( GU, GERS(2*I) ) 470 5 CONTINUE 471* Compute global Gerschgorin bounds and spectral diameter 472 TNORM = MAX( ABS( GL ), ABS( GU ) ) 473 GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN 474 GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN 475* [JAN/28/2009] remove the line below since SPDIAM variable not use 476* SPDIAM = GU - GL 477* Input arguments for DLAEBZ: 478* The relative tolerance. An interval (a,b] lies within 479* "relative tolerance" if b-a < RELTOL*max(|a|,|b|), 480 RTOLI = RELTOL 481* Set the absolute tolerance for interval convergence to zero to force 482* interval convergence based on relative size of the interval. 483* This is dangerous because intervals might not converge when RELTOL is 484* small. But at least a very small number should be selected so that for 485* strongly graded matrices, the code can get relatively accurate 486* eigenvalues. 487 ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN 488 489 IF( IRANGE.EQ.INDRNG ) THEN 490 491* RANGE='I': Compute an interval containing eigenvalues 492* IL through IU. The initial interval [GL,GU] from the global 493* Gerschgorin bounds GL and GU is refined by DLAEBZ. 494 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 495 $ LOG( TWO ) ) + 2 496 WORK( N+1 ) = GL 497 WORK( N+2 ) = GL 498 WORK( N+3 ) = GU 499 WORK( N+4 ) = GU 500 WORK( N+5 ) = GL 501 WORK( N+6 ) = GU 502 IWORK( 1 ) = -1 503 IWORK( 2 ) = -1 504 IWORK( 3 ) = N + 1 505 IWORK( 4 ) = N + 1 506 IWORK( 5 ) = IL - 1 507 IWORK( 6 ) = IU 508* 509 CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, 510 $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, 511 $ IWORK, W, IBLOCK, IINFO ) 512 IF( IINFO .NE. 0 ) THEN 513 INFO = IINFO 514 RETURN 515 END IF 516* On exit, output intervals may not be ordered by ascending negcount 517 IF( IWORK( 6 ).EQ.IU ) THEN 518 WL = WORK( N+1 ) 519 WLU = WORK( N+3 ) 520 NWL = IWORK( 1 ) 521 WU = WORK( N+4 ) 522 WUL = WORK( N+2 ) 523 NWU = IWORK( 4 ) 524 ELSE 525 WL = WORK( N+2 ) 526 WLU = WORK( N+4 ) 527 NWL = IWORK( 2 ) 528 WU = WORK( N+3 ) 529 WUL = WORK( N+1 ) 530 NWU = IWORK( 3 ) 531 END IF 532* On exit, the interval [WL, WLU] contains a value with negcount NWL, 533* and [WUL, WU] contains a value with negcount NWU. 534 IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN 535 INFO = 4 536 RETURN 537 END IF 538 539 ELSEIF( IRANGE.EQ.VALRNG ) THEN 540 WL = VL 541 WU = VU 542 543 ELSEIF( IRANGE.EQ.ALLRNG ) THEN 544 WL = GL 545 WU = GU 546 ENDIF 547 548 549 550* Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. 551* NWL accumulates the number of eigenvalues .le. WL, 552* NWU accumulates the number of eigenvalues .le. WU 553 M = 0 554 IEND = 0 555 INFO = 0 556 NWL = 0 557 NWU = 0 558* 559 DO 70 JBLK = 1, NSPLIT 560 IOFF = IEND 561 IBEGIN = IOFF + 1 562 IEND = ISPLIT( JBLK ) 563 IN = IEND - IOFF 564* 565 IF( IN.EQ.1 ) THEN 566* 1x1 block 567 IF( WL.GE.D( IBEGIN )-PIVMIN ) 568 $ NWL = NWL + 1 569 IF( WU.GE.D( IBEGIN )-PIVMIN ) 570 $ NWU = NWU + 1 571 IF( IRANGE.EQ.ALLRNG .OR. 572 $ ( WL.LT.D( IBEGIN )-PIVMIN 573 $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN 574 M = M + 1 575 W( M ) = D( IBEGIN ) 576 WERR(M) = ZERO 577* The gap for a single block doesn't matter for the later 578* algorithm and is assigned an arbitrary large value 579 IBLOCK( M ) = JBLK 580 INDEXW( M ) = 1 581 END IF 582 583* Disabled 2x2 case because of a failure on the following matrix 584* RANGE = 'I', IL = IU = 4 585* Original Tridiagonal, d = [ 586* -0.150102010615740E+00 587* -0.849897989384260E+00 588* -0.128208148052635E-15 589* 0.128257718286320E-15 590* ]; 591* e = [ 592* -0.357171383266986E+00 593* -0.180411241501588E-15 594* -0.175152352710251E-15 595* ]; 596* 597* ELSE IF( IN.EQ.2 ) THEN 598** 2x2 block 599* DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) 600* TMP1 = HALF*(D(IBEGIN)+D(IEND)) 601* L1 = TMP1 - DISC 602* IF( WL.GE. L1-PIVMIN ) 603* $ NWL = NWL + 1 604* IF( WU.GE. L1-PIVMIN ) 605* $ NWU = NWU + 1 606* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. 607* $ L1-PIVMIN ) ) THEN 608* M = M + 1 609* W( M ) = L1 610** The uncertainty of eigenvalues of a 2x2 matrix is very small 611* WERR( M ) = EPS * ABS( W( M ) ) * TWO 612* IBLOCK( M ) = JBLK 613* INDEXW( M ) = 1 614* ENDIF 615* L2 = TMP1 + DISC 616* IF( WL.GE. L2-PIVMIN ) 617* $ NWL = NWL + 1 618* IF( WU.GE. L2-PIVMIN ) 619* $ NWU = NWU + 1 620* IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. 621* $ L2-PIVMIN ) ) THEN 622* M = M + 1 623* W( M ) = L2 624** The uncertainty of eigenvalues of a 2x2 matrix is very small 625* WERR( M ) = EPS * ABS( W( M ) ) * TWO 626* IBLOCK( M ) = JBLK 627* INDEXW( M ) = 2 628* ENDIF 629 ELSE 630* General Case - block of size IN >= 2 631* Compute local Gerschgorin interval and use it as the initial 632* interval for DLAEBZ 633 GU = D( IBEGIN ) 634 GL = D( IBEGIN ) 635 TMP1 = ZERO 636 637 DO 40 J = IBEGIN, IEND 638 GL = MIN( GL, GERS( 2*J - 1)) 639 GU = MAX( GU, GERS(2*J) ) 640 40 CONTINUE 641* [JAN/28/2009] 642* change SPDIAM by TNORM in lines 2 and 3 thereafter 643* line 1: remove computation of SPDIAM (not useful anymore) 644* SPDIAM = GU - GL 645* GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN 646* GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN 647 GL = GL - FUDGE*TNORM*EPS*IN - FUDGE*PIVMIN 648 GU = GU + FUDGE*TNORM*EPS*IN + FUDGE*PIVMIN 649* 650 IF( IRANGE.GT.1 ) THEN 651 IF( GU.LT.WL ) THEN 652* the local block contains none of the wanted eigenvalues 653 NWL = NWL + IN 654 NWU = NWU + IN 655 GO TO 70 656 END IF 657* refine search interval if possible, only range (WL,WU] matters 658 GL = MAX( GL, WL ) 659 GU = MIN( GU, WU ) 660 IF( GL.GE.GU ) 661 $ GO TO 70 662 END IF 663 664* Find negcount of initial interval boundaries GL and GU 665 WORK( N+1 ) = GL 666 WORK( N+IN+1 ) = GU 667 CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 668 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), 669 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, 670 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 671 IF( IINFO .NE. 0 ) THEN 672 INFO = IINFO 673 RETURN 674 END IF 675* 676 NWL = NWL + IWORK( 1 ) 677 NWU = NWU + IWORK( IN+1 ) 678 IWOFF = M - IWORK( 1 ) 679 680* Compute Eigenvalues 681 ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / 682 $ LOG( TWO ) ) + 2 683 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, 684 $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), 685 $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, 686 $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) 687 IF( IINFO .NE. 0 ) THEN 688 INFO = IINFO 689 RETURN 690 END IF 691* 692* Copy eigenvalues into W and IBLOCK 693* Use -JBLK for block number for unconverged eigenvalues. 694* Loop over the number of output intervals from DLAEBZ 695 DO 60 J = 1, IOUT 696* eigenvalue approximation is middle point of interval 697 TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) 698* semi length of error interval 699 TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) ) 700 IF( J.GT.IOUT-IINFO ) THEN 701* Flag non-convergence. 702 NCNVRG = .TRUE. 703 IB = -JBLK 704 ELSE 705 IB = JBLK 706 END IF 707 DO 50 JE = IWORK( J ) + 1 + IWOFF, 708 $ IWORK( J+IN ) + IWOFF 709 W( JE ) = TMP1 710 WERR( JE ) = TMP2 711 INDEXW( JE ) = JE - IWOFF 712 IBLOCK( JE ) = IB 713 50 CONTINUE 714 60 CONTINUE 715* 716 M = M + IM 717 END IF 718 70 CONTINUE 719 720* If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU 721* If NWL+1 < IL or NWU > IU, discard extra eigenvalues. 722 IF( IRANGE.EQ.INDRNG ) THEN 723 IDISCL = IL - 1 - NWL 724 IDISCU = NWU - IU 725* 726 IF( IDISCL.GT.0 ) THEN 727 IM = 0 728 DO 80 JE = 1, M 729* Remove some of the smallest eigenvalues from the left so that 730* at the end IDISCL =0. Move all eigenvalues up to the left. 731 IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN 732 IDISCL = IDISCL - 1 733 ELSE 734 IM = IM + 1 735 W( IM ) = W( JE ) 736 WERR( IM ) = WERR( JE ) 737 INDEXW( IM ) = INDEXW( JE ) 738 IBLOCK( IM ) = IBLOCK( JE ) 739 END IF 740 80 CONTINUE 741 M = IM 742 END IF 743 IF( IDISCU.GT.0 ) THEN 744* Remove some of the largest eigenvalues from the right so that 745* at the end IDISCU =0. Move all eigenvalues up to the left. 746 IM=M+1 747 DO 81 JE = M, 1, -1 748 IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN 749 IDISCU = IDISCU - 1 750 ELSE 751 IM = IM - 1 752 W( IM ) = W( JE ) 753 WERR( IM ) = WERR( JE ) 754 INDEXW( IM ) = INDEXW( JE ) 755 IBLOCK( IM ) = IBLOCK( JE ) 756 END IF 757 81 CONTINUE 758 JEE = 0 759 DO 82 JE = IM, M 760 JEE = JEE + 1 761 W( JEE ) = W( JE ) 762 WERR( JEE ) = WERR( JE ) 763 INDEXW( JEE ) = INDEXW( JE ) 764 IBLOCK( JEE ) = IBLOCK( JE ) 765 82 CONTINUE 766 M = M-IM+1 767 END IF 768 769 IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN 770* Code to deal with effects of bad arithmetic. (If N(w) is 771* monotone non-decreasing, this should never happen.) 772* Some low eigenvalues to be discarded are not in (WL,WLU], 773* or high eigenvalues to be discarded are not in (WUL,WU] 774* so just kill off the smallest IDISCL/largest IDISCU 775* eigenvalues, by marking the corresponding IBLOCK = 0 776 IF( IDISCL.GT.0 ) THEN 777 WKILL = WU 778 DO 100 JDISC = 1, IDISCL 779 IW = 0 780 DO 90 JE = 1, M 781 IF( IBLOCK( JE ).NE.0 .AND. 782 $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN 783 IW = JE 784 WKILL = W( JE ) 785 END IF 786 90 CONTINUE 787 IBLOCK( IW ) = 0 788 100 CONTINUE 789 END IF 790 IF( IDISCU.GT.0 ) THEN 791 WKILL = WL 792 DO 120 JDISC = 1, IDISCU 793 IW = 0 794 DO 110 JE = 1, M 795 IF( IBLOCK( JE ).NE.0 .AND. 796 $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN 797 IW = JE 798 WKILL = W( JE ) 799 END IF 800 110 CONTINUE 801 IBLOCK( IW ) = 0 802 120 CONTINUE 803 END IF 804* Now erase all eigenvalues with IBLOCK set to zero 805 IM = 0 806 DO 130 JE = 1, M 807 IF( IBLOCK( JE ).NE.0 ) THEN 808 IM = IM + 1 809 W( IM ) = W( JE ) 810 WERR( IM ) = WERR( JE ) 811 INDEXW( IM ) = INDEXW( JE ) 812 IBLOCK( IM ) = IBLOCK( JE ) 813 END IF 814 130 CONTINUE 815 M = IM 816 END IF 817 IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN 818 TOOFEW = .TRUE. 819 END IF 820 END IF 821* 822 IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR. 823 $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN 824 TOOFEW = .TRUE. 825 END IF 826 827* If ORDER='B', do nothing the eigenvalues are already sorted by 828* block. 829* If ORDER='E', sort the eigenvalues from smallest to largest 830 831 IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN 832 DO 150 JE = 1, M - 1 833 IE = 0 834 TMP1 = W( JE ) 835 DO 140 J = JE + 1, M 836 IF( W( J ).LT.TMP1 ) THEN 837 IE = J 838 TMP1 = W( J ) 839 END IF 840 140 CONTINUE 841 IF( IE.NE.0 ) THEN 842 TMP2 = WERR( IE ) 843 ITMP1 = IBLOCK( IE ) 844 ITMP2 = INDEXW( IE ) 845 W( IE ) = W( JE ) 846 WERR( IE ) = WERR( JE ) 847 IBLOCK( IE ) = IBLOCK( JE ) 848 INDEXW( IE ) = INDEXW( JE ) 849 W( JE ) = TMP1 850 WERR( JE ) = TMP2 851 IBLOCK( JE ) = ITMP1 852 INDEXW( JE ) = ITMP2 853 END IF 854 150 CONTINUE 855 END IF 856* 857 INFO = 0 858 IF( NCNVRG ) 859 $ INFO = INFO + 1 860 IF( TOOFEW ) 861 $ INFO = INFO + 2 862 RETURN 863* 864* End of DLARRD 865* 866 END 867