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