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