Annotation of rpl/lapack/lapack/zlar1v.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
! 2: $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
! 3: $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
! 4: *
! 5: * -- LAPACK auxiliary routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * .. Scalar Arguments ..
! 11: LOGICAL WANTNC
! 12: INTEGER B1, BN, N, NEGCNT, R
! 13: DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
! 14: $ RQCORR, ZTZ
! 15: * ..
! 16: * .. Array Arguments ..
! 17: INTEGER ISUPPZ( * )
! 18: DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
! 19: $ WORK( * )
! 20: COMPLEX*16 Z( * )
! 21: * ..
! 22: *
! 23: * Purpose
! 24: * =======
! 25: *
! 26: * ZLAR1V computes the (scaled) r-th column of the inverse of
! 27: * the sumbmatrix in rows B1 through BN of the tridiagonal matrix
! 28: * L D L^T - sigma I. When sigma is close to an eigenvalue, the
! 29: * computed vector is an accurate eigenvector. Usually, r corresponds
! 30: * to the index where the eigenvector is largest in magnitude.
! 31: * The following steps accomplish this computation :
! 32: * (a) Stationary qd transform, L D L^T - sigma I = L(+) D(+) L(+)^T,
! 33: * (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
! 34: * (c) Computation of the diagonal elements of the inverse of
! 35: * L D L^T - sigma I by combining the above transforms, and choosing
! 36: * r as the index where the diagonal of the inverse is (one of the)
! 37: * largest in magnitude.
! 38: * (d) Computation of the (scaled) r-th column of the inverse using the
! 39: * twisted factorization obtained by combining the top part of the
! 40: * the stationary and the bottom part of the progressive transform.
! 41: *
! 42: * Arguments
! 43: * =========
! 44: *
! 45: * N (input) INTEGER
! 46: * The order of the matrix L D L^T.
! 47: *
! 48: * B1 (input) INTEGER
! 49: * First index of the submatrix of L D L^T.
! 50: *
! 51: * BN (input) INTEGER
! 52: * Last index of the submatrix of L D L^T.
! 53: *
! 54: * LAMBDA (input) DOUBLE PRECISION
! 55: * The shift. In order to compute an accurate eigenvector,
! 56: * LAMBDA should be a good approximation to an eigenvalue
! 57: * of L D L^T.
! 58: *
! 59: * L (input) DOUBLE PRECISION array, dimension (N-1)
! 60: * The (n-1) subdiagonal elements of the unit bidiagonal matrix
! 61: * L, in elements 1 to N-1.
! 62: *
! 63: * D (input) DOUBLE PRECISION array, dimension (N)
! 64: * The n diagonal elements of the diagonal matrix D.
! 65: *
! 66: * LD (input) DOUBLE PRECISION array, dimension (N-1)
! 67: * The n-1 elements L(i)*D(i).
! 68: *
! 69: * LLD (input) DOUBLE PRECISION array, dimension (N-1)
! 70: * The n-1 elements L(i)*L(i)*D(i).
! 71: *
! 72: * PIVMIN (input) DOUBLE PRECISION
! 73: * The minimum pivot in the Sturm sequence.
! 74: *
! 75: * GAPTOL (input) DOUBLE PRECISION
! 76: * Tolerance that indicates when eigenvector entries are negligible
! 77: * w.r.t. their contribution to the residual.
! 78: *
! 79: * Z (input/output) COMPLEX*16 array, dimension (N)
! 80: * On input, all entries of Z must be set to 0.
! 81: * On output, Z contains the (scaled) r-th column of the
! 82: * inverse. The scaling is such that Z(R) equals 1.
! 83: *
! 84: * WANTNC (input) LOGICAL
! 85: * Specifies whether NEGCNT has to be computed.
! 86: *
! 87: * NEGCNT (output) INTEGER
! 88: * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
! 89: * in the matrix factorization L D L^T, and NEGCNT = -1 otherwise.
! 90: *
! 91: * ZTZ (output) DOUBLE PRECISION
! 92: * The square of the 2-norm of Z.
! 93: *
! 94: * MINGMA (output) DOUBLE PRECISION
! 95: * The reciprocal of the largest (in magnitude) diagonal
! 96: * element of the inverse of L D L^T - sigma I.
! 97: *
! 98: * R (input/output) INTEGER
! 99: * The twist index for the twisted factorization used to
! 100: * compute Z.
! 101: * On input, 0 <= R <= N. If R is input as 0, R is set to
! 102: * the index where (L D L^T - sigma I)^{-1} is largest
! 103: * in magnitude. If 1 <= R <= N, R is unchanged.
! 104: * On output, R contains the twist index used to compute Z.
! 105: * Ideally, R designates the position of the maximum entry in the
! 106: * eigenvector.
! 107: *
! 108: * ISUPPZ (output) INTEGER array, dimension (2)
! 109: * The support of the vector in Z, i.e., the vector Z is
! 110: * nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
! 111: *
! 112: * NRMINV (output) DOUBLE PRECISION
! 113: * NRMINV = 1/SQRT( ZTZ )
! 114: *
! 115: * RESID (output) DOUBLE PRECISION
! 116: * The residual of the FP vector.
! 117: * RESID = ABS( MINGMA )/SQRT( ZTZ )
! 118: *
! 119: * RQCORR (output) DOUBLE PRECISION
! 120: * The Rayleigh Quotient correction to LAMBDA.
! 121: * RQCORR = MINGMA*TMP
! 122: *
! 123: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
! 124: *
! 125: * Further Details
! 126: * ===============
! 127: *
! 128: * Based on contributions by
! 129: * Beresford Parlett, University of California, Berkeley, USA
! 130: * Jim Demmel, University of California, Berkeley, USA
! 131: * Inderjit Dhillon, University of Texas, Austin, USA
! 132: * Osni Marques, LBNL/NERSC, USA
! 133: * Christof Voemel, University of California, Berkeley, USA
! 134: *
! 135: * =====================================================================
! 136: *
! 137: * .. Parameters ..
! 138: DOUBLE PRECISION ZERO, ONE
! 139: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 140: COMPLEX*16 CONE
! 141: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
! 142:
! 143: * ..
! 144: * .. Local Scalars ..
! 145: LOGICAL SAWNAN1, SAWNAN2
! 146: INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
! 147: $ R2
! 148: DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
! 149: * ..
! 150: * .. External Functions ..
! 151: LOGICAL DISNAN
! 152: DOUBLE PRECISION DLAMCH
! 153: EXTERNAL DISNAN, DLAMCH
! 154: * ..
! 155: * .. Intrinsic Functions ..
! 156: INTRINSIC ABS, DBLE
! 157: * ..
! 158: * .. Executable Statements ..
! 159: *
! 160: EPS = DLAMCH( 'Precision' )
! 161:
! 162:
! 163: IF( R.EQ.0 ) THEN
! 164: R1 = B1
! 165: R2 = BN
! 166: ELSE
! 167: R1 = R
! 168: R2 = R
! 169: END IF
! 170:
! 171: * Storage for LPLUS
! 172: INDLPL = 0
! 173: * Storage for UMINUS
! 174: INDUMN = N
! 175: INDS = 2*N + 1
! 176: INDP = 3*N + 1
! 177:
! 178: IF( B1.EQ.1 ) THEN
! 179: WORK( INDS ) = ZERO
! 180: ELSE
! 181: WORK( INDS+B1-1 ) = LLD( B1-1 )
! 182: END IF
! 183:
! 184: *
! 185: * Compute the stationary transform (using the differential form)
! 186: * until the index R2.
! 187: *
! 188: SAWNAN1 = .FALSE.
! 189: NEG1 = 0
! 190: S = WORK( INDS+B1-1 ) - LAMBDA
! 191: DO 50 I = B1, R1 - 1
! 192: DPLUS = D( I ) + S
! 193: WORK( INDLPL+I ) = LD( I ) / DPLUS
! 194: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
! 195: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
! 196: S = WORK( INDS+I ) - LAMBDA
! 197: 50 CONTINUE
! 198: SAWNAN1 = DISNAN( S )
! 199: IF( SAWNAN1 ) GOTO 60
! 200: DO 51 I = R1, R2 - 1
! 201: DPLUS = D( I ) + S
! 202: WORK( INDLPL+I ) = LD( I ) / DPLUS
! 203: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
! 204: S = WORK( INDS+I ) - LAMBDA
! 205: 51 CONTINUE
! 206: SAWNAN1 = DISNAN( S )
! 207: *
! 208: 60 CONTINUE
! 209: IF( SAWNAN1 ) THEN
! 210: * Runs a slower version of the above loop if a NaN is detected
! 211: NEG1 = 0
! 212: S = WORK( INDS+B1-1 ) - LAMBDA
! 213: DO 70 I = B1, R1 - 1
! 214: DPLUS = D( I ) + S
! 215: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
! 216: WORK( INDLPL+I ) = LD( I ) / DPLUS
! 217: IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
! 218: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
! 219: IF( WORK( INDLPL+I ).EQ.ZERO )
! 220: $ WORK( INDS+I ) = LLD( I )
! 221: S = WORK( INDS+I ) - LAMBDA
! 222: 70 CONTINUE
! 223: DO 71 I = R1, R2 - 1
! 224: DPLUS = D( I ) + S
! 225: IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
! 226: WORK( INDLPL+I ) = LD( I ) / DPLUS
! 227: WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
! 228: IF( WORK( INDLPL+I ).EQ.ZERO )
! 229: $ WORK( INDS+I ) = LLD( I )
! 230: S = WORK( INDS+I ) - LAMBDA
! 231: 71 CONTINUE
! 232: END IF
! 233: *
! 234: * Compute the progressive transform (using the differential form)
! 235: * until the index R1
! 236: *
! 237: SAWNAN2 = .FALSE.
! 238: NEG2 = 0
! 239: WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
! 240: DO 80 I = BN - 1, R1, -1
! 241: DMINUS = LLD( I ) + WORK( INDP+I )
! 242: TMP = D( I ) / DMINUS
! 243: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
! 244: WORK( INDUMN+I ) = L( I )*TMP
! 245: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
! 246: 80 CONTINUE
! 247: TMP = WORK( INDP+R1-1 )
! 248: SAWNAN2 = DISNAN( TMP )
! 249:
! 250: IF( SAWNAN2 ) THEN
! 251: * Runs a slower version of the above loop if a NaN is detected
! 252: NEG2 = 0
! 253: DO 100 I = BN-1, R1, -1
! 254: DMINUS = LLD( I ) + WORK( INDP+I )
! 255: IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
! 256: TMP = D( I ) / DMINUS
! 257: IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
! 258: WORK( INDUMN+I ) = L( I )*TMP
! 259: WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
! 260: IF( TMP.EQ.ZERO )
! 261: $ WORK( INDP+I-1 ) = D( I ) - LAMBDA
! 262: 100 CONTINUE
! 263: END IF
! 264: *
! 265: * Find the index (from R1 to R2) of the largest (in magnitude)
! 266: * diagonal element of the inverse
! 267: *
! 268: MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
! 269: IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
! 270: IF( WANTNC ) THEN
! 271: NEGCNT = NEG1 + NEG2
! 272: ELSE
! 273: NEGCNT = -1
! 274: ENDIF
! 275: IF( ABS(MINGMA).EQ.ZERO )
! 276: $ MINGMA = EPS*WORK( INDS+R1-1 )
! 277: R = R1
! 278: DO 110 I = R1, R2 - 1
! 279: TMP = WORK( INDS+I ) + WORK( INDP+I )
! 280: IF( TMP.EQ.ZERO )
! 281: $ TMP = EPS*WORK( INDS+I )
! 282: IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
! 283: MINGMA = TMP
! 284: R = I + 1
! 285: END IF
! 286: 110 CONTINUE
! 287: *
! 288: * Compute the FP vector: solve N^T v = e_r
! 289: *
! 290: ISUPPZ( 1 ) = B1
! 291: ISUPPZ( 2 ) = BN
! 292: Z( R ) = CONE
! 293: ZTZ = ONE
! 294: *
! 295: * Compute the FP vector upwards from R
! 296: *
! 297: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
! 298: DO 210 I = R-1, B1, -1
! 299: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
! 300: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
! 301: $ THEN
! 302: Z( I ) = ZERO
! 303: ISUPPZ( 1 ) = I + 1
! 304: GOTO 220
! 305: ENDIF
! 306: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
! 307: 210 CONTINUE
! 308: 220 CONTINUE
! 309: ELSE
! 310: * Run slower loop if NaN occurred.
! 311: DO 230 I = R - 1, B1, -1
! 312: IF( Z( I+1 ).EQ.ZERO ) THEN
! 313: Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
! 314: ELSE
! 315: Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
! 316: END IF
! 317: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
! 318: $ THEN
! 319: Z( I ) = ZERO
! 320: ISUPPZ( 1 ) = I + 1
! 321: GO TO 240
! 322: END IF
! 323: ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
! 324: 230 CONTINUE
! 325: 240 CONTINUE
! 326: ENDIF
! 327:
! 328: * Compute the FP vector downwards from R in blocks of size BLKSIZ
! 329: IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
! 330: DO 250 I = R, BN-1
! 331: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
! 332: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
! 333: $ THEN
! 334: Z( I+1 ) = ZERO
! 335: ISUPPZ( 2 ) = I
! 336: GO TO 260
! 337: END IF
! 338: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
! 339: 250 CONTINUE
! 340: 260 CONTINUE
! 341: ELSE
! 342: * Run slower loop if NaN occurred.
! 343: DO 270 I = R, BN - 1
! 344: IF( Z( I ).EQ.ZERO ) THEN
! 345: Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
! 346: ELSE
! 347: Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
! 348: END IF
! 349: IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
! 350: $ THEN
! 351: Z( I+1 ) = ZERO
! 352: ISUPPZ( 2 ) = I
! 353: GO TO 280
! 354: END IF
! 355: ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
! 356: 270 CONTINUE
! 357: 280 CONTINUE
! 358: END IF
! 359: *
! 360: * Compute quantities for convergence test
! 361: *
! 362: TMP = ONE / ZTZ
! 363: NRMINV = SQRT( TMP )
! 364: RESID = ABS( MINGMA )*NRMINV
! 365: RQCORR = MINGMA*TMP
! 366: *
! 367: *
! 368: RETURN
! 369: *
! 370: * End of ZLAR1V
! 371: *
! 372: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlar1v.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack