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