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