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