![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 2: $ LDVR, MM, M, WORK, RWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER HOWMNY, SIDE 11: INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 12: * .. 13: * .. Array Arguments .. 14: LOGICAL SELECT( * ) 15: DOUBLE PRECISION RWORK( * ) 16: COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 17: $ WORK( * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * ZTREVC computes some or all of the right and/or left eigenvectors of 24: * a complex upper triangular matrix T. 25: * Matrices of this type are produced by the Schur factorization of 26: * a complex general matrix: A = Q*T*Q**H, as computed by ZHSEQR. 27: * 28: * The right eigenvector x and the left eigenvector y of T corresponding 29: * to an eigenvalue w are defined by: 30: * 31: * T*x = w*x, (y**H)*T = w*(y**H) 32: * 33: * where y**H denotes the conjugate transpose of the vector y. 34: * The eigenvalues are not input to this routine, but are read directly 35: * from the diagonal of T. 36: * 37: * This routine returns the matrices X and/or Y of right and left 38: * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 39: * input matrix. If Q is the unitary factor that reduces a matrix A to 40: * Schur form T, then Q*X and Q*Y are the matrices of right and left 41: * eigenvectors of A. 42: * 43: * Arguments 44: * ========= 45: * 46: * SIDE (input) CHARACTER*1 47: * = 'R': compute right eigenvectors only; 48: * = 'L': compute left eigenvectors only; 49: * = 'B': compute both right and left eigenvectors. 50: * 51: * HOWMNY (input) CHARACTER*1 52: * = 'A': compute all right and/or left eigenvectors; 53: * = 'B': compute all right and/or left eigenvectors, 54: * backtransformed using the matrices supplied in 55: * VR and/or VL; 56: * = 'S': compute selected right and/or left eigenvectors, 57: * as indicated by the logical array SELECT. 58: * 59: * SELECT (input) LOGICAL array, dimension (N) 60: * If HOWMNY = 'S', SELECT specifies the eigenvectors to be 61: * computed. 62: * The eigenvector corresponding to the j-th eigenvalue is 63: * computed if SELECT(j) = .TRUE.. 64: * Not referenced if HOWMNY = 'A' or 'B'. 65: * 66: * N (input) INTEGER 67: * The order of the matrix T. N >= 0. 68: * 69: * T (input/output) COMPLEX*16 array, dimension (LDT,N) 70: * The upper triangular matrix T. T is modified, but restored 71: * on exit. 72: * 73: * LDT (input) INTEGER 74: * The leading dimension of the array T. LDT >= max(1,N). 75: * 76: * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) 77: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 78: * contain an N-by-N matrix Q (usually the unitary matrix Q of 79: * Schur vectors returned by ZHSEQR). 80: * On exit, if SIDE = 'L' or 'B', VL contains: 81: * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 82: * if HOWMNY = 'B', the matrix Q*Y; 83: * if HOWMNY = 'S', the left eigenvectors of T specified by 84: * SELECT, stored consecutively in the columns 85: * of VL, in the same order as their 86: * eigenvalues. 87: * Not referenced if SIDE = 'R'. 88: * 89: * LDVL (input) INTEGER 90: * The leading dimension of the array VL. LDVL >= 1, and if 91: * SIDE = 'L' or 'B', LDVL >= N. 92: * 93: * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) 94: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 95: * contain an N-by-N matrix Q (usually the unitary matrix Q of 96: * Schur vectors returned by ZHSEQR). 97: * On exit, if SIDE = 'R' or 'B', VR contains: 98: * if HOWMNY = 'A', the matrix X of right eigenvectors of T; 99: * if HOWMNY = 'B', the matrix Q*X; 100: * if HOWMNY = 'S', the right eigenvectors of T specified by 101: * SELECT, stored consecutively in the columns 102: * of VR, in the same order as their 103: * eigenvalues. 104: * Not referenced if SIDE = 'L'. 105: * 106: * LDVR (input) INTEGER 107: * The leading dimension of the array VR. LDVR >= 1, and if 108: * SIDE = 'R' or 'B'; LDVR >= N. 109: * 110: * MM (input) INTEGER 111: * The number of columns in the arrays VL and/or VR. MM >= M. 112: * 113: * M (output) INTEGER 114: * The number of columns in the arrays VL and/or VR actually 115: * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 116: * is set to N. Each selected eigenvector occupies one 117: * column. 118: * 119: * WORK (workspace) COMPLEX*16 array, dimension (2*N) 120: * 121: * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 122: * 123: * INFO (output) INTEGER 124: * = 0: successful exit 125: * < 0: if INFO = -i, the i-th argument had an illegal value 126: * 127: * Further Details 128: * =============== 129: * 130: * The algorithm used in this program is basically backward (forward) 131: * substitution, with scaling to make the the code robust against 132: * possible overflow. 133: * 134: * Each eigenvector is normalized so that the element of largest 135: * magnitude has magnitude 1; here the magnitude of a complex number 136: * (x,y) is taken to be |x| + |y|. 137: * 138: * ===================================================================== 139: * 140: * .. Parameters .. 141: DOUBLE PRECISION ZERO, ONE 142: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 143: COMPLEX*16 CMZERO, CMONE 144: PARAMETER ( CMZERO = ( 0.0D+0, 0.0D+0 ), 145: $ CMONE = ( 1.0D+0, 0.0D+0 ) ) 146: * .. 147: * .. Local Scalars .. 148: LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV 149: INTEGER I, II, IS, J, K, KI 150: DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL 151: COMPLEX*16 CDUM 152: * .. 153: * .. External Functions .. 154: LOGICAL LSAME 155: INTEGER IZAMAX 156: DOUBLE PRECISION DLAMCH, DZASUM 157: EXTERNAL LSAME, IZAMAX, DLAMCH, DZASUM 158: * .. 159: * .. External Subroutines .. 160: EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS 161: * .. 162: * .. Intrinsic Functions .. 163: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX 164: * .. 165: * .. Statement Functions .. 166: DOUBLE PRECISION CABS1 167: * .. 168: * .. Statement Function definitions .. 169: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 170: * .. 171: * .. Executable Statements .. 172: * 173: * Decode and test the input parameters 174: * 175: BOTHV = LSAME( SIDE, 'B' ) 176: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 177: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 178: * 179: ALLV = LSAME( HOWMNY, 'A' ) 180: OVER = LSAME( HOWMNY, 'B' ) 181: SOMEV = LSAME( HOWMNY, 'S' ) 182: * 183: * Set M to the number of columns required to store the selected 184: * eigenvectors. 185: * 186: IF( SOMEV ) THEN 187: M = 0 188: DO 10 J = 1, N 189: IF( SELECT( J ) ) 190: $ M = M + 1 191: 10 CONTINUE 192: ELSE 193: M = N 194: END IF 195: * 196: INFO = 0 197: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 198: INFO = -1 199: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 200: INFO = -2 201: ELSE IF( N.LT.0 ) THEN 202: INFO = -4 203: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 204: INFO = -6 205: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 206: INFO = -8 207: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 208: INFO = -10 209: ELSE IF( MM.LT.M ) THEN 210: INFO = -11 211: END IF 212: IF( INFO.NE.0 ) THEN 213: CALL XERBLA( 'ZTREVC', -INFO ) 214: RETURN 215: END IF 216: * 217: * Quick return if possible. 218: * 219: IF( N.EQ.0 ) 220: $ RETURN 221: * 222: * Set the constants to control overflow. 223: * 224: UNFL = DLAMCH( 'Safe minimum' ) 225: OVFL = ONE / UNFL 226: CALL DLABAD( UNFL, OVFL ) 227: ULP = DLAMCH( 'Precision' ) 228: SMLNUM = UNFL*( N / ULP ) 229: * 230: * Store the diagonal elements of T in working array WORK. 231: * 232: DO 20 I = 1, N 233: WORK( I+N ) = T( I, I ) 234: 20 CONTINUE 235: * 236: * Compute 1-norm of each column of strictly upper triangular 237: * part of T to control overflow in triangular solver. 238: * 239: RWORK( 1 ) = ZERO 240: DO 30 J = 2, N 241: RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 ) 242: 30 CONTINUE 243: * 244: IF( RIGHTV ) THEN 245: * 246: * Compute right eigenvectors. 247: * 248: IS = M 249: DO 80 KI = N, 1, -1 250: * 251: IF( SOMEV ) THEN 252: IF( .NOT.SELECT( KI ) ) 253: $ GO TO 80 254: END IF 255: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 256: * 257: WORK( 1 ) = CMONE 258: * 259: * Form right-hand side. 260: * 261: DO 40 K = 1, KI - 1 262: WORK( K ) = -T( K, KI ) 263: 40 CONTINUE 264: * 265: * Solve the triangular system: 266: * (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. 267: * 268: DO 50 K = 1, KI - 1 269: T( K, K ) = T( K, K ) - T( KI, KI ) 270: IF( CABS1( T( K, K ) ).LT.SMIN ) 271: $ T( K, K ) = SMIN 272: 50 CONTINUE 273: * 274: IF( KI.GT.1 ) THEN 275: CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', 276: $ KI-1, T, LDT, WORK( 1 ), SCALE, RWORK, 277: $ INFO ) 278: WORK( KI ) = SCALE 279: END IF 280: * 281: * Copy the vector x or Q*x to VR and normalize. 282: * 283: IF( .NOT.OVER ) THEN 284: CALL ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 ) 285: * 286: II = IZAMAX( KI, VR( 1, IS ), 1 ) 287: REMAX = ONE / CABS1( VR( II, IS ) ) 288: CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 ) 289: * 290: DO 60 K = KI + 1, N 291: VR( K, IS ) = CMZERO 292: 60 CONTINUE 293: ELSE 294: IF( KI.GT.1 ) 295: $ CALL ZGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ), 296: $ 1, DCMPLX( SCALE ), VR( 1, KI ), 1 ) 297: * 298: II = IZAMAX( N, VR( 1, KI ), 1 ) 299: REMAX = ONE / CABS1( VR( II, KI ) ) 300: CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 ) 301: END IF 302: * 303: * Set back the original diagonal elements of T. 304: * 305: DO 70 K = 1, KI - 1 306: T( K, K ) = WORK( K+N ) 307: 70 CONTINUE 308: * 309: IS = IS - 1 310: 80 CONTINUE 311: END IF 312: * 313: IF( LEFTV ) THEN 314: * 315: * Compute left eigenvectors. 316: * 317: IS = 1 318: DO 130 KI = 1, N 319: * 320: IF( SOMEV ) THEN 321: IF( .NOT.SELECT( KI ) ) 322: $ GO TO 130 323: END IF 324: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 325: * 326: WORK( N ) = CMONE 327: * 328: * Form right-hand side. 329: * 330: DO 90 K = KI + 1, N 331: WORK( K ) = -DCONJG( T( KI, K ) ) 332: 90 CONTINUE 333: * 334: * Solve the triangular system: 335: * (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. 336: * 337: DO 100 K = KI + 1, N 338: T( K, K ) = T( K, K ) - T( KI, KI ) 339: IF( CABS1( T( K, K ) ).LT.SMIN ) 340: $ T( K, K ) = SMIN 341: 100 CONTINUE 342: * 343: IF( KI.LT.N ) THEN 344: CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 345: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, 346: $ WORK( KI+1 ), SCALE, RWORK, INFO ) 347: WORK( KI ) = SCALE 348: END IF 349: * 350: * Copy the vector x or Q*x to VL and normalize. 351: * 352: IF( .NOT.OVER ) THEN 353: CALL ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 ) 354: * 355: II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 356: REMAX = ONE / CABS1( VL( II, IS ) ) 357: CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 358: * 359: DO 110 K = 1, KI - 1 360: VL( K, IS ) = CMZERO 361: 110 CONTINUE 362: ELSE 363: IF( KI.LT.N ) 364: $ CALL ZGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL, 365: $ WORK( KI+1 ), 1, DCMPLX( SCALE ), 366: $ VL( 1, KI ), 1 ) 367: * 368: II = IZAMAX( N, VL( 1, KI ), 1 ) 369: REMAX = ONE / CABS1( VL( II, KI ) ) 370: CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 ) 371: END IF 372: * 373: * Set back the original diagonal elements of T. 374: * 375: DO 120 K = KI + 1, N 376: T( K, K ) = WORK( K+N ) 377: 120 CONTINUE 378: * 379: IS = IS + 1 380: 130 CONTINUE 381: END IF 382: * 383: RETURN 384: * 385: * End of ZTREVC 386: * 387: END