001: SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 002: $ LDVR, MM, M, WORK, RWORK, INFO ) 003: * 004: * -- LAPACK routine (version 3.2) -- 005: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 006: * November 2006 007: * 008: * .. Scalar Arguments .. 009: CHARACTER HOWMNY, SIDE 010: INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 011: * .. 012: * .. Array Arguments .. 013: LOGICAL SELECT( * ) 014: REAL RWORK( * ) 015: COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 016: $ WORK( * ) 017: * .. 018: * 019: * Purpose 020: * ======= 021: * 022: * CTREVC computes some or all of the right and/or left eigenvectors of 023: * a complex upper triangular matrix T. 024: * Matrices of this type are produced by the Schur factorization of 025: * a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR. 026: * 027: * The right eigenvector x and the left eigenvector y of T corresponding 028: * to an eigenvalue w are defined by: 029: * 030: * T*x = w*x, (y**H)*T = w*(y**H) 031: * 032: * where y**H denotes the conjugate transpose of the vector y. 033: * The eigenvalues are not input to this routine, but are read directly 034: * from the diagonal of T. 035: * 036: * This routine returns the matrices X and/or Y of right and left 037: * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 038: * input matrix. If Q is the unitary factor that reduces a matrix A to 039: * Schur form T, then Q*X and Q*Y are the matrices of right and left 040: * eigenvectors of A. 041: * 042: * Arguments 043: * ========= 044: * 045: * SIDE (input) CHARACTER*1 046: * = 'R': compute right eigenvectors only; 047: * = 'L': compute left eigenvectors only; 048: * = 'B': compute both right and left eigenvectors. 049: * 050: * HOWMNY (input) CHARACTER*1 051: * = 'A': compute all right and/or left eigenvectors; 052: * = 'B': compute all right and/or left eigenvectors, 053: * backtransformed using the matrices supplied in 054: * VR and/or VL; 055: * = 'S': compute selected right and/or left eigenvectors, 056: * as indicated by the logical array SELECT. 057: * 058: * SELECT (input) LOGICAL array, dimension (N) 059: * If HOWMNY = 'S', SELECT specifies the eigenvectors to be 060: * computed. 061: * The eigenvector corresponding to the j-th eigenvalue is 062: * computed if SELECT(j) = .TRUE.. 063: * Not referenced if HOWMNY = 'A' or 'B'. 064: * 065: * N (input) INTEGER 066: * The order of the matrix T. N >= 0. 067: * 068: * T (input/output) COMPLEX array, dimension (LDT,N) 069: * The upper triangular matrix T. T is modified, but restored 070: * on exit. 071: * 072: * LDT (input) INTEGER 073: * The leading dimension of the array T. LDT >= max(1,N). 074: * 075: * VL (input/output) COMPLEX array, dimension (LDVL,MM) 076: * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 077: * contain an N-by-N matrix Q (usually the unitary matrix Q of 078: * Schur vectors returned by CHSEQR). 079: * On exit, if SIDE = 'L' or 'B', VL contains: 080: * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 081: * if HOWMNY = 'B', the matrix Q*Y; 082: * if HOWMNY = 'S', the left eigenvectors of T specified by 083: * SELECT, stored consecutively in the columns 084: * of VL, in the same order as their 085: * eigenvalues. 086: * Not referenced if SIDE = 'R'. 087: * 088: * LDVL (input) INTEGER 089: * The leading dimension of the array VL. LDVL >= 1, and if 090: * SIDE = 'L' or 'B', LDVL >= N. 091: * 092: * VR (input/output) COMPLEX array, dimension (LDVR,MM) 093: * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 094: * contain an N-by-N matrix Q (usually the unitary matrix Q of 095: * Schur vectors returned by CHSEQR). 096: * On exit, if SIDE = 'R' or 'B', VR contains: 097: * if HOWMNY = 'A', the matrix X of right eigenvectors of T; 098: * if HOWMNY = 'B', the matrix Q*X; 099: * if HOWMNY = 'S', the right eigenvectors of T specified by 100: * SELECT, stored consecutively in the columns 101: * of VR, in the same order as their 102: * eigenvalues. 103: * Not referenced if SIDE = 'L'. 104: * 105: * LDVR (input) INTEGER 106: * The leading dimension of the array VR. LDVR >= 1, and if 107: * SIDE = 'R' or 'B'; LDVR >= N. 108: * 109: * MM (input) INTEGER 110: * The number of columns in the arrays VL and/or VR. MM >= M. 111: * 112: * M (output) INTEGER 113: * The number of columns in the arrays VL and/or VR actually 114: * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 115: * is set to N. Each selected eigenvector occupies one 116: * column. 117: * 118: * WORK (workspace) COMPLEX array, dimension (2*N) 119: * 120: * RWORK (workspace) REAL array, dimension (N) 121: * 122: * INFO (output) INTEGER 123: * = 0: successful exit 124: * < 0: if INFO = -i, the i-th argument had an illegal value 125: * 126: * Further Details 127: * =============== 128: * 129: * The algorithm used in this program is basically backward (forward) 130: * substitution, with scaling to make the the code robust against 131: * possible overflow. 132: * 133: * Each eigenvector is normalized so that the element of largest 134: * magnitude has magnitude 1; here the magnitude of a complex number 135: * (x,y) is taken to be |x| + |y|. 136: * 137: * ===================================================================== 138: * 139: * .. Parameters .. 140: REAL ZERO, ONE 141: PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 142: COMPLEX CMZERO, CMONE 143: PARAMETER ( CMZERO = ( 0.0E+0, 0.0E+0 ), 144: $ CMONE = ( 1.0E+0, 0.0E+0 ) ) 145: * .. 146: * .. Local Scalars .. 147: LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV 148: INTEGER I, II, IS, J, K, KI 149: REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL 150: COMPLEX CDUM 151: * .. 152: * .. External Functions .. 153: LOGICAL LSAME 154: INTEGER ICAMAX 155: REAL SCASUM, SLAMCH 156: EXTERNAL LSAME, ICAMAX, SCASUM, SLAMCH 157: * .. 158: * .. External Subroutines .. 159: EXTERNAL CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA 160: * .. 161: * .. Intrinsic Functions .. 162: INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, REAL 163: * .. 164: * .. Statement Functions .. 165: REAL CABS1 166: * .. 167: * .. Statement Function definitions .. 168: CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 169: * .. 170: * .. Executable Statements .. 171: * 172: * Decode and test the input parameters 173: * 174: BOTHV = LSAME( SIDE, 'B' ) 175: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 176: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 177: * 178: ALLV = LSAME( HOWMNY, 'A' ) 179: OVER = LSAME( HOWMNY, 'B' ) 180: SOMEV = LSAME( HOWMNY, 'S' ) 181: * 182: * Set M to the number of columns required to store the selected 183: * eigenvectors. 184: * 185: IF( SOMEV ) THEN 186: M = 0 187: DO 10 J = 1, N 188: IF( SELECT( J ) ) 189: $ M = M + 1 190: 10 CONTINUE 191: ELSE 192: M = N 193: END IF 194: * 195: INFO = 0 196: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 197: INFO = -1 198: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 199: INFO = -2 200: ELSE IF( N.LT.0 ) THEN 201: INFO = -4 202: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 203: INFO = -6 204: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 205: INFO = -8 206: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 207: INFO = -10 208: ELSE IF( MM.LT.M ) THEN 209: INFO = -11 210: END IF 211: IF( INFO.NE.0 ) THEN 212: CALL XERBLA( 'CTREVC', -INFO ) 213: RETURN 214: END IF 215: * 216: * Quick return if possible. 217: * 218: IF( N.EQ.0 ) 219: $ RETURN 220: * 221: * Set the constants to control overflow. 222: * 223: UNFL = SLAMCH( 'Safe minimum' ) 224: OVFL = ONE / UNFL 225: CALL SLABAD( UNFL, OVFL ) 226: ULP = SLAMCH( 'Precision' ) 227: SMLNUM = UNFL*( N / ULP ) 228: * 229: * Store the diagonal elements of T in working array WORK. 230: * 231: DO 20 I = 1, N 232: WORK( I+N ) = T( I, I ) 233: 20 CONTINUE 234: * 235: * Compute 1-norm of each column of strictly upper triangular 236: * part of T to control overflow in triangular solver. 237: * 238: RWORK( 1 ) = ZERO 239: DO 30 J = 2, N 240: RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 ) 241: 30 CONTINUE 242: * 243: IF( RIGHTV ) THEN 244: * 245: * Compute right eigenvectors. 246: * 247: IS = M 248: DO 80 KI = N, 1, -1 249: * 250: IF( SOMEV ) THEN 251: IF( .NOT.SELECT( KI ) ) 252: $ GO TO 80 253: END IF 254: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 255: * 256: WORK( 1 ) = CMONE 257: * 258: * Form right-hand side. 259: * 260: DO 40 K = 1, KI - 1 261: WORK( K ) = -T( K, KI ) 262: 40 CONTINUE 263: * 264: * Solve the triangular system: 265: * (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. 266: * 267: DO 50 K = 1, KI - 1 268: T( K, K ) = T( K, K ) - T( KI, KI ) 269: IF( CABS1( T( K, K ) ).LT.SMIN ) 270: $ T( K, K ) = SMIN 271: 50 CONTINUE 272: * 273: IF( KI.GT.1 ) THEN 274: CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', 275: $ KI-1, T, LDT, WORK( 1 ), SCALE, RWORK, 276: $ INFO ) 277: WORK( KI ) = SCALE 278: END IF 279: * 280: * Copy the vector x or Q*x to VR and normalize. 281: * 282: IF( .NOT.OVER ) THEN 283: CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 ) 284: * 285: II = ICAMAX( KI, VR( 1, IS ), 1 ) 286: REMAX = ONE / CABS1( VR( II, IS ) ) 287: CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 ) 288: * 289: DO 60 K = KI + 1, N 290: VR( K, IS ) = CMZERO 291: 60 CONTINUE 292: ELSE 293: IF( KI.GT.1 ) 294: $ CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ), 295: $ 1, CMPLX( SCALE ), VR( 1, KI ), 1 ) 296: * 297: II = ICAMAX( N, VR( 1, KI ), 1 ) 298: REMAX = ONE / CABS1( VR( II, KI ) ) 299: CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 ) 300: END IF 301: * 302: * Set back the original diagonal elements of T. 303: * 304: DO 70 K = 1, KI - 1 305: T( K, K ) = WORK( K+N ) 306: 70 CONTINUE 307: * 308: IS = IS - 1 309: 80 CONTINUE 310: END IF 311: * 312: IF( LEFTV ) THEN 313: * 314: * Compute left eigenvectors. 315: * 316: IS = 1 317: DO 130 KI = 1, N 318: * 319: IF( SOMEV ) THEN 320: IF( .NOT.SELECT( KI ) ) 321: $ GO TO 130 322: END IF 323: SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 324: * 325: WORK( N ) = CMONE 326: * 327: * Form right-hand side. 328: * 329: DO 90 K = KI + 1, N 330: WORK( K ) = -CONJG( T( KI, K ) ) 331: 90 CONTINUE 332: * 333: * Solve the triangular system: 334: * (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. 335: * 336: DO 100 K = KI + 1, N 337: T( K, K ) = T( K, K ) - T( KI, KI ) 338: IF( CABS1( T( K, K ) ).LT.SMIN ) 339: $ T( K, K ) = SMIN 340: 100 CONTINUE 341: * 342: IF( KI.LT.N ) THEN 343: CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 344: $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, 345: $ WORK( KI+1 ), SCALE, RWORK, INFO ) 346: WORK( KI ) = SCALE 347: END IF 348: * 349: * Copy the vector x or Q*x to VL and normalize. 350: * 351: IF( .NOT.OVER ) THEN 352: CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 ) 353: * 354: II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 355: REMAX = ONE / CABS1( VL( II, IS ) ) 356: CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 357: * 358: DO 110 K = 1, KI - 1 359: VL( K, IS ) = CMZERO 360: 110 CONTINUE 361: ELSE 362: IF( KI.LT.N ) 363: $ CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL, 364: $ WORK( KI+1 ), 1, CMPLX( SCALE ), 365: $ VL( 1, KI ), 1 ) 366: * 367: II = ICAMAX( N, VL( 1, KI ), 1 ) 368: REMAX = ONE / CABS1( VL( II, KI ) ) 369: CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 ) 370: END IF 371: * 372: * Set back the original diagonal elements of T. 373: * 374: DO 120 K = KI + 1, N 375: T( K, K ) = WORK( K+N ) 376: 120 CONTINUE 377: * 378: IS = IS + 1 379: 130 CONTINUE 380: END IF 381: * 382: RETURN 383: * 384: * End of CTREVC 385: * 386: END 387: