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*> \date June 2016 270* 271*> \ingroup complex16OTHERauxiliary 272* 273*> \par Contributors: 274* ================== 275*> 276*> Beresford Parlett, University of California, Berkeley, USA \n 277*> Jim Demmel, University of California, Berkeley, USA \n 278*> Inderjit Dhillon, University of Texas, Austin, USA \n 279*> Osni Marques, LBNL/NERSC, USA \n 280*> Christof Voemel, University of California, Berkeley, USA 281* 282* ===================================================================== 283 SUBROUTINE ZLARRV( N, VL, VU, D, L, PIVMIN, 284 $ ISPLIT, M, DOL, DOU, MINRGP, 285 $ RTOL1, RTOL2, W, WERR, WGAP, 286 $ IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 287 $ WORK, IWORK, INFO ) 288* 289* -- LAPACK auxiliary routine (version 3.7.1) -- 290* -- LAPACK is a software package provided by Univ. of Tennessee, -- 291* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 292* June 2016 293* 294* .. Scalar Arguments .. 295 INTEGER DOL, DOU, INFO, LDZ, M, N 296 DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 297* .. 298* .. Array Arguments .. 299 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 300 $ ISUPPZ( * ), IWORK( * ) 301 DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 302 $ WGAP( * ), WORK( * ) 303 COMPLEX*16 Z( LDZ, * ) 304* .. 305* 306* ===================================================================== 307* 308* .. Parameters .. 309 INTEGER MAXITR 310 PARAMETER ( MAXITR = 10 ) 311 COMPLEX*16 CZERO 312 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) ) 313 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, HALF 314 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 315 $ TWO = 2.0D0, THREE = 3.0D0, 316 $ FOUR = 4.0D0, HALF = 0.5D0) 317* .. 318* .. Local Scalars .. 319 LOGICAL ESKIP, NEEDBS, STP2II, TRYRQC, USEDBS, USEDRQ 320 INTEGER DONE, I, IBEGIN, IDONE, IEND, II, IINDC1, 321 $ IINDC2, IINDR, IINDWK, IINFO, IM, IN, INDEIG, 322 $ INDLD, INDLLD, INDWRK, ISUPMN, ISUPMX, ITER, 323 $ ITMP1, J, JBLK, K, MINIWSIZE, MINWSIZE, NCLUS, 324 $ NDEPTH, NEGCNT, NEWCLS, NEWFST, NEWFTT, NEWLST, 325 $ NEWSIZ, OFFSET, OLDCLS, OLDFST, OLDIEN, OLDLST, 326 $ OLDNCL, P, PARITY, Q, WBEGIN, WEND, WINDEX, 327 $ WINDMN, WINDPL, ZFROM, ZTO, ZUSEDL, ZUSEDU, 328 $ ZUSEDW 329 INTEGER INDIN1, INDIN2 330 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 337 EXTERNAL DLAMCH 338* .. 339* .. External Subroutines .. 340 EXTERNAL DCOPY, DLARRB, DLARRF, ZDSCAL, ZLAR1V, 341 $ ZLASET 342* .. 343* .. Intrinsic Functions .. 344 INTRINSIC ABS, DBLE, MAX, MIN 345 INTRINSIC DCMPLX 346* .. 347* .. Executable Statements .. 348* .. 349 350 INFO = 0 351* 352* Quick return if possible 353* 354 IF( (N.LE.0).OR.(M.LE.0) ) THEN 355 RETURN 356 END IF 357* 358* The first N entries of WORK are reserved for the eigenvalues 359 INDLD = N+1 360 INDLLD= 2*N+1 361 INDIN1 = 3*N + 1 362 INDIN2 = 4*N + 1 363 INDWRK = 5*N + 1 364 MINWSIZE = 12 * N 365 366 DO 5 I= 1,MINWSIZE 367 WORK( I ) = ZERO 368 5 CONTINUE 369 370* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the 371* factorization used to compute the FP vector 372 IINDR = 0 373* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current 374* layer and the one above. 375 IINDC1 = N 376 IINDC2 = 2*N 377 IINDWK = 3*N + 1 378 379 MINIWSIZE = 7 * N 380 DO 10 I= 1,MINIWSIZE 381 IWORK( I ) = 0 382 10 CONTINUE 383 384 ZUSEDL = 1 385 IF(DOL.GT.1) THEN 386* Set lower bound for use of Z 387 ZUSEDL = DOL-1 388 ENDIF 389 ZUSEDU = M 390 IF(DOU.LT.M) THEN 391* Set lower bound for use of Z 392 ZUSEDU = DOU+1 393 ENDIF 394* The width of the part of Z that is used 395 ZUSEDW = ZUSEDU - ZUSEDL + 1 396 397 398 CALL ZLASET( 'Full', N, ZUSEDW, CZERO, CZERO, 399 $ Z(1,ZUSEDL), LDZ ) 400 401 EPS = DLAMCH( 'Precision' ) 402 RQTOL = TWO * EPS 403* 404* Set expert flags for standard code. 405 TRYRQC = .TRUE. 406 407 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 408 ELSE 409* Only selected eigenpairs are computed. Since the other evalues 410* are not refined by RQ iteration, bisection has to compute to full 411* accuracy. 412 RTOL1 = FOUR * EPS 413 RTOL2 = FOUR * EPS 414 ENDIF 415 416* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the 417* desired eigenvalues. The support of the nonzero eigenvector 418* entries is contained in the interval IBEGIN:IEND. 419* Remark that if k eigenpairs are desired, then the eigenvectors 420* are stored in k contiguous columns of Z. 421 422* DONE is the number of eigenvectors already computed 423 DONE = 0 424 IBEGIN = 1 425 WBEGIN = 1 426 DO 170 JBLK = 1, IBLOCK( M ) 427 IEND = ISPLIT( JBLK ) 428 SIGMA = L( IEND ) 429* Find the eigenvectors of the submatrix indexed IBEGIN 430* through IEND. 431 WEND = WBEGIN - 1 432 15 CONTINUE 433 IF( WEND.LT.M ) THEN 434 IF( IBLOCK( WEND+1 ).EQ.JBLK ) THEN 435 WEND = WEND + 1 436 GO TO 15 437 END IF 438 END IF 439 IF( WEND.LT.WBEGIN ) THEN 440 IBEGIN = IEND + 1 441 GO TO 170 442 ELSEIF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN 443 IBEGIN = IEND + 1 444 WBEGIN = WEND + 1 445 GO TO 170 446 END IF 447 448* Find local spectral diameter of the block 449 GL = GERS( 2*IBEGIN-1 ) 450 GU = GERS( 2*IBEGIN ) 451 DO 20 I = IBEGIN+1 , IEND 452 GL = MIN( GERS( 2*I-1 ), GL ) 453 GU = MAX( GERS( 2*I ), GU ) 454 20 CONTINUE 455 SPDIAM = GU - GL 456 457* OLDIEN is the last index of the previous block 458 OLDIEN = IBEGIN - 1 459* Calculate the size of the current block 460 IN = IEND - IBEGIN + 1 461* The number of eigenvalues in the current block 462 IM = WEND - WBEGIN + 1 463 464* This is for a 1x1 block 465 IF( IBEGIN.EQ.IEND ) THEN 466 DONE = DONE+1 467 Z( IBEGIN, WBEGIN ) = DCMPLX( ONE, ZERO ) 468 ISUPPZ( 2*WBEGIN-1 ) = IBEGIN 469 ISUPPZ( 2*WBEGIN ) = IBEGIN 470 W( WBEGIN ) = W( WBEGIN ) + SIGMA 471 WORK( WBEGIN ) = W( WBEGIN ) 472 IBEGIN = IEND + 1 473 WBEGIN = WBEGIN + 1 474 GO TO 170 475 END IF 476 477* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) 478* Note that these can be approximations, in this case, the corresp. 479* entries of WERR give the size of the uncertainty interval. 480* The eigenvalue approximations will be refined when necessary as 481* high relative accuracy is required for the computation of the 482* corresponding eigenvectors. 483 CALL DCOPY( IM, W( WBEGIN ), 1, 484 $ WORK( WBEGIN ), 1 ) 485 486* We store in W the eigenvalue approximations w.r.t. the original 487* matrix T. 488 DO 30 I=1,IM 489 W(WBEGIN+I-1) = W(WBEGIN+I-1)+SIGMA 490 30 CONTINUE 491 492 493* NDEPTH is the current depth of the representation tree 494 NDEPTH = 0 495* PARITY is either 1 or 0 496 PARITY = 1 497* NCLUS is the number of clusters for the next level of the 498* representation tree, we start with NCLUS = 1 for the root 499 NCLUS = 1 500 IWORK( IINDC1+1 ) = 1 501 IWORK( IINDC1+2 ) = IM 502 503* IDONE is the number of eigenvectors already computed in the current 504* block 505 IDONE = 0 506* loop while( IDONE.LT.IM ) 507* generate the representation tree for the current block and 508* compute the eigenvectors 509 40 CONTINUE 510 IF( IDONE.LT.IM ) THEN 511* This is a crude protection against infinitely deep trees 512 IF( NDEPTH.GT.M ) THEN 513 INFO = -2 514 RETURN 515 ENDIF 516* breadth first processing of the current level of the representation 517* tree: OLDNCL = number of clusters on current level 518 OLDNCL = NCLUS 519* reset NCLUS to count the number of child clusters 520 NCLUS = 0 521* 522 PARITY = 1 - PARITY 523 IF( PARITY.EQ.0 ) THEN 524 OLDCLS = IINDC1 525 NEWCLS = IINDC2 526 ELSE 527 OLDCLS = IINDC2 528 NEWCLS = IINDC1 529 END IF 530* Process the clusters on the current level 531 DO 150 I = 1, OLDNCL 532 J = OLDCLS + 2*I 533* OLDFST, OLDLST = first, last index of current cluster. 534* cluster indices start with 1 and are relative 535* to WBEGIN when accessing W, WGAP, WERR, Z 536 OLDFST = IWORK( J-1 ) 537 OLDLST = IWORK( J ) 538 IF( NDEPTH.GT.0 ) THEN 539* Retrieve relatively robust representation (RRR) of cluster 540* that has been computed at the previous level 541* The RRR is stored in Z and overwritten once the eigenvectors 542* have been computed or when the cluster is refined 543 544 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 545* Get representation from location of the leftmost evalue 546* of the cluster 547 J = WBEGIN + OLDFST - 1 548 ELSE 549 IF(WBEGIN+OLDFST-1.LT.DOL) THEN 550* Get representation from the left end of Z array 551 J = DOL - 1 552 ELSEIF(WBEGIN+OLDFST-1.GT.DOU) THEN 553* Get representation from the right end of Z array 554 J = DOU 555 ELSE 556 J = WBEGIN + OLDFST - 1 557 ENDIF 558 ENDIF 559 DO 45 K = 1, IN - 1 560 D( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1, 561 $ J ) ) 562 L( IBEGIN+K-1 ) = DBLE( Z( IBEGIN+K-1, 563 $ J+1 ) ) 564 45 CONTINUE 565 D( IEND ) = DBLE( Z( IEND, J ) ) 566 SIGMA = DBLE( Z( IEND, J+1 ) ) 567 568* Set the corresponding entries in Z to zero 569 CALL ZLASET( 'Full', IN, 2, CZERO, CZERO, 570 $ Z( IBEGIN, J), LDZ ) 571 END IF 572 573* Compute DL and DLL of current RRR 574 DO 50 J = IBEGIN, IEND-1 575 TMP = D( J )*L( J ) 576 WORK( INDLD-1+J ) = TMP 577 WORK( INDLLD-1+J ) = TMP*L( J ) 578 50 CONTINUE 579 580 IF( NDEPTH.GT.0 ) THEN 581* P and Q are index of the first and last eigenvalue to compute 582* within the current block 583 P = INDEXW( WBEGIN-1+OLDFST ) 584 Q = INDEXW( WBEGIN-1+OLDLST ) 585* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 586* through the Q-OFFSET elements of these arrays are to be used. 587* OFFSET = P-OLDFST 588 OFFSET = INDEXW( WBEGIN ) - 1 589* perform limited bisection (if necessary) to get approximate 590* eigenvalues to the precision needed. 591 CALL DLARRB( IN, D( IBEGIN ), 592 $ WORK(INDLLD+IBEGIN-1), 593 $ P, Q, RTOL1, RTOL2, OFFSET, 594 $ WORK(WBEGIN),WGAP(WBEGIN),WERR(WBEGIN), 595 $ WORK( INDWRK ), IWORK( IINDWK ), 596 $ PIVMIN, SPDIAM, IN, IINFO ) 597 IF( IINFO.NE.0 ) THEN 598 INFO = -1 599 RETURN 600 ENDIF 601* We also recompute the extremal gaps. W holds all eigenvalues 602* of the unshifted matrix and must be used for computation 603* of WGAP, the entries of WORK might stem from RRRs with 604* different shifts. The gaps from WBEGIN-1+OLDFST to 605* WBEGIN-1+OLDLST are correctly computed in DLARRB. 606* However, we only allow the gaps to become greater since 607* this is what should happen when we decrease WERR 608 IF( OLDFST.GT.1) THEN 609 WGAP( WBEGIN+OLDFST-2 ) = 610 $ MAX(WGAP(WBEGIN+OLDFST-2), 611 $ W(WBEGIN+OLDFST-1)-WERR(WBEGIN+OLDFST-1) 612 $ - W(WBEGIN+OLDFST-2)-WERR(WBEGIN+OLDFST-2) ) 613 ENDIF 614 IF( WBEGIN + OLDLST -1 .LT. WEND ) THEN 615 WGAP( WBEGIN+OLDLST-1 ) = 616 $ MAX(WGAP(WBEGIN+OLDLST-1), 617 $ W(WBEGIN+OLDLST)-WERR(WBEGIN+OLDLST) 618 $ - W(WBEGIN+OLDLST-1)-WERR(WBEGIN+OLDLST-1) ) 619 ENDIF 620* Each time the eigenvalues in WORK get refined, we store 621* the newly found approximation with all shifts applied in W 622 DO 53 J=OLDFST,OLDLST 623 W(WBEGIN+J-1) = WORK(WBEGIN+J-1)+SIGMA 624 53 CONTINUE 625 END IF 626 627* Process the current node. 628 NEWFST = OLDFST 629 DO 140 J = OLDFST, OLDLST 630 IF( J.EQ.OLDLST ) THEN 631* we are at the right end of the cluster, this is also the 632* boundary of the child cluster 633 NEWLST = J 634 ELSE IF ( WGAP( WBEGIN + J -1).GE. 635 $ MINRGP* ABS( WORK(WBEGIN + J -1) ) ) THEN 636* the right relative gap is big enough, the child cluster 637* (NEWFST,..,NEWLST) is well separated from the following 638 NEWLST = J 639 ELSE 640* inside a child cluster, the relative gap is not 641* big enough. 642 GOTO 140 643 END IF 644 645* Compute size of child cluster found 646 NEWSIZ = NEWLST - NEWFST + 1 647 648* NEWFTT is the place in Z where the new RRR or the computed 649* eigenvector is to be stored 650 IF((DOL.EQ.1).AND.(DOU.EQ.M)) THEN 651* Store representation at location of the leftmost evalue 652* of the cluster 653 NEWFTT = WBEGIN + NEWFST - 1 654 ELSE 655 IF(WBEGIN+NEWFST-1.LT.DOL) THEN 656* Store representation at the left end of Z array 657 NEWFTT = DOL - 1 658 ELSEIF(WBEGIN+NEWFST-1.GT.DOU) THEN 659* Store representation at the right end of Z array 660 NEWFTT = DOU 661 ELSE 662 NEWFTT = WBEGIN + NEWFST - 1 663 ENDIF 664 ENDIF 665 666 IF( NEWSIZ.GT.1) THEN 667* 668* Current child is not a singleton but a cluster. 669* Compute and store new representation of child. 670* 671* 672* Compute left and right cluster gap. 673* 674* LGAP and RGAP are not computed from WORK because 675* the eigenvalue approximations may stem from RRRs 676* different shifts. However, W hold all eigenvalues 677* of the unshifted matrix. Still, the entries in WGAP 678* have to be computed from WORK since the entries 679* in W might be of the same order so that gaps are not 680* exhibited correctly for very close eigenvalues. 681 IF( NEWFST.EQ.1 ) THEN 682 LGAP = MAX( ZERO, 683 $ W(WBEGIN)-WERR(WBEGIN) - VL ) 684 ELSE 685 LGAP = WGAP( WBEGIN+NEWFST-2 ) 686 ENDIF 687 RGAP = WGAP( WBEGIN+NEWLST-1 ) 688* 689* Compute left- and rightmost eigenvalue of child 690* to high precision in order to shift as close 691* as possible and obtain as large relative gaps 692* as possible 693* 694 DO 55 K =1,2 695 IF(K.EQ.1) THEN 696 P = INDEXW( WBEGIN-1+NEWFST ) 697 ELSE 698 P = INDEXW( WBEGIN-1+NEWLST ) 699 ENDIF 700 OFFSET = INDEXW( WBEGIN ) - 1 701 CALL DLARRB( IN, D(IBEGIN), 702 $ WORK( INDLLD+IBEGIN-1 ),P,P, 703 $ RQTOL, RQTOL, OFFSET, 704 $ WORK(WBEGIN),WGAP(WBEGIN), 705 $ WERR(WBEGIN),WORK( INDWRK ), 706 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 707 $ IN, IINFO ) 708 55 CONTINUE 709* 710 IF((WBEGIN+NEWLST-1.LT.DOL).OR. 711 $ (WBEGIN+NEWFST-1.GT.DOU)) THEN 712* if the cluster contains no desired eigenvalues 713* skip the computation of that branch of the rep. tree 714* 715* We could skip before the refinement of the extremal 716* eigenvalues of the child, but then the representation 717* tree could be different from the one when nothing is 718* skipped. For this reason we skip at this place. 719 IDONE = IDONE + NEWLST - NEWFST + 1 720 GOTO 139 721 ENDIF 722* 723* Compute RRR of child cluster. 724* Note that the new RRR is stored in Z 725* 726* DLARRF needs LWORK = 2*N 727 CALL DLARRF( IN, D( IBEGIN ), L( IBEGIN ), 728 $ WORK(INDLD+IBEGIN-1), 729 $ NEWFST, NEWLST, WORK(WBEGIN), 730 $ WGAP(WBEGIN), WERR(WBEGIN), 731 $ SPDIAM, LGAP, RGAP, PIVMIN, TAU, 732 $ WORK( INDIN1 ), WORK( INDIN2 ), 733 $ WORK( INDWRK ), IINFO ) 734* In the complex case, DLARRF cannot write 735* the new RRR directly into Z and needs an intermediate 736* workspace 737 DO 56 K = 1, IN-1 738 Z( IBEGIN+K-1, NEWFTT ) = 739 $ DCMPLX( WORK( INDIN1+K-1 ), ZERO ) 740 Z( IBEGIN+K-1, NEWFTT+1 ) = 741 $ DCMPLX( WORK( INDIN2+K-1 ), ZERO ) 742 56 CONTINUE 743 Z( IEND, NEWFTT ) = 744 $ DCMPLX( WORK( INDIN1+IN-1 ), ZERO ) 745 IF( IINFO.EQ.0 ) THEN 746* a new RRR for the cluster was found by DLARRF 747* update shift and store it 748 SSIGMA = SIGMA + TAU 749 Z( IEND, NEWFTT+1 ) = DCMPLX( SSIGMA, ZERO ) 750* WORK() are the midpoints and WERR() the semi-width 751* Note that the entries in W are unchanged. 752 DO 116 K = NEWFST, NEWLST 753 FUDGE = 754 $ THREE*EPS*ABS(WORK(WBEGIN+K-1)) 755 WORK( WBEGIN + K - 1 ) = 756 $ WORK( WBEGIN + K - 1) - TAU 757 FUDGE = FUDGE + 758 $ FOUR*EPS*ABS(WORK(WBEGIN+K-1)) 759* Fudge errors 760 WERR( WBEGIN + K - 1 ) = 761 $ WERR( WBEGIN + K - 1 ) + FUDGE 762* Gaps are not fudged. Provided that WERR is small 763* when eigenvalues are close, a zero gap indicates 764* that a new representation is needed for resolving 765* the cluster. A fudge could lead to a wrong decision 766* of judging eigenvalues 'separated' which in 767* reality are not. This could have a negative impact 768* on the orthogonality of the computed eigenvectors. 769 116 CONTINUE 770 771 NCLUS = NCLUS + 1 772 K = NEWCLS + 2*NCLUS 773 IWORK( K-1 ) = NEWFST 774 IWORK( K ) = NEWLST 775 ELSE 776 INFO = -2 777 RETURN 778 ENDIF 779 ELSE 780* 781* Compute eigenvector of singleton 782* 783 ITER = 0 784* 785 TOL = FOUR * LOG(DBLE(IN)) * EPS 786* 787 K = NEWFST 788 WINDEX = WBEGIN + K - 1 789 WINDMN = MAX(WINDEX - 1,1) 790 WINDPL = MIN(WINDEX + 1,M) 791 LAMBDA = WORK( WINDEX ) 792 DONE = DONE + 1 793* Check if eigenvector computation is to be skipped 794 IF((WINDEX.LT.DOL).OR. 795 $ (WINDEX.GT.DOU)) THEN 796 ESKIP = .TRUE. 797 GOTO 125 798 ELSE 799 ESKIP = .FALSE. 800 ENDIF 801 LEFT = WORK( WINDEX ) - WERR( WINDEX ) 802 RIGHT = WORK( WINDEX ) + WERR( WINDEX ) 803 INDEIG = INDEXW( WINDEX ) 804* Note that since we compute the eigenpairs for a child, 805* all eigenvalue approximations are w.r.t the same shift. 806* In this case, the entries in WORK should be used for 807* computing the gaps since they exhibit even very small 808* differences in the eigenvalues, as opposed to the 809* entries in W which might "look" the same. 810 811 IF( K .EQ. 1) THEN 812* In the case RANGE='I' and with not much initial 813* accuracy in LAMBDA and VL, the formula 814* LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 815* can lead to an overestimation of the left gap and 816* thus to inadequately early RQI 'convergence'. 817* Prevent this by forcing a small left gap. 818 LGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 819 ELSE 820 LGAP = WGAP(WINDMN) 821 ENDIF 822 IF( K .EQ. IM) THEN 823* In the case RANGE='I' and with not much initial 824* accuracy in LAMBDA and VU, the formula 825* can lead to an overestimation of the right gap and 826* thus to inadequately early RQI 'convergence'. 827* Prevent this by forcing a small right gap. 828 RGAP = EPS*MAX(ABS(LEFT),ABS(RIGHT)) 829 ELSE 830 RGAP = WGAP(WINDEX) 831 ENDIF 832 GAP = MIN( LGAP, RGAP ) 833 IF(( K .EQ. 1).OR.(K .EQ. IM)) THEN 834* The eigenvector support can become wrong 835* because significant entries could be cut off due to a 836* large GAPTOL parameter in LAR1V. Prevent this. 837 GAPTOL = ZERO 838 ELSE 839 GAPTOL = GAP * EPS 840 ENDIF 841 ISUPMN = IN 842 ISUPMX = 1 843* Update WGAP so that it holds the minimum gap 844* to the left or the right. This is crucial in the 845* case where bisection is used to ensure that the 846* eigenvalue is refined up to the required precision. 847* The correct value is restored afterwards. 848 SAVGAP = WGAP(WINDEX) 849 WGAP(WINDEX) = GAP 850* We want to use the Rayleigh Quotient Correction 851* as often as possible since it converges quadratically 852* when we are close enough to the desired eigenvalue. 853* However, the Rayleigh Quotient can have the wrong sign 854* and lead us away from the desired eigenvalue. In this 855* case, the best we can do is to use bisection. 856 USEDBS = .FALSE. 857 USEDRQ = .FALSE. 858* Bisection is initially turned off unless it is forced 859 NEEDBS = .NOT.TRYRQC 860 120 CONTINUE 861* Check if bisection should be used to refine eigenvalue 862 IF(NEEDBS) THEN 863* Take the bisection as new iterate 864 USEDBS = .TRUE. 865 ITMP1 = IWORK( IINDR+WINDEX ) 866 OFFSET = INDEXW( WBEGIN ) - 1 867 CALL DLARRB( IN, D(IBEGIN), 868 $ WORK(INDLLD+IBEGIN-1),INDEIG,INDEIG, 869 $ ZERO, TWO*EPS, OFFSET, 870 $ WORK(WBEGIN),WGAP(WBEGIN), 871 $ WERR(WBEGIN),WORK( INDWRK ), 872 $ IWORK( IINDWK ), PIVMIN, SPDIAM, 873 $ ITMP1, IINFO ) 874 IF( IINFO.NE.0 ) THEN 875 INFO = -3 876 RETURN 877 ENDIF 878 LAMBDA = WORK( WINDEX ) 879* Reset twist index from inaccurate LAMBDA to 880* force computation of true MINGMA 881 IWORK( IINDR+WINDEX ) = 0 882 ENDIF 883* Given LAMBDA, compute the eigenvector. 884 CALL ZLAR1V( IN, 1, IN, LAMBDA, D( IBEGIN ), 885 $ L( IBEGIN ), WORK(INDLD+IBEGIN-1), 886 $ WORK(INDLLD+IBEGIN-1), 887 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 888 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 889 $ IWORK( IINDR+WINDEX ), ISUPPZ( 2*WINDEX-1 ), 890 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 891 IF(ITER .EQ. 0) THEN 892 BSTRES = RESID 893 BSTW = LAMBDA 894 ELSEIF(RESID.LT.BSTRES) THEN 895 BSTRES = RESID 896 BSTW = LAMBDA 897 ENDIF 898 ISUPMN = MIN(ISUPMN,ISUPPZ( 2*WINDEX-1 )) 899 ISUPMX = MAX(ISUPMX,ISUPPZ( 2*WINDEX )) 900 ITER = ITER + 1 901 902* sin alpha <= |resid|/gap 903* Note that both the residual and the gap are 904* proportional to the matrix, so ||T|| doesn't play 905* a role in the quotient 906 907* 908* Convergence test for Rayleigh-Quotient iteration 909* (omitted when Bisection has been used) 910* 911 IF( RESID.GT.TOL*GAP .AND. ABS( RQCORR ).GT. 912 $ RQTOL*ABS( LAMBDA ) .AND. .NOT. USEDBS) 913 $ THEN 914* We need to check that the RQCORR update doesn't 915* move the eigenvalue away from the desired one and 916* towards a neighbor. -> protection with bisection 917 IF(INDEIG.LE.NEGCNT) THEN 918* The wanted eigenvalue lies to the left 919 SGNDEF = -ONE 920 ELSE 921* The wanted eigenvalue lies to the right 922 SGNDEF = ONE 923 ENDIF 924* We only use the RQCORR if it improves the 925* the iterate reasonably. 926 IF( ( RQCORR*SGNDEF.GE.ZERO ) 927 $ .AND.( LAMBDA + RQCORR.LE. RIGHT) 928 $ .AND.( LAMBDA + RQCORR.GE. LEFT) 929 $ ) THEN 930 USEDRQ = .TRUE. 931* Store new midpoint of bisection interval in WORK 932 IF(SGNDEF.EQ.ONE) THEN 933* The current LAMBDA is on the left of the true 934* eigenvalue 935 LEFT = LAMBDA 936* We prefer to assume that the error estimate 937* is correct. We could make the interval not 938* as a bracket but to be modified if the RQCORR 939* chooses to. In this case, the RIGHT side should 940* be modified as follows: 941* RIGHT = MAX(RIGHT, LAMBDA + RQCORR) 942 ELSE 943* The current LAMBDA is on the right of the true 944* eigenvalue 945 RIGHT = LAMBDA 946* See comment about assuming the error estimate is 947* correct above. 948* LEFT = MIN(LEFT, LAMBDA + RQCORR) 949 ENDIF 950 WORK( WINDEX ) = 951 $ HALF * (RIGHT + LEFT) 952* Take RQCORR since it has the correct sign and 953* improves the iterate reasonably 954 LAMBDA = LAMBDA + RQCORR 955* Update width of error interval 956 WERR( WINDEX ) = 957 $ HALF * (RIGHT-LEFT) 958 ELSE 959 NEEDBS = .TRUE. 960 ENDIF 961 IF(RIGHT-LEFT.LT.RQTOL*ABS(LAMBDA)) THEN 962* The eigenvalue is computed to bisection accuracy 963* compute eigenvector and stop 964 USEDBS = .TRUE. 965 GOTO 120 966 ELSEIF( ITER.LT.MAXITR ) THEN 967 GOTO 120 968 ELSEIF( ITER.EQ.MAXITR ) THEN 969 NEEDBS = .TRUE. 970 GOTO 120 971 ELSE 972 INFO = 5 973 RETURN 974 END IF 975 ELSE 976 STP2II = .FALSE. 977 IF(USEDRQ .AND. USEDBS .AND. 978 $ BSTRES.LE.RESID) THEN 979 LAMBDA = BSTW 980 STP2II = .TRUE. 981 ENDIF 982 IF (STP2II) THEN 983* improve error angle by second step 984 CALL ZLAR1V( IN, 1, IN, LAMBDA, 985 $ D( IBEGIN ), L( IBEGIN ), 986 $ WORK(INDLD+IBEGIN-1), 987 $ WORK(INDLLD+IBEGIN-1), 988 $ PIVMIN, GAPTOL, Z( IBEGIN, WINDEX ), 989 $ .NOT.USEDBS, NEGCNT, ZTZ, MINGMA, 990 $ IWORK( IINDR+WINDEX ), 991 $ ISUPPZ( 2*WINDEX-1 ), 992 $ NRMINV, RESID, RQCORR, WORK( INDWRK ) ) 993 ENDIF 994 WORK( WINDEX ) = LAMBDA 995 END IF 996* 997* Compute FP-vector support w.r.t. whole matrix 998* 999 ISUPPZ( 2*WINDEX-1 ) = ISUPPZ( 2*WINDEX-1 )+OLDIEN 1000 ISUPPZ( 2*WINDEX ) = ISUPPZ( 2*WINDEX )+OLDIEN 1001 ZFROM = ISUPPZ( 2*WINDEX-1 ) 1002 ZTO = ISUPPZ( 2*WINDEX ) 1003 ISUPMN = ISUPMN + OLDIEN 1004 ISUPMX = ISUPMX + OLDIEN 1005* Ensure vector is ok if support in the RQI has changed 1006 IF(ISUPMN.LT.ZFROM) THEN 1007 DO 122 II = ISUPMN,ZFROM-1 1008 Z( II, WINDEX ) = ZERO 1009 122 CONTINUE 1010 ENDIF 1011 IF(ISUPMX.GT.ZTO) THEN 1012 DO 123 II = ZTO+1,ISUPMX 1013 Z( II, WINDEX ) = ZERO 1014 123 CONTINUE 1015 ENDIF 1016 CALL ZDSCAL( ZTO-ZFROM+1, NRMINV, 1017 $ Z( ZFROM, WINDEX ), 1 ) 1018 125 CONTINUE 1019* Update W 1020 W( WINDEX ) = LAMBDA+SIGMA 1021* Recompute the gaps on the left and right 1022* But only allow them to become larger and not 1023* smaller (which can only happen through "bad" 1024* cancellation and doesn't reflect the theory 1025* where the initial gaps are underestimated due 1026* to WERR being too crude.) 1027 IF(.NOT.ESKIP) THEN 1028 IF( K.GT.1) THEN 1029 WGAP( WINDMN ) = MAX( WGAP(WINDMN), 1030 $ W(WINDEX)-WERR(WINDEX) 1031 $ - W(WINDMN)-WERR(WINDMN) ) 1032 ENDIF 1033 IF( WINDEX.LT.WEND ) THEN 1034 WGAP( WINDEX ) = MAX( SAVGAP, 1035 $ W( WINDPL )-WERR( WINDPL ) 1036 $ - W( WINDEX )-WERR( WINDEX) ) 1037 ENDIF 1038 ENDIF 1039 IDONE = IDONE + 1 1040 ENDIF 1041* here ends the code for the current child 1042* 1043 139 CONTINUE 1044* Proceed to any remaining child nodes 1045 NEWFST = J + 1 1046 140 CONTINUE 1047 150 CONTINUE 1048 NDEPTH = NDEPTH + 1 1049 GO TO 40 1050 END IF 1051 IBEGIN = IEND + 1 1052 WBEGIN = WEND + 1 1053 170 CONTINUE 1054* 1055 1056 RETURN 1057* 1058* End of ZLARRV 1059* 1060 END 1061