1*> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARRE + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2, 22* RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, 23* W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, 24* WORK, IWORK, INFO ) 25* 26* .. Scalar Arguments .. 27* CHARACTER RANGE 28* INTEGER IL, INFO, IU, M, N, NSPLIT 29* DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU 30* .. 31* .. Array Arguments .. 32* INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), 33* $ INDEXW( * ) 34* DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ), 35* $ W( * ),WERR( * ), WGAP( * ), WORK( * ) 36* .. 37* 38* 39*> \par Purpose: 40* ============= 41*> 42*> \verbatim 43*> 44*> To find the desired eigenvalues of a given real symmetric 45*> tridiagonal matrix T, DLARRE sets any "small" off-diagonal 46*> elements to zero, and for each unreduced block T_i, it finds 47*> (a) a suitable shift at one end of the block's spectrum, 48*> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and 49*> (c) eigenvalues of each L_i D_i L_i^T. 50*> The representations and eigenvalues found are then used by 51*> DSTEMR to compute the eigenvectors of T. 52*> The accuracy varies depending on whether bisection is used to 53*> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to 54*> conpute all and then discard any unwanted one. 55*> As an added benefit, DLARRE also outputs the n 56*> Gerschgorin intervals for the matrices L_i D_i L_i^T. 57*> \endverbatim 58* 59* Arguments: 60* ========== 61* 62*> \param[in] RANGE 63*> \verbatim 64*> RANGE is CHARACTER*1 65*> = 'A': ("All") all eigenvalues will be found. 66*> = 'V': ("Value") all eigenvalues in the half-open interval 67*> (VL, VU] will be found. 68*> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the 69*> entire matrix) will be found. 70*> \endverbatim 71*> 72*> \param[in] N 73*> \verbatim 74*> N is INTEGER 75*> The order of the matrix. N > 0. 76*> \endverbatim 77*> 78*> \param[in,out] VL 79*> \verbatim 80*> VL is DOUBLE PRECISION 81*> \endverbatim 82*> 83*> \param[in,out] VU 84*> \verbatim 85*> VU is DOUBLE PRECISION 86*> If RANGE='V', the lower and upper bounds for the eigenvalues. 87*> Eigenvalues less than or equal to VL, or greater than VU, 88*> will not be returned. VL < VU. 89*> If RANGE='I' or ='A', DLARRE computes bounds on the desired 90*> part of the spectrum. 91*> \endverbatim 92*> 93*> \param[in] IL 94*> \verbatim 95*> IL is INTEGER 96*> \endverbatim 97*> 98*> \param[in] IU 99*> \verbatim 100*> IU is INTEGER 101*> If RANGE='I', the indices (in ascending order) of the 102*> smallest and largest eigenvalues to be returned. 103*> 1 <= IL <= IU <= N. 104*> \endverbatim 105*> 106*> \param[in,out] D 107*> \verbatim 108*> D is DOUBLE PRECISION array, dimension (N) 109*> On entry, the N diagonal elements of the tridiagonal 110*> matrix T. 111*> On exit, the N diagonal elements of the diagonal 112*> matrices D_i. 113*> \endverbatim 114*> 115*> \param[in,out] E 116*> \verbatim 117*> E is DOUBLE PRECISION array, dimension (N) 118*> On entry, the first (N-1) entries contain the subdiagonal 119*> elements of the tridiagonal matrix T; E(N) need not be set. 120*> On exit, E contains the subdiagonal elements of the unit 121*> bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), 122*> 1 <= I <= NSPLIT, contain the base points sigma_i on output. 123*> \endverbatim 124*> 125*> \param[in,out] E2 126*> \verbatim 127*> E2 is DOUBLE PRECISION array, dimension (N) 128*> On entry, the first (N-1) entries contain the SQUARES of the 129*> subdiagonal elements of the tridiagonal matrix T; 130*> E2(N) need not be set. 131*> On exit, the entries E2( ISPLIT( I ) ), 132*> 1 <= I <= NSPLIT, have been set to zero 133*> \endverbatim 134*> 135*> \param[in] RTOL1 136*> \verbatim 137*> RTOL1 is DOUBLE PRECISION 138*> \endverbatim 139*> 140*> \param[in] RTOL2 141*> \verbatim 142*> RTOL2 is DOUBLE PRECISION 143*> Parameters for bisection. 144*> An interval [LEFT,RIGHT] has converged if 145*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) 146*> \endverbatim 147*> 148*> \param[in] SPLTOL 149*> \verbatim 150*> SPLTOL is DOUBLE PRECISION 151*> The threshold for splitting. 152*> \endverbatim 153*> 154*> \param[out] NSPLIT 155*> \verbatim 156*> NSPLIT is INTEGER 157*> The number of blocks T splits into. 1 <= NSPLIT <= N. 158*> \endverbatim 159*> 160*> \param[out] ISPLIT 161*> \verbatim 162*> ISPLIT is INTEGER array, dimension (N) 163*> The splitting points, at which T breaks up into blocks. 164*> The first block consists of rows/columns 1 to ISPLIT(1), 165*> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), 166*> etc., and the NSPLIT-th consists of rows/columns 167*> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. 168*> \endverbatim 169*> 170*> \param[out] M 171*> \verbatim 172*> M is INTEGER 173*> The total number of eigenvalues (of all L_i D_i L_i^T) 174*> found. 175*> \endverbatim 176*> 177*> \param[out] W 178*> \verbatim 179*> W is DOUBLE PRECISION array, dimension (N) 180*> The first M elements contain the eigenvalues. The 181*> eigenvalues of each of the blocks, L_i D_i L_i^T, are 182*> sorted in ascending order ( DLARRE may use the 183*> remaining N-M elements as workspace). 184*> \endverbatim 185*> 186*> \param[out] WERR 187*> \verbatim 188*> WERR is DOUBLE PRECISION array, dimension (N) 189*> The error bound on the corresponding eigenvalue in W. 190*> \endverbatim 191*> 192*> \param[out] WGAP 193*> \verbatim 194*> WGAP is DOUBLE PRECISION array, dimension (N) 195*> The separation from the right neighbor eigenvalue in W. 196*> The gap is only with respect to the eigenvalues of the same block 197*> as each block has its own representation tree. 198*> Exception: at the right end of a block we store the left gap 199*> \endverbatim 200*> 201*> \param[out] IBLOCK 202*> \verbatim 203*> IBLOCK is INTEGER array, dimension (N) 204*> The indices of the blocks (submatrices) associated with the 205*> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue 206*> W(i) belongs to the first block from the top, =2 if W(i) 207*> belongs to the second block, etc. 208*> \endverbatim 209*> 210*> \param[out] INDEXW 211*> \verbatim 212*> INDEXW is INTEGER array, dimension (N) 213*> The indices of the eigenvalues within each block (submatrix); 214*> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the 215*> i-th eigenvalue W(i) is the 10-th eigenvalue in block 2 216*> \endverbatim 217*> 218*> \param[out] GERS 219*> \verbatim 220*> GERS is DOUBLE PRECISION array, dimension (2*N) 221*> The N Gerschgorin intervals (the i-th Gerschgorin interval 222*> is (GERS(2*i-1), GERS(2*i)). 223*> \endverbatim 224*> 225*> \param[out] PIVMIN 226*> \verbatim 227*> PIVMIN is DOUBLE PRECISION 228*> The minimum pivot in the Sturm sequence for T. 229*> \endverbatim 230*> 231*> \param[out] WORK 232*> \verbatim 233*> WORK is DOUBLE PRECISION array, dimension (6*N) 234*> Workspace. 235*> \endverbatim 236*> 237*> \param[out] IWORK 238*> \verbatim 239*> IWORK is INTEGER array, dimension (5*N) 240*> Workspace. 241*> \endverbatim 242*> 243*> \param[out] INFO 244*> \verbatim 245*> INFO is INTEGER 246*> = 0: successful exit 247*> > 0: A problem occured in DLARRE. 248*> < 0: One of the called subroutines signaled an internal problem. 249*> Needs inspection of the corresponding parameter IINFO 250*> for further information. 251*> 252*> =-1: Problem in DLARRD. 253*> = 2: No base representation could be found in MAXTRY iterations. 254*> Increasing MAXTRY and recompilation might be a remedy. 255*> =-3: Problem in DLARRB when computing the refined root 256*> representation for DLASQ2. 257*> =-4: Problem in DLARRB when preforming bisection on the 258*> desired part of the spectrum. 259*> =-5: Problem in DLASQ2. 260*> =-6: Problem in DLASQ2. 261*> \endverbatim 262* 263* Authors: 264* ======== 265* 266*> \author Univ. of Tennessee 267*> \author Univ. of California Berkeley 268*> \author Univ. of Colorado Denver 269*> \author NAG Ltd. 270* 271*> \date September 2012 272* 273*> \ingroup auxOTHERauxiliary 274* 275*> \par Further Details: 276* ===================== 277*> 278*> \verbatim 279*> 280*> The base representations are required to suffer very little 281*> element growth and consequently define all their eigenvalues to 282*> high relative accuracy. 283*> \endverbatim 284* 285*> \par Contributors: 286* ================== 287*> 288*> Beresford Parlett, University of California, Berkeley, USA \n 289*> Jim Demmel, University of California, Berkeley, USA \n 290*> Inderjit Dhillon, University of Texas, Austin, USA \n 291*> Osni Marques, LBNL/NERSC, USA \n 292*> Christof Voemel, University of California, Berkeley, USA \n 293*> 294* ===================================================================== 295 SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2, 296 $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, 297 $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, 298 $ WORK, IWORK, INFO ) 299* 300* -- LAPACK auxiliary routine (version 3.4.2) -- 301* -- LAPACK is a software package provided by Univ. of Tennessee, -- 302* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 303* September 2012 304* 305* .. Scalar Arguments .. 306 CHARACTER RANGE 307 INTEGER IL, INFO, IU, M, N, NSPLIT 308 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU 309* .. 310* .. Array Arguments .. 311 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), 312 $ INDEXW( * ) 313 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ), 314 $ W( * ),WERR( * ), WGAP( * ), WORK( * ) 315* .. 316* 317* ===================================================================== 318* 319* .. Parameters .. 320 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD, 321 $ MAXGROWTH, ONE, PERT, TWO, ZERO 322 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 323 $ TWO = 2.0D0, FOUR=4.0D0, 324 $ HNDRD = 100.0D0, 325 $ PERT = 8.0D0, 326 $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF, 327 $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 ) 328 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG 329 PARAMETER ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2, 330 $ VALRNG = 3 ) 331* .. 332* .. Local Scalars .. 333 LOGICAL FORCEB, NOREP, USEDQD 334 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO, 335 $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM, 336 $ WBEGIN, WEND 337 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS, 338 $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL, 339 $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM, 340 $ TAU, TMP, TMP1 341 342 343* .. 344* .. Local Arrays .. 345 INTEGER ISEED( 4 ) 346* .. 347* .. External Functions .. 348 LOGICAL LSAME 349 DOUBLE PRECISION DLAMCH 350 EXTERNAL DLAMCH, LSAME 351 352* .. 353* .. External Subroutines .. 354 EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB, DLARRC, DLARRD, 355 $ DLASQ2 356* .. 357* .. Intrinsic Functions .. 358 INTRINSIC ABS, MAX, MIN 359 360* .. 361* .. Executable Statements .. 362* 363 364 INFO = 0 365 366* 367* Decode RANGE 368* 369 IF( LSAME( RANGE, 'A' ) ) THEN 370 IRANGE = ALLRNG 371 ELSE IF( LSAME( RANGE, 'V' ) ) THEN 372 IRANGE = VALRNG 373 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 374 IRANGE = INDRNG 375 END IF 376 377 M = 0 378 379* Get machine constants 380 SAFMIN = DLAMCH( 'S' ) 381 EPS = DLAMCH( 'P' ) 382 383* Set parameters 384 RTL = SQRT(EPS) 385 BSRTOL = SQRT(EPS) 386 387* Treat case of 1x1 matrix for quick return 388 IF( N.EQ.1 ) THEN 389 IF( (IRANGE.EQ.ALLRNG).OR. 390 $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. 391 $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN 392 M = 1 393 W(1) = D(1) 394* The computation error of the eigenvalue is zero 395 WERR(1) = ZERO 396 WGAP(1) = ZERO 397 IBLOCK( 1 ) = 1 398 INDEXW( 1 ) = 1 399 GERS(1) = D( 1 ) 400 GERS(2) = D( 1 ) 401 ENDIF 402* store the shift for the initial RRR, which is zero in this case 403 E(1) = ZERO 404 RETURN 405 END IF 406 407* General case: tridiagonal matrix of order > 1 408* 409* Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter. 410* Compute maximum off-diagonal entry and pivmin. 411 GL = D(1) 412 GU = D(1) 413 EOLD = ZERO 414 EMAX = ZERO 415 E(N) = ZERO 416 DO 5 I = 1,N 417 WERR(I) = ZERO 418 WGAP(I) = ZERO 419 EABS = ABS( E(I) ) 420 IF( EABS .GE. EMAX ) THEN 421 EMAX = EABS 422 END IF 423 TMP1 = EABS + EOLD 424 GERS( 2*I-1) = D(I) - TMP1 425 GL = MIN( GL, GERS( 2*I - 1)) 426 GERS( 2*I ) = D(I) + TMP1 427 GU = MAX( GU, GERS(2*I) ) 428 EOLD = EABS 429 5 CONTINUE 430* The minimum pivot allowed in the Sturm sequence for T 431 PIVMIN = SAFMIN * MAX( ONE, EMAX**2 ) 432* Compute spectral diameter. The Gerschgorin bounds give an 433* estimate that is wrong by at most a factor of SQRT(2) 434 SPDIAM = GU - GL 435 436* Compute splitting points 437 CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM, 438 $ NSPLIT, ISPLIT, IINFO ) 439 440* Can force use of bisection instead of faster DQDS. 441* Option left in the code for future multisection work. 442 FORCEB = .FALSE. 443 444* Initialize USEDQD, DQDS should be used for ALLRNG unless someone 445* explicitly wants bisection. 446 USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB)) 447 448 IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN 449* Set interval [VL,VU] that contains all eigenvalues 450 VL = GL 451 VU = GU 452 ELSE 453* We call DLARRD to find crude approximations to the eigenvalues 454* in the desired range. In case IRANGE = INDRNG, we also obtain the 455* interval (VL,VU] that contains all the wanted eigenvalues. 456* An interval [LEFT,RIGHT] has converged if 457* RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT)) 458* DLARRD needs a WORK of size 4*N, IWORK of size 3*N 459 CALL DLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS, 460 $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, 461 $ MM, W, WERR, VL, VU, IBLOCK, INDEXW, 462 $ WORK, IWORK, IINFO ) 463 IF( IINFO.NE.0 ) THEN 464 INFO = -1 465 RETURN 466 ENDIF 467* Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0 468 DO 14 I = MM+1,N 469 W( I ) = ZERO 470 WERR( I ) = ZERO 471 IBLOCK( I ) = 0 472 INDEXW( I ) = 0 473 14 CONTINUE 474 END IF 475 476 477*** 478* Loop over unreduced blocks 479 IBEGIN = 1 480 WBEGIN = 1 481 DO 170 JBLK = 1, NSPLIT 482 IEND = ISPLIT( JBLK ) 483 IN = IEND - IBEGIN + 1 484 485* 1 X 1 block 486 IF( IN.EQ.1 ) THEN 487 IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND. 488 $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) ) 489 $ .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK)) 490 $ ) THEN 491 M = M + 1 492 W( M ) = D( IBEGIN ) 493 WERR(M) = ZERO 494* The gap for a single block doesn't matter for the later 495* algorithm and is assigned an arbitrary large value 496 WGAP(M) = ZERO 497 IBLOCK( M ) = JBLK 498 INDEXW( M ) = 1 499 WBEGIN = WBEGIN + 1 500 ENDIF 501* E( IEND ) holds the shift for the initial RRR 502 E( IEND ) = ZERO 503 IBEGIN = IEND + 1 504 GO TO 170 505 END IF 506* 507* Blocks of size larger than 1x1 508* 509* E( IEND ) will hold the shift for the initial RRR, for now set it =0 510 E( IEND ) = ZERO 511* 512* Find local outer bounds GL,GU for the block 513 GL = D(IBEGIN) 514 GU = D(IBEGIN) 515 DO 15 I = IBEGIN , IEND 516 GL = MIN( GERS( 2*I-1 ), GL ) 517 GU = MAX( GERS( 2*I ), GU ) 518 15 CONTINUE 519 SPDIAM = GU - GL 520 521 IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN 522* Count the number of eigenvalues in the current block. 523 MB = 0 524 DO 20 I = WBEGIN,MM 525 IF( IBLOCK(I).EQ.JBLK ) THEN 526 MB = MB+1 527 ELSE 528 GOTO 21 529 ENDIF 530 20 CONTINUE 531 21 CONTINUE 532 533 IF( MB.EQ.0) THEN 534* No eigenvalue in the current block lies in the desired range 535* E( IEND ) holds the shift for the initial RRR 536 E( IEND ) = ZERO 537 IBEGIN = IEND + 1 538 GO TO 170 539 ELSE 540 541* Decide whether dqds or bisection is more efficient 542 USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) ) 543 WEND = WBEGIN + MB - 1 544* Calculate gaps for the current block 545* In later stages, when representations for individual 546* eigenvalues are different, we use SIGMA = E( IEND ). 547 SIGMA = ZERO 548 DO 30 I = WBEGIN, WEND - 1 549 WGAP( I ) = MAX( ZERO, 550 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 551 30 CONTINUE 552 WGAP( WEND ) = MAX( ZERO, 553 $ VU - SIGMA - (W( WEND )+WERR( WEND ))) 554* Find local index of the first and last desired evalue. 555 INDL = INDEXW(WBEGIN) 556 INDU = INDEXW( WEND ) 557 ENDIF 558 ENDIF 559 IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN 560* Case of DQDS 561* Find approximations to the extremal eigenvalues of the block 562 CALL DLARRK( IN, 1, GL, GU, D(IBEGIN), 563 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 564 IF( IINFO.NE.0 ) THEN 565 INFO = -1 566 RETURN 567 ENDIF 568 ISLEFT = MAX(GL, TMP - TMP1 569 $ - HNDRD * EPS* ABS(TMP - TMP1)) 570 571 CALL DLARRK( IN, IN, GL, GU, D(IBEGIN), 572 $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) 573 IF( IINFO.NE.0 ) THEN 574 INFO = -1 575 RETURN 576 ENDIF 577 ISRGHT = MIN(GU, TMP + TMP1 578 $ + HNDRD * EPS * ABS(TMP + TMP1)) 579* Improve the estimate of the spectral diameter 580 SPDIAM = ISRGHT - ISLEFT 581 ELSE 582* Case of bisection 583* Find approximations to the wanted extremal eigenvalues 584 ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN) 585 $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) )) 586 ISRGHT = MIN(GU,W(WEND) + WERR(WEND) 587 $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND))) 588 ENDIF 589 590 591* Decide whether the base representation for the current block 592* L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I 593* should be on the left or the right end of the current block. 594* The strategy is to shift to the end which is "more populated" 595* Furthermore, decide whether to use DQDS for the computation of 596* the eigenvalue approximations at the end of DLARRE or bisection. 597* dqds is chosen if all eigenvalues are desired or the number of 598* eigenvalues to be computed is large compared to the blocksize. 599 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 600* If all the eigenvalues have to be computed, we use dqd 601 USEDQD = .TRUE. 602* INDL is the local index of the first eigenvalue to compute 603 INDL = 1 604 INDU = IN 605* MB = number of eigenvalues to compute 606 MB = IN 607 WEND = WBEGIN + MB - 1 608* Define 1/4 and 3/4 points of the spectrum 609 S1 = ISLEFT + FOURTH * SPDIAM 610 S2 = ISRGHT - FOURTH * SPDIAM 611 ELSE 612* DLARRD has computed IBLOCK and INDEXW for each eigenvalue 613* approximation. 614* choose sigma 615 IF( USEDQD ) THEN 616 S1 = ISLEFT + FOURTH * SPDIAM 617 S2 = ISRGHT - FOURTH * SPDIAM 618 ELSE 619 TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL) 620 S1 = MAX(ISLEFT,VL) + FOURTH * TMP 621 S2 = MIN(ISRGHT,VU) - FOURTH * TMP 622 ENDIF 623 ENDIF 624 625* Compute the negcount at the 1/4 and 3/4 points 626 IF(MB.GT.1) THEN 627 CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN), 628 $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO) 629 ENDIF 630 631 IF(MB.EQ.1) THEN 632 SIGMA = GL 633 SGNDEF = ONE 634 ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN 635 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 636 SIGMA = MAX(ISLEFT,GL) 637 ELSEIF( USEDQD ) THEN 638* use Gerschgorin bound as shift to get pos def matrix 639* for dqds 640 SIGMA = ISLEFT 641 ELSE 642* use approximation of the first desired eigenvalue of the 643* block as shift 644 SIGMA = MAX(ISLEFT,VL) 645 ENDIF 646 SGNDEF = ONE 647 ELSE 648 IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN 649 SIGMA = MIN(ISRGHT,GU) 650 ELSEIF( USEDQD ) THEN 651* use Gerschgorin bound as shift to get neg def matrix 652* for dqds 653 SIGMA = ISRGHT 654 ELSE 655* use approximation of the first desired eigenvalue of the 656* block as shift 657 SIGMA = MIN(ISRGHT,VU) 658 ENDIF 659 SGNDEF = -ONE 660 ENDIF 661 662 663* An initial SIGMA has been chosen that will be used for computing 664* T - SIGMA I = L D L^T 665* Define the increment TAU of the shift in case the initial shift 666* needs to be refined to obtain a factorization with not too much 667* element growth. 668 IF( USEDQD ) THEN 669* The initial SIGMA was to the outer end of the spectrum 670* the matrix is definite and we need not retreat. 671 TAU = SPDIAM*EPS*N + TWO*PIVMIN 672 TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) ) 673 ELSE 674 IF(MB.GT.1) THEN 675 CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) 676 AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN)) 677 IF( SGNDEF.EQ.ONE ) THEN 678 TAU = HALF*MAX(WGAP(WBEGIN),AVGAP) 679 TAU = MAX(TAU,WERR(WBEGIN)) 680 ELSE 681 TAU = HALF*MAX(WGAP(WEND-1),AVGAP) 682 TAU = MAX(TAU,WERR(WEND)) 683 ENDIF 684 ELSE 685 TAU = WERR(WBEGIN) 686 ENDIF 687 ENDIF 688* 689 DO 80 IDUM = 1, MAXTRY 690* Compute L D L^T factorization of tridiagonal matrix T - sigma I. 691* Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of 692* pivots in WORK(2*IN+1:3*IN) 693 DPIVOT = D( IBEGIN ) - SIGMA 694 WORK( 1 ) = DPIVOT 695 DMAX = ABS( WORK(1) ) 696 J = IBEGIN 697 DO 70 I = 1, IN - 1 698 WORK( 2*IN+I ) = ONE / WORK( I ) 699 TMP = E( J )*WORK( 2*IN+I ) 700 WORK( IN+I ) = TMP 701 DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J ) 702 WORK( I+1 ) = DPIVOT 703 DMAX = MAX( DMAX, ABS(DPIVOT) ) 704 J = J + 1 705 70 CONTINUE 706* check for element growth 707 IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN 708 NOREP = .TRUE. 709 ELSE 710 NOREP = .FALSE. 711 ENDIF 712 IF( USEDQD .AND. .NOT.NOREP ) THEN 713* Ensure the definiteness of the representation 714* All entries of D (of L D L^T) must have the same sign 715 DO 71 I = 1, IN 716 TMP = SGNDEF*WORK( I ) 717 IF( TMP.LT.ZERO ) NOREP = .TRUE. 718 71 CONTINUE 719 ENDIF 720 IF(NOREP) THEN 721* Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin 722* shift which makes the matrix definite. So we should end up 723* here really only in the case of IRANGE = VALRNG or INDRNG. 724 IF( IDUM.EQ.MAXTRY-1 ) THEN 725 IF( SGNDEF.EQ.ONE ) THEN 726* The fudged Gerschgorin shift should succeed 727 SIGMA = 728 $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN 729 ELSE 730 SIGMA = 731 $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN 732 END IF 733 ELSE 734 SIGMA = SIGMA - SGNDEF * TAU 735 TAU = TWO * TAU 736 END IF 737 ELSE 738* an initial RRR is found 739 GO TO 83 740 END IF 741 80 CONTINUE 742* if the program reaches this point, no base representation could be 743* found in MAXTRY iterations. 744 INFO = 2 745 RETURN 746 747 83 CONTINUE 748* At this point, we have found an initial base representation 749* T - SIGMA I = L D L^T with not too much element growth. 750* Store the shift. 751 E( IEND ) = SIGMA 752* Store D and L. 753 CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 ) 754 CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 ) 755 756 757 IF(MB.GT.1 ) THEN 758* 759* Perturb each entry of the base representation by a small 760* (but random) relative amount to overcome difficulties with 761* glued matrices. 762* 763 DO 122 I = 1, 4 764 ISEED( I ) = 1 765 122 CONTINUE 766 767 CALL DLARNV(2, ISEED, 2*IN-1, WORK(1)) 768 DO 125 I = 1,IN-1 769 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I)) 770 E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I)) 771 125 CONTINUE 772 D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN)) 773* 774 ENDIF 775* 776* Don't update the Gerschgorin intervals because keeping track 777* of the updates would be too much work in DLARRV. 778* We update W instead and use it to locate the proper Gerschgorin 779* intervals. 780 781* Compute the required eigenvalues of L D L' by bisection or dqds 782 IF ( .NOT.USEDQD ) THEN 783* If DLARRD has been used, shift the eigenvalue approximations 784* according to their representation. This is necessary for 785* a uniform DLARRV since dqds computes eigenvalues of the 786* shifted representation. In DLARRV, W will always hold the 787* UNshifted eigenvalue approximation. 788 DO 134 J=WBEGIN,WEND 789 W(J) = W(J) - SIGMA 790 WERR(J) = WERR(J) + ABS(W(J)) * EPS 791 134 CONTINUE 792* call DLARRB to reduce eigenvalue error of the approximations 793* from DLARRD 794 DO 135 I = IBEGIN, IEND-1 795 WORK( I ) = D( I ) * E( I )**2 796 135 CONTINUE 797* use bisection to find EV from INDL to INDU 798 CALL DLARRB(IN, D(IBEGIN), WORK(IBEGIN), 799 $ INDL, INDU, RTOL1, RTOL2, INDL-1, 800 $ W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN), 801 $ WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM, 802 $ IN, IINFO ) 803 IF( IINFO .NE. 0 ) THEN 804 INFO = -4 805 RETURN 806 END IF 807* DLARRB computes all gaps correctly except for the last one 808* Record distance to VU/GU 809 WGAP( WEND ) = MAX( ZERO, 810 $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) ) 811 DO 138 I = INDL, INDU 812 M = M + 1 813 IBLOCK(M) = JBLK 814 INDEXW(M) = I 815 138 CONTINUE 816 ELSE 817* Call dqds to get all eigs (and then possibly delete unwanted 818* eigenvalues). 819* Note that dqds finds the eigenvalues of the L D L^T representation 820* of T to high relative accuracy. High relative accuracy 821* might be lost when the shift of the RRR is subtracted to obtain 822* the eigenvalues of T. However, T is not guaranteed to define its 823* eigenvalues to high relative accuracy anyway. 824* Set RTOL to the order of the tolerance used in DLASQ2 825* This is an ESTIMATED error, the worst case bound is 4*N*EPS 826* which is usually too large and requires unnecessary work to be 827* done by bisection when computing the eigenvectors 828 RTOL = LOG(DBLE(IN)) * FOUR * EPS 829 J = IBEGIN 830 DO 140 I = 1, IN - 1 831 WORK( 2*I-1 ) = ABS( D( J ) ) 832 WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 ) 833 J = J + 1 834 140 CONTINUE 835 WORK( 2*IN-1 ) = ABS( D( IEND ) ) 836 WORK( 2*IN ) = ZERO 837 CALL DLASQ2( IN, WORK, IINFO ) 838 IF( IINFO .NE. 0 ) THEN 839* If IINFO = -5 then an index is part of a tight cluster 840* and should be changed. The index is in IWORK(1) and the 841* gap is in WORK(N+1) 842 INFO = -5 843 RETURN 844 ELSE 845* Test that all eigenvalues are positive as expected 846 DO 149 I = 1, IN 847 IF( WORK( I ).LT.ZERO ) THEN 848 INFO = -6 849 RETURN 850 ENDIF 851 149 CONTINUE 852 END IF 853 IF( SGNDEF.GT.ZERO ) THEN 854 DO 150 I = INDL, INDU 855 M = M + 1 856 W( M ) = WORK( IN-I+1 ) 857 IBLOCK( M ) = JBLK 858 INDEXW( M ) = I 859 150 CONTINUE 860 ELSE 861 DO 160 I = INDL, INDU 862 M = M + 1 863 W( M ) = -WORK( I ) 864 IBLOCK( M ) = JBLK 865 INDEXW( M ) = I 866 160 CONTINUE 867 END IF 868 869 DO 165 I = M - MB + 1, M 870* the value of RTOL below should be the tolerance in DLASQ2 871 WERR( I ) = RTOL * ABS( W(I) ) 872 165 CONTINUE 873 DO 166 I = M - MB + 1, M - 1 874* compute the right gap between the intervals 875 WGAP( I ) = MAX( ZERO, 876 $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 877 166 CONTINUE 878 WGAP( M ) = MAX( ZERO, 879 $ ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) ) 880 END IF 881* proceed with next block 882 IBEGIN = IEND + 1 883 WBEGIN = WEND + 1 884 170 CONTINUE 885* 886 887 RETURN 888* 889* end of DLARRE 890* 891 END 892