1*> \brief \b SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues of L D LT. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLARRV + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrv.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrv.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrv.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN, 22* ISPLIT, M, DOL, DOU, MINRGP, 23* RTOL1, RTOL2, W, WERR, WGAP, 24* IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 25* WORK, IWORK, INFO ) 26* 27* .. Scalar Arguments .. 28* INTEGER DOL, DOU, INFO, LDZ, M, N 29* REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 30* .. 31* .. Array Arguments .. 32* INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 33* $ ISUPPZ( * ), IWORK( * ) 34* REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 35* $ WGAP( * ), WORK( * ) 36* REAL Z( LDZ, * ) 37* .. 38* 39* 40*> \par Purpose: 41* ============= 42*> 43*> \verbatim 44*> 45*> SLARRV computes the eigenvectors of the tridiagonal matrix 46*> T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. 47*> The input eigenvalues should have been computed by SLARRE. 48*> \endverbatim 49* 50* Arguments: 51* ========== 52* 53*> \param[in] N 54*> \verbatim 55*> N is INTEGER 56*> The order of the matrix. N >= 0. 57*> \endverbatim 58*> 59*> \param[in] VL 60*> \verbatim 61*> VL is REAL 62*> Lower bound of the interval that contains the desired 63*> eigenvalues. VL < VU. Needed to compute gaps on the left or right 64*> end of the extremal eigenvalues in the desired RANGE. 65*> \endverbatim 66*> 67*> \param[in] VU 68*> \verbatim 69*> VU is REAL 70*> Upper bound of the interval that contains the desired 71*> eigenvalues. VL < VU. 72*> Note: VU is currently not used by this implementation of SLARRV, VU is 73*> passed to SLARRV because it could be used compute gaps on the right end 74*> of the extremal eigenvalues. However, with not much initial accuracy in 75*> LAMBDA and VU, the formula can lead to an overestimation of the right gap 76*> and thus to inadequately early RQI 'convergence'. This is currently 77*> prevented this by forcing a small right gap. And so it turns out that VU 78*> is currently not used by this implementation of SLARRV. 79*> \endverbatim 80*> 81*> \param[in,out] D 82*> \verbatim 83*> D is REAL array, dimension (N) 84*> On entry, the N diagonal elements of the diagonal matrix D. 85*> On exit, D may be overwritten. 86*> \endverbatim 87*> 88*> \param[in,out] L 89*> \verbatim 90*> L is REAL array, dimension (N) 91*> On entry, the (N-1) subdiagonal elements of the unit 92*> bidiagonal matrix L are in elements 1 to N-1 of L 93*> (if the matrix is not split.) At the end of each block 94*> is stored the corresponding shift as given by SLARRE. 95*> On exit, L is overwritten. 96*> \endverbatim 97*> 98*> \param[in] PIVMIN 99*> \verbatim 100*> PIVMIN is REAL 101*> The minimum pivot allowed in the Sturm sequence. 102*> \endverbatim 103*> 104*> \param[in] ISPLIT 105*> \verbatim 106*> ISPLIT is 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 109*> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 110*> through ISPLIT( 2 ), etc. 111*> \endverbatim 112*> 113*> \param[in] M 114*> \verbatim 115*> M is INTEGER 116*> The total number of input eigenvalues. 0 <= M <= N. 117*> \endverbatim 118*> 119*> \param[in] DOL 120*> \verbatim 121*> DOL is INTEGER 122*> \endverbatim 123*> 124*> \param[in] DOU 125*> \verbatim 126*> DOU is INTEGER 127*> If the user wants to compute only selected eigenvectors from all 128*> the eigenvalues supplied, he can specify an index range DOL:DOU. 129*> Or else the setting DOL=1, DOU=M should be applied. 130*> Note that DOL and DOU refer to the order in which the eigenvalues 131*> are stored in W. 132*> If the user wants to compute only selected eigenpairs, then 133*> the columns DOL-1 to DOU+1 of the eigenvector space Z contain the 134*> computed eigenvectors. All other columns of Z are set to zero. 135*> \endverbatim 136*> 137*> \param[in] MINRGP 138*> \verbatim 139*> MINRGP is REAL 140*> \endverbatim 141*> 142*> \param[in] RTOL1 143*> \verbatim 144*> RTOL1 is REAL 145*> \endverbatim 146*> 147*> \param[in] RTOL2 148*> \verbatim 149*> RTOL2 is REAL 150*> Parameters for bisection. 151*> An interval [LEFT,RIGHT] has converged if 152*> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) 153*> \endverbatim 154*> 155*> \param[in,out] W 156*> \verbatim 157*> W is REAL array, dimension (N) 158*> The first M elements of W contain the APPROXIMATE eigenvalues for 159*> which eigenvectors are to be computed. The eigenvalues 160*> should be grouped by split-off block and ordered from 161*> smallest to largest within the block ( The output array 162*> W from SLARRE is expected here ). Furthermore, they are with 163*> respect to the shift of the corresponding root representation 164*> for their block. On exit, W holds the eigenvalues of the 165*> UNshifted matrix. 166*> \endverbatim 167*> 168*> \param[in,out] WERR 169*> \verbatim 170*> WERR is REAL array, dimension (N) 171*> The first M elements contain the semiwidth of the uncertainty 172*> interval of the corresponding eigenvalue in W 173*> \endverbatim 174*> 175*> \param[in,out] WGAP 176*> \verbatim 177*> WGAP is REAL array, dimension (N) 178*> The separation from the right neighbor eigenvalue in W. 179*> \endverbatim 180*> 181*> \param[in] IBLOCK 182*> \verbatim 183*> IBLOCK is INTEGER array, dimension (N) 184*> The indices of the blocks (submatrices) associated with the 185*> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue 186*> W(i) belongs to the first block from the top, =2 if W(i) 187*> belongs to the second block, etc. 188*> \endverbatim 189*> 190*> \param[in] INDEXW 191*> \verbatim 192*> INDEXW is INTEGER array, dimension (N) 193*> The indices of the eigenvalues within each block (submatrix); 194*> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the 195*> i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. 196*> \endverbatim 197*> 198*> \param[in] GERS 199*> \verbatim 200*> GERS is REAL array, dimension (2*N) 201*> The N Gerschgorin intervals (the i-th Gerschgorin interval 202*> is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should 203*> be computed from the original UNshifted matrix. 204*> \endverbatim 205*> 206*> \param[out] Z 207*> \verbatim 208*> Z is REAL array, dimension (LDZ, max(1,M) ) 209*> If INFO = 0, the first M columns of Z contain the 210*> orthonormal eigenvectors of the matrix T 211*> corresponding to the input eigenvalues, with the i-th 212*> column of Z holding the eigenvector associated with W(i). 213*> Note: the user must ensure that at least max(1,M) columns are 214*> supplied in the array Z. 215*> \endverbatim 216*> 217*> \param[in] LDZ 218*> \verbatim 219*> LDZ is INTEGER 220*> The leading dimension of the array Z. LDZ >= 1, and if 221*> JOBZ = 'V', LDZ >= max(1,N). 222*> \endverbatim 223*> 224*> \param[out] ISUPPZ 225*> \verbatim 226*> ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) 227*> The support of the eigenvectors in Z, i.e., the indices 228*> indicating the nonzero elements in Z. The I-th eigenvector 229*> is nonzero only in elements ISUPPZ( 2*I-1 ) through 230*> ISUPPZ( 2*I ). 231*> \endverbatim 232*> 233*> \param[out] WORK 234*> \verbatim 235*> WORK is REAL array, dimension (12*N) 236*> \endverbatim 237*> 238*> \param[out] IWORK 239*> \verbatim 240*> IWORK is INTEGER array, dimension (7*N) 241*> \endverbatim 242*> 243*> \param[out] INFO 244*> \verbatim 245*> INFO is INTEGER 246*> = 0: successful exit 247*> 248*> > 0: A problem occurred in SLARRV. 249*> < 0: One of the called subroutines signaled an internal problem. 250*> Needs inspection of the corresponding parameter IINFO 251*> for further information. 252*> 253*> =-1: Problem in SLARRB when refining a child's eigenvalues. 254*> =-2: Problem in SLARRF when computing the RRR of a child. 255*> When a child is inside a tight cluster, it can be difficult 256*> to find an RRR. A partial remedy from the user's point of 257*> view is to make the parameter MINRGP smaller and recompile. 258*> However, as the orthogonality of the computed vectors is 259*> proportional to 1/MINRGP, the user should be aware that 260*> he might be trading in precision when he decreases MINRGP. 261*> =-3: Problem in SLARRB when refining a single eigenvalue 262*> after the Rayleigh correction was rejected. 263*> = 5: The Rayleigh Quotient Iteration failed to converge to 264*> full accuracy in MAXITR steps. 265*> \endverbatim 266* 267* Authors: 268* ======== 269* 270*> \author Univ. of Tennessee 271*> \author Univ. of California Berkeley 272*> \author Univ. of Colorado Denver 273*> \author NAG Ltd. 274* 275*> \ingroup realOTHERauxiliary 276* 277*> \par Contributors: 278* ================== 279*> 280*> Beresford Parlett, University of California, Berkeley, USA \n 281*> Jim Demmel, University of California, Berkeley, USA \n 282*> Inderjit Dhillon, University of Texas, Austin, USA \n 283*> Osni Marques, LBNL/NERSC, USA \n 284*> Christof Voemel, University of California, Berkeley, USA 285* 286* ===================================================================== 287 SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN, 288 $ ISPLIT, M, DOL, DOU, MINRGP, 289 $ RTOL1, RTOL2, W, WERR, WGAP, 290 $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 291 $ WORK, IWORK, INFO ) 292* 293* -- LAPACK auxiliary routine -- 294* -- LAPACK is a software package provided by Univ. of Tennessee, -- 295* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 296* 297* .. Scalar Arguments .. 298 INTEGER DOL, DOU, INFO, LDZ, M, N 299 REAL MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 300* .. 301* .. Array Arguments .. 302 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 303 $ ISUPPZ( * ), IWORK( * ) 304 REAL D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 305 $ WGAP( * ), WORK( * ) 306 REAL Z( LDZ, * ) 307* .. 308* 309* ===================================================================== 310* 311* .. Parameters .. 312 INTEGER MAXITR 313 PARAMETER ( MAXITR = 10 ) 314 REAL ZERO, ONE, TWO, THREE, FOUR, HALF 315 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, 316 $ TWO = 2.0E0, THREE = 3.0E0, 317 $ FOUR = 4.0E0, HALF = 0.5E0) 318* .. 319* .. Local Scalars .. 320 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ 321 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1, 322 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG, 323 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER, 324 $ ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS, 325 $ NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST, 326 $ NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST, 327 $ OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX, 328 $ WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU, 329 $ ZUSEDW 330 REAL BSTRES, BSTW, EPS, FUDGE, GAP, GAPTOL, GL, GU, 331 $ LAMBDA, LEFT, LGAP, MINGMA, NRMINV, RESID, 332 $ RGAP, RIGHT, RQCORR, RQTOL, SAVGAP, SGNDEF, 333 $ SIGMA, SPDIAM, SSIGMA, TAU, TMP, TOL, ZTZ 334* .. 335* .. External Functions .. 336 REAL SLAMCH 337 EXTERNAL SLAMCH 338* .. 339* .. External Subroutines .. 340 EXTERNAL SCOPY, SLAR1V, SLARRB, SLARRF, SLASET, 341 $ SSCAL 342* .. 343* .. Intrinsic Functions .. 344 INTRINSIC ABS, REAL, MAX, MIN 345* .. 346* .. Executable Statements .. 347* .. 348 349 INFO = 0 350* 351* Quick return if possible 352* 353 IF( N.LE.0 ) THEN 354 RETURN 355 END IF 356* 357* The first N entries of WORK are reserved for the eigenvalues 358 INDLD = N+1 359 INDLLD= 2*N+1 360 INDWRK= 3*N+1 361 MINWSIZE = 12 * N 362 363 DO 5 I= 1,MINWSIZE 364 WORK( I ) = ZERO 365 5 CONTINUE 366 367* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the 368* factorization used to compute the FP vector 369 IINDR = 0 370* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current 371* layer and the one above. 372 IINDC1 = N 373 IINDC2 = 2*N 374 IINDWK = 3*N + 1 375 376 MINIWSIZE = 7 * N 377 DO 10 I= 1,MINIWSIZE 378 IWORK( I ) = 0 379 10 CONTINUE 380 381 ZUSEDL = 1 382 IF(DOL.GT.1) THEN 383* Set lower bound for use of Z 384 ZUSEDL = DOL-1 385 ENDIF 386 ZUSEDU = M 387 IF(DOU.LT.M) THEN 388* Set lower bound for use of Z 389 ZUSEDU = DOU+1 390 ENDIF 391* The width of the part of Z that is used 392 ZUSEDW = ZUSEDU - ZUSEDL + 1 393 394 395 CALL SLASET( 'Full', N, ZUSEDW, ZERO, ZERO, 396 $ Z(1,ZUSEDL), LDZ ) 397 398 EPS = SLAMCH( 'Precision' ) 399 RQTOL = TWO * EPS 400* 401* Set expert flags for standard code. 402 TRYRQC = .TRUE. 403 404 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 405 ELSE 406* Only selected eigenpairs are computed. Since the other evalues 407* are not refined by RQ iteration, bisection has to compute to full 408* accuracy. 409 RTOL1 = FOUR * EPS 410 RTOL2 = FOUR * EPS 411 ENDIF 412 413* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the 414* desired eigenvalues. The support of the nonzero eigenvector 415* entries is contained in the interval IBEGIN:IEND. 416* Remark that if k eigenpairs are desired, then the eigenvectors 417* are stored in k contiguous columns of Z. 418 419* DONE is the number of eigenvectors already computed 420 DONE = 0 421 IBEGIN = 1 422 WBEGIN = 1 423 DO 170 JBLK = 1, IBLOCK( M ) 424 IEND = ISPLIT( JBLK ) 425 SIGMA = L( IEND ) 426* Find the eigenvectors of the submatrix indexed IBEGIN 427* through IEND. 428 WEND = WBEGIN - 1 429 15 CONTINUE 430 IF( WEND.LT.M ) THEN 431 IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN 432 WEND = WEND + 1 433 GO TO 15 434 END IF 435 END IF 436 IF( WEND.LT.WBEGIN ) THEN 437 IBEGIN = IEND + 1 438 GO TO 170 439 ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN 440 IBEGIN = IEND + 1 441 WBEGIN = WEND + 1 442 GO TO 170 443 END IF 444 445* Find local spectral diameter of the block 446 GL = GERS( 2*IBEGIN-1 ) 447 GU = GERS( 2*IBEGIN ) 448 DO 20 I = IBEGIN+1 , IEND 449 GL = MIN( GERS( 2*I-1 ), GL ) 450 GU = MAX( GERS( 2*I ), GU ) 451 20 CONTINUE 452 SPDIAM = GU - GL 453 454* OLDIEN is the last index of the previous block 455 OLDIEN = IBEGIN - 1 456* Calculate the size of the current block 457 IN = IEND - IBEGIN + 1 458* The number of eigenvalues in the current block 459 IM = WEND - WBEGIN + 1 460 461* This is for a 1x1 block 462 IF( IBEGIN.EQ.IEND ) THEN 463 DONE = DONE+1 464 Z( IBEGIN, WBEGIN ) = ONE 465 ISUPPZ( 2*WBEGIN-1 ) = IBEGIN 466 ISUPPZ( 2*WBEGIN ) = IBEGIN 467 W( WBEGIN ) = W( WBEGIN ) + SIGMA 468 WORK( WBEGIN ) = W( WBEGIN ) 469 IBEGIN = IEND + 1 470 WBEGIN = WBEGIN + 1 471 GO TO 170 472 END IF 473 474* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) 475* Note that these can be approximations, in this case, the corresp. 476* entries of WERR give the size of the uncertainty interval. 477* The eigenvalue approximations will be refined when necessary as 478* high relative accuracy is required for the computation of the 479* corresponding eigenvectors. 480 CALL SCOPY( IM, W( WBEGIN ), 1, 481 $ WORK( WBEGIN ), 1 ) 482 483* We store in W the eigenvalue approximations w.r.t. the original 484* matrix T. 485 DO 30 I=1,IM 486 W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA 487 30 CONTINUE 488 489 490* NDEPTH is the current depth of the representation tree 491 NDEPTH = 0 492* PARITY is either 1 or 0 493 PARITY = 1 494* NCLUS is the number of clusters for the next level of the 495* representation tree, we start with NCLUS = 1 for the root 496 NCLUS = 1 497 IWORK( IINDC1+1 ) = 1 498 IWORK( IINDC1+2 ) = IM 499 500* IDONE is the number of eigenvectors already computed in the current 501* block 502 IDONE = 0 503* loop while( IDONE.LT.IM ) 504* generate the representation tree for the current block and 505* compute the eigenvectors 506 40 CONTINUE 507 IF( IDONE.LT.IM ) THEN 508* This is a crude protection against infinitely deep trees 509 IF( NDEPTH.GT.M ) THEN 510 INFO = -2 511 RETURN 512 ENDIF 513* breadth first processing of the current level of the representation 514* tree: OLDNCL = number of clusters on current level 515 OLDNCL = NCLUS 516* reset NCLUS to count the number of child clusters 517 NCLUS = 0 518* 519 PARITY = 1 - PARITY 520 IF( PARITY.EQ.0 ) THEN 521 OLDCLS = IINDC1 522 NEWCLS = IINDC2 523 ELSE 524 OLDCLS = IINDC2 525 NEWCLS = IINDC1 526 END IF 527* Process the clusters on the current level 528 DO 150 I = 1, OLDNCL 529 J = OLDCLS + 2*I 530* OLDFST, OLDLST = first, last index of current cluster. 531* cluster indices start with 1 and are relative 532* to WBEGIN when accessing W, WGAP, WERR, Z 533 OLDFST = IWORK( J-1 ) 534 OLDLST = IWORK( J ) 535 IF( NDEPTH.GT.0 ) THEN 536* Retrieve relatively robust representation (RRR) of cluster 537* that has been computed at the previous level 538* The RRR is stored in Z and overwritten once the eigenvectors 539* have been computed or when the cluster is refined 540 541 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 542* Get representation from location of the leftmost evalue 543* of the cluster 544 J = WBEGIN + OLDFST - 1 545 ELSE 546 IF(WBEGIN+OLDFST-1.LT.DOL) THEN 547* Get representation from the left end of Z array 548 J = DOL - 1 549 ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN 550* Get representation from the right end of Z array 551 J = DOU 552 ELSE 553 J = WBEGIN + OLDFST - 1 554 ENDIF 555 ENDIF 556 CALL SCOPY( IN, Z( IBEGIN, J ), 1, D( IBEGIN ), 1 ) 557 CALL SCOPY( IN-1, Z( IBEGIN, J+1 ), 1, L( IBEGIN ), 558 $ 1 ) 559 SIGMA = Z( IEND, J+1 ) 560 561* Set the corresponding entries in Z to zero 562 CALL SLASET( 'Full', IN, 2, ZERO, ZERO, 563 $ Z( IBEGIN, J), LDZ ) 564 END IF 565 566* Compute DL and DLL of current RRR 567 DO 50 J = IBEGIN, IEND-1 568 TMP = D( J )*L( J ) 569 WORK( INDLD-1+J ) = TMP 570 WORK( INDLLD-1+J ) = TMP*L( J ) 571 50 CONTINUE 572 573 IF( NDEPTH.GT.0 ) THEN 574* P and Q are index of the first and last eigenvalue to compute 575* within the current block 576 P = INDEXW( WBEGIN-1+OLDFST ) 577 Q = INDEXW( WBEGIN-1+OLDLST ) 578* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 579* through the Q-OFFSET elements of these arrays are to be used. 580* OFFSET = P-OLDFST 581 OFFSET = INDEXW( WBEGIN ) - 1 582* perform limited bisection (if necessary) to get approximate 583* eigenvalues to the precision needed. 584 CALL SLARRB( IN, D( IBEGIN ), 585 $ WORK(INDLLD+IBEGIN-1), 586 $ P, Q, RTOL1, RTOL2, OFFSET, 587 $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN), 588 $ WORK( INDWRK ), IWORK( IINDWK ), 589 $ PIVMIN, SPDIAM, IN, IINFO ) 590 IF( IINFO.NE.0 ) THEN 591 INFO = -1 592 RETURN 593 ENDIF 594* We also recompute the extremal gaps. W holds all eigenvalues 595* of the unshifted matrix and must be used for computation 596* of WGAP, the entries of WORK might stem from RRRs with 597* different shifts. The gaps from WBEGIN-1+OLDFST to 598* WBEGIN-1+OLDLST are correctly computed in SLARRB. 599* However, we only allow the gaps to become greater since 600* this is what should happen when we decrease WERR 601 IF( OLDFST.GT.1) THEN 602 WGAP( WBEGIN+OLDFST-2 ) = 603 $ MAX(WGAP(WBEGIN+OLDFST-2), 604 $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1) 605 $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) ) 606 ENDIF 607 IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN 608 WGAP( WBEGIN+OLDLST-1 ) = 609 $ MAX(WGAP(WBEGIN+OLDLST-1), 610 $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST) 611 $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) ) 612 ENDIF 613* Each time the eigenvalues in WORK get refined, we store 614* the newly found approximation with all shifts applied in W 615 DO 53 J=OLDFST,OLDLST 616 W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA 617 53 CONTINUE 618 END IF 619 620* Process the current node. 621 NEWFST = OLDFST 622 DO 140 J = OLDFST, OLDLST 623 IF( J.EQ.OLDLST ) THEN 624* we are at the right end of the cluster, this is also the 625* boundary of the child cluster 626 NEWLST = J 627 ELSE IF ( WGAP( WBEGIN + J -1).GE. 628 $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN 629* the right relative gap is big enough, the child cluster 630* (NEWFST,..,NEWLST) is well separated from the following 631 NEWLST = J 632 ELSE 633* inside a child cluster, the relative gap is not 634* big enough. 635 GOTO 140 636 END IF 637 638* Compute size of child cluster found 639 NEWSIZ = NEWLST - NEWFST + 1 640 641* NEWFTT is the place in Z where the new RRR or the computed 642* eigenvector is to be stored 643 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 644* Store representation at location of the leftmost evalue 645* of the cluster 646 NEWFTT = WBEGIN + NEWFST - 1 647 ELSE 648 IF(WBEGIN+NEWFST-1.LT.DOL) THEN 649* Store representation at the left end of Z array 650 NEWFTT = DOL - 1 651 ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN 652* Store representation at the right end of Z array 653 NEWFTT = DOU 654 ELSE 655 NEWFTT = WBEGIN + NEWFST - 1 656 ENDIF 657 ENDIF 658 659 IF( NEWSIZ.GT.1) THEN 660* 661* Current child is not a singleton but a cluster. 662* Compute and store new representation of child. 663* 664* 665* Compute left and right cluster gap. 666* 667* LGAP and RGAP are not computed from WORK because 668* the eigenvalue approximations may stem from RRRs 669* different shifts. However, W hold all eigenvalues 670* of the unshifted matrix. Still, the entries in WGAP 671* have to be computed from WORK since the entries 672* in W might be of the same order so that gaps are not 673* exhibited correctly for very close eigenvalues. 674 IF( NEWFST.EQ.1 ) THEN 675 LGAP = MAX( ZERO, 676 $ W(WBEGIN)-WERR(WBEGIN) - VL ) 677 ELSE 678 LGAP = WGAP( WBEGIN+NEWFST-2 ) 679 ENDIF 680 RGAP = WGAP( WBEGIN+NEWLST-1 ) 681* 682* Compute left- and rightmost eigenvalue of child 683* to high precision in order to shift as close 684* as possible and obtain as large relative gaps 685* as possible 686* 687 DO 55 K =1,2 688 IF(K.EQ.1) THEN 689 P = INDEXW( WBEGIN-1+NEWFST ) 690 ELSE 691 P = INDEXW( WBEGIN-1+NEWLST ) 692 ENDIF 693 OFFSET = INDEXW( WBEGIN ) - 1 694 CALL SLARRB( IN, D(IBEGIN), 695 $ WORK( INDLLD+IBEGIN-1 ),P,P, 696 $ RQTOL, RQTOL, OFFSET, 697 $ WORK(WBEGIN),WGAP(WBEGIN), 698 $ WERR(WBEGIN),WORK( INDWRK ), 699 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 700 $ IN, IINFO ) 701 55 CONTINUE 702* 703 IF((WBEGIN+NEWLST-1.LT.DOL).OR. 704 $ (WBEGIN+NEWFST-1.GT.DOU)) THEN 705* if the cluster contains no desired eigenvalues 706* skip the computation of that branch of the rep. tree 707* 708* We could skip before the refinement of the extremal 709* eigenvalues of the child, but then the representation 710* tree could be different from the one when nothing is 711* skipped. For this reason we skip at this place. 712 IDONE = IDONE + NEWLST - NEWFST + 1 713 GOTO 139 714 ENDIF 715* 716* Compute RRR of child cluster. 717* Note that the new RRR is stored in Z 718* 719* SLARRF needs LWORK = 2*N 720 CALL SLARRF( IN, D( IBEGIN ), L( IBEGIN ), 721 $ WORK(INDLD+IBEGIN-1), 722 $ NEWFST, NEWLST, WORK(WBEGIN), 723 $ WGAP(WBEGIN), WERR(WBEGIN), 724 $ SPDIAM, LGAP, RGAP, PIVMIN, TAU, 725 $ Z(IBEGIN, NEWFTT),Z(IBEGIN, NEWFTT+1), 726 $ WORK( INDWRK ), IINFO ) 727 IF( IINFO.EQ.0 ) THEN 728* a new RRR for the cluster was found by SLARRF 729* update shift and store it 730 SSIGMA = SIGMA + TAU 731 Z( IEND, NEWFTT+1 ) = SSIGMA 732* WORK() are the midpoints and WERR() the semi-width 733* Note that the entries in W are unchanged. 734 DO 116 K = NEWFST, NEWLST 735 FUDGE = 736 $ THREE*EPS*ABS(WORK(WBEGIN+K-1)) 737 WORK( WBEGIN + K - 1 ) = 738 $ WORK( WBEGIN + K - 1) - TAU 739 FUDGE = FUDGE + 740 $ FOUR*EPS*ABS(WORK(WBEGIN+K-1)) 741* Fudge errors 742 WERR( WBEGIN + K - 1 ) = 743 $ WERR( WBEGIN + K - 1 ) + FUDGE 744* Gaps are not fudged. Provided that WERR is small 745* when eigenvalues are close, a zero gap indicates 746* that a new representation is needed for resolving 747* the cluster. A fudge could lead to a wrong decision 748* of judging eigenvalues 'separated' which in 749* reality are not. This could have a negative impact 750* on the orthogonality of the computed eigenvectors. 751 116 CONTINUE 752 753 NCLUS = NCLUS + 1 754 K = NEWCLS + 2*NCLUS 755 IWORK( K-1 ) = NEWFST 756 IWORK( K ) = NEWLST 757 ELSE 758 INFO = -2 759 RETURN 760 ENDIF 761 ELSE 762* 763* Compute eigenvector of singleton 764* 765 ITER = 0 766* 767 TOL = FOUR * LOG(REAL(IN)) * EPS 768* 769 K = NEWFST 770 WINDEX = WBEGIN + K - 1 771 WINDMN = MAX(WINDEX - 1,1) 772 WINDPL = MIN(WINDEX + 1,M) 773 LAMBDA = WORK( WINDEX ) 774 DONE = DONE + 1 775* Check if eigenvector computation is to be skipped 776 IF((WINDEX.LT.DOL).OR. 777 $ (WINDEX.GT.DOU)) THEN 778 ESKIP = .TRUE. 779 GOTO 125 780 ELSE 781 ESKIP = .FALSE. 782 ENDIF 783 LEFT = WORK( WINDEX ) - WERR( WINDEX ) 784 RIGHT = WORK( WINDEX ) + WERR( WINDEX ) 785 INDEIG = INDEXW( WINDEX ) 786* Note that since we compute the eigenpairs for a child, 787* all eigenvalue approximations are w.r.t the same shift. 788* In this case, the entries in WORK should be used for 789* computing the gaps since they exhibit even very small 790* differences in the eigenvalues, as opposed to the 791* entries in W which might "look" the same. 792 793 IF( K .EQ. 1) THEN 794* In the case RANGE='I' and with not much initial 795* accuracy in LAMBDA and VL, the formula 796* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 797* can lead to an overestimation of the left gap and 798* thus to inadequately early RQI 'convergence'. 799* Prevent this by forcing a small left gap. 800 LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 801 ELSE 802 LGAP = WGAP(WINDMN) 803 ENDIF 804 IF( K .EQ. IM) THEN 805* In the case RANGE='I' and with not much initial 806* accuracy in LAMBDA and VU, the formula 807* can lead to an overestimation of the right gap and 808* thus to inadequately early RQI 'convergence'. 809* Prevent this by forcing a small right gap. 810 RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 811 ELSE 812 RGAP = WGAP(WINDEX) 813 ENDIF 814 GAP = MIN( LGAP, RGAP ) 815 IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN 816* The eigenvector support can become wrong 817* because significant entries could be cut off due to a 818* large GAPTOL parameter in LAR1V. Prevent this. 819 GAPTOL = ZERO 820 ELSE 821 GAPTOL = GAP * EPS 822 ENDIF 823 ISUPMN = IN 824 ISUPMX = 1 825* Update WGAP so that it holds the minimum gap 826* to the left or the right. This is crucial in the 827* case where bisection is used to ensure that the 828* eigenvalue is refined up to the required precision. 829* The correct value is restored afterwards. 830 SAVGAP = WGAP(WINDEX) 831 WGAP(WINDEX) = GAP 832* We want to use the Rayleigh Quotient Correction 833* as often as possible since it converges quadratically 834* when we are close enough to the desired eigenvalue. 835* However, the Rayleigh Quotient can have the wrong sign 836* and lead us away from the desired eigenvalue. In this 837* case, the best we can do is to use bisection. 838 USEDBS = .FALSE. 839 USEDRQ = .FALSE. 840* Bisection is initially turned off unless it is forced 841 NEEDBS = .NOT.TRYRQC 842 120 CONTINUE 843* Check if bisection should be used to refine eigenvalue 844 IF(NEEDBS) THEN 845* Take the bisection as new iterate 846 USEDBS = .TRUE. 847 ITMP1 = IWORK( IINDR+WINDEX ) 848 OFFSET = INDEXW( WBEGIN ) - 1 849 CALL SLARRB( IN, D(IBEGIN), 850 $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG, 851 $ ZERO, TWO*EPS, OFFSET, 852 $ WORK(WBEGIN),WGAP(WBEGIN), 853 $ WERR(WBEGIN),WORK( INDWRK ), 854 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 855 $ ITMP1, IINFO ) 856 IF( IINFO.NE.0 ) THEN 857 INFO = -3 858 RETURN 859 ENDIF 860 LAMBDA = WORK( WINDEX ) 861* Reset twist index from inaccurate LAMBDA to 862* force computation of true MINGMA 863 IWORK( IINDR+WINDEX ) = 0 864 ENDIF 865* Given LAMBDA, compute the eigenvector. 866 CALL SLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 867 $ L( IBEGIN ), WORK(INDLD+IBEGIN-1), 868 $ WORK(INDLLD+IBEGIN-1), 869 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 870 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 871 $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ), 872 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 873 IF(ITER .EQ. 0) THEN 874 BSTRES = RESID 875 BSTW = LAMBDA 876 ELSEIF(RESID.LT.BSTRES) THEN 877 BSTRES = RESID 878 BSTW = LAMBDA 879 ENDIF 880 ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 )) 881 ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX )) 882 ITER = ITER + 1 883 884* sin alpha <= |resid|/gap 885* Note that both the residual and the gap are 886* proportional to the matrix, so ||T|| doesn't play 887* a role in the quotient 888 889* 890* Convergence test for Rayleigh-Quotient iteration 891* (omitted when Bisection has been used) 892* 893 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 894 $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS) 895 $ THEN 896* We need to check that the RQCORR update doesn't 897* move the eigenvalue away from the desired one and 898* towards a neighbor. -> protection with bisection 899 IF(INDEIG.LE.NEGCNT) THEN 900* The wanted eigenvalue lies to the left 901 SGNDEF = -ONE 902 ELSE 903* The wanted eigenvalue lies to the right 904 SGNDEF = ONE 905 ENDIF 906* We only use the RQCORR if it improves the 907* the iterate reasonably. 908 IF( ( RQCORR*SGNDEF.GE.ZERO ) 909 $ .AND.( LAMBDA + RQCORR.LE. RIGHT) 910 $ .AND.( LAMBDA + RQCORR.GE. LEFT) 911 $ ) THEN 912 USEDRQ = .TRUE. 913* Store new midpoint of bisection interval in WORK 914 IF(SGNDEF.EQ.ONE) THEN 915* The current LAMBDA is on the left of the true 916* eigenvalue 917 LEFT = LAMBDA 918* We prefer to assume that the error estimate 919* is correct. We could make the interval not 920* as a bracket but to be modified if the RQCORR 921* chooses to. In this case, the RIGHT side should 922* be modified as follows: 923* RIGHT = MAX(RIGHT, LAMBDA + RQCORR) 924 ELSE 925* The current LAMBDA is on the right of the true 926* eigenvalue 927 RIGHT = LAMBDA 928* See comment about assuming the error estimate is 929* correct above. 930* LEFT = MIN(LEFT, LAMBDA + RQCORR) 931 ENDIF 932 WORK( WINDEX ) = 933 $ HALF * (RIGHT + LEFT) 934* Take RQCORR since it has the correct sign and 935* improves the iterate reasonably 936 LAMBDA = LAMBDA + RQCORR 937* Update width of error interval 938 WERR( WINDEX ) = 939 $ HALF * (RIGHT-LEFT) 940 ELSE 941 NEEDBS = .TRUE. 942 ENDIF 943 IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN 944* The eigenvalue is computed to bisection accuracy 945* compute eigenvector and stop 946 USEDBS = .TRUE. 947 GOTO 120 948 ELSEIF( ITER.LT.MAXITR ) THEN 949 GOTO 120 950 ELSEIF( ITER.EQ.MAXITR ) THEN 951 NEEDBS = .TRUE. 952 GOTO 120 953 ELSE 954 INFO = 5 955 RETURN 956 END IF 957 ELSE 958 STP2II = .FALSE. 959 IF(USEDRQ .AND. USEDBS .AND. 960 $ BSTRES.LE.RESID) THEN 961 LAMBDA = BSTW 962 STP2II = .TRUE. 963 ENDIF 964 IF (STP2II) THEN 965* improve error angle by second step 966 CALL SLAR1V( IN, 1, IN, LAMBDA, 967 $ D( IBEGIN ), L( IBEGIN ), 968 $ WORK(INDLD+IBEGIN-1), 969 $ WORK(INDLLD+IBEGIN-1), 970 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 971 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 972 $ IWORK( IINDR+WINDEX ), 973 $ ISUPPZ( 2*WINDEX-1 ), 974 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 975 ENDIF 976 WORK( WINDEX ) = LAMBDA 977 END IF 978* 979* Compute FP-vector support w.r.t. whole matrix 980* 981 ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN 982 ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN 983 ZFROM = ISUPPZ( 2*WINDEX-1 ) 984 ZTO = ISUPPZ( 2*WINDEX ) 985 ISUPMN = ISUPMN + OLDIEN 986 ISUPMX = ISUPMX + OLDIEN 987* Ensure vector is ok if support in the RQI has changed 988 IF(ISUPMN.LT.ZFROM) THEN 989 DO 122 II = ISUPMN,ZFROM-1 990 Z( II, WINDEX ) = ZERO 991 122 CONTINUE 992 ENDIF 993 IF(ISUPMX.GT.ZTO) THEN 994 DO 123 II = ZTO+1,ISUPMX 995 Z( II, WINDEX ) = ZERO 996 123 CONTINUE 997 ENDIF 998 CALL SSCAL( ZTO-ZFROM+1, NRMINV, 999 $ Z( ZFROM, WINDEX ), 1 ) 1000 125 CONTINUE 1001* Update W 1002 W( WINDEX ) = LAMBDA+SIGMA 1003* Recompute the gaps on the left and right 1004* But only allow them to become larger and not 1005* smaller (which can only happen through "bad" 1006* cancellation and doesn't reflect the theory 1007* where the initial gaps are underestimated due 1008* to WERR being too crude.) 1009 IF(.NOT.ESKIP) THEN 1010 IF( K.GT.1) THEN 1011 WGAP( WINDMN ) = MAX( WGAP(WINDMN), 1012 $ W(WINDEX)-WERR(WINDEX) 1013 $ - W(WINDMN)-WERR(WINDMN) ) 1014 ENDIF 1015 IF( WINDEX.LT.WEND ) THEN 1016 WGAP( WINDEX ) = MAX( SAVGAP, 1017 $ W( WINDPL )-WERR( WINDPL ) 1018 $ - W( WINDEX )-WERR( WINDEX) ) 1019 ENDIF 1020 ENDIF 1021 IDONE = IDONE + 1 1022 ENDIF 1023* here ends the code for the current child 1024* 1025 139 CONTINUE 1026* Proceed to any remaining child nodes 1027 NEWFST = J + 1 1028 140 CONTINUE 1029 150 CONTINUE 1030 NDEPTH = NDEPTH + 1 1031 GO TO 40 1032 END IF 1033 IBEGIN = IEND + 1 1034 WBEGIN = WEND + 1 1035 170 CONTINUE 1036* 1037 1038 RETURN 1039* 1040* End of SLARRV 1041* 1042 END 1043