Annotation of rpl/lapack/lapack/zlarrv.f, revision 1.1
1.1 ! bertrand 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) --
! 8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 10: * November 2006
! 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., th P-OFFSET
! 443: * thru' 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlarrv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack