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