![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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