 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00002 $ LDVR, MM, M, WORK, RWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER HOWMNY, SIDE 00011 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 00012 * .. 00013 * .. Array Arguments .. 00014 LOGICAL SELECT( * ) 00015 REAL RWORK( * ) 00016 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00017 $ WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * CTREVC computes some or all of the right and/or left eigenvectors of 00024 * a complex upper triangular matrix T. 00025 * Matrices of this type are produced by the Schur factorization of 00026 * a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR. 00027 * 00028 * The right eigenvector x and the left eigenvector y of T corresponding 00029 * to an eigenvalue w are defined by: 00030 * 00031 * T*x = w*x, (y**H)*T = w*(y**H) 00032 * 00033 * where y**H denotes the conjugate transpose of the vector y. 00034 * The eigenvalues are not input to this routine, but are read directly 00035 * from the diagonal of T. 00036 * 00037 * This routine returns the matrices X and/or Y of right and left 00038 * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 00039 * input matrix. If Q is the unitary factor that reduces a matrix A to 00040 * Schur form T, then Q*X and Q*Y are the matrices of right and left 00041 * eigenvectors of A. 00042 * 00043 * Arguments 00044 * ========= 00045 * 00046 * SIDE (input) CHARACTER*1 00047 * = 'R': compute right eigenvectors only; 00048 * = 'L': compute left eigenvectors only; 00049 * = 'B': compute both right and left eigenvectors. 00050 * 00051 * HOWMNY (input) CHARACTER*1 00052 * = 'A': compute all right and/or left eigenvectors; 00053 * = 'B': compute all right and/or left eigenvectors, 00054 * backtransformed using the matrices supplied in 00055 * VR and/or VL; 00056 * = 'S': compute selected right and/or left eigenvectors, 00057 * as indicated by the logical array SELECT. 00058 * 00059 * SELECT (input) LOGICAL array, dimension (N) 00060 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be 00061 * computed. 00062 * The eigenvector corresponding to the j-th eigenvalue is 00063 * computed if SELECT(j) = .TRUE.. 00064 * Not referenced if HOWMNY = 'A' or 'B'. 00065 * 00066 * N (input) INTEGER 00067 * The order of the matrix T. N >= 0. 00068 * 00069 * T (input/output) COMPLEX array, dimension (LDT,N) 00070 * The upper triangular matrix T. T is modified, but restored 00071 * on exit. 00072 * 00073 * LDT (input) INTEGER 00074 * The leading dimension of the array T. LDT >= max(1,N). 00075 * 00076 * VL (input/output) COMPLEX array, dimension (LDVL,MM) 00077 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00078 * contain an N-by-N matrix Q (usually the unitary matrix Q of 00079 * Schur vectors returned by CHSEQR). 00080 * On exit, if SIDE = 'L' or 'B', VL contains: 00081 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 00082 * if HOWMNY = 'B', the matrix Q*Y; 00083 * if HOWMNY = 'S', the left eigenvectors of T specified by 00084 * SELECT, stored consecutively in the columns 00085 * of VL, in the same order as their 00086 * eigenvalues. 00087 * Not referenced if SIDE = 'R'. 00088 * 00089 * LDVL (input) INTEGER 00090 * The leading dimension of the array VL. LDVL >= 1, and if 00091 * SIDE = 'L' or 'B', LDVL >= N. 00092 * 00093 * VR (input/output) COMPLEX array, dimension (LDVR,MM) 00094 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00095 * contain an N-by-N matrix Q (usually the unitary matrix Q of 00096 * Schur vectors returned by CHSEQR). 00097 * On exit, if SIDE = 'R' or 'B', VR contains: 00098 * if HOWMNY = 'A', the matrix X of right eigenvectors of T; 00099 * if HOWMNY = 'B', the matrix Q*X; 00100 * if HOWMNY = 'S', the right eigenvectors of T specified by 00101 * SELECT, stored consecutively in the columns 00102 * of VR, in the same order as their 00103 * eigenvalues. 00104 * Not referenced if SIDE = 'L'. 00105 * 00106 * LDVR (input) INTEGER 00107 * The leading dimension of the array VR. LDVR >= 1, and if 00108 * SIDE = 'R' or 'B'; LDVR >= N. 00109 * 00110 * MM (input) INTEGER 00111 * The number of columns in the arrays VL and/or VR. MM >= M. 00112 * 00113 * M (output) INTEGER 00114 * The number of columns in the arrays VL and/or VR actually 00115 * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 00116 * is set to N. Each selected eigenvector occupies one 00117 * column. 00118 * 00119 * WORK (workspace) COMPLEX array, dimension (2*N) 00120 * 00121 * RWORK (workspace) REAL array, dimension (N) 00122 * 00123 * INFO (output) INTEGER 00124 * = 0: successful exit 00125 * < 0: if INFO = -i, the i-th argument had an illegal value 00126 * 00127 * Further Details 00128 * =============== 00129 * 00130 * The algorithm used in this program is basically backward (forward) 00131 * substitution, with scaling to make the the code robust against 00132 * possible overflow. 00133 * 00134 * Each eigenvector is normalized so that the element of largest 00135 * magnitude has magnitude 1; here the magnitude of a complex number 00136 * (x,y) is taken to be |x| + |y|. 00137 * 00138 * ===================================================================== 00139 * 00140 * .. Parameters .. 00141 REAL ZERO, ONE 00142 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00143 COMPLEX CMZERO, CMONE 00144 PARAMETER ( CMZERO = ( 0.0E+0, 0.0E+0 ), 00145 $ CMONE = ( 1.0E+0, 0.0E+0 ) ) 00146 * .. 00147 * .. Local Scalars .. 00148 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV 00149 INTEGER I, II, IS, J, K, KI 00150 REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL 00151 COMPLEX CDUM 00152 * .. 00153 * .. External Functions .. 00154 LOGICAL LSAME 00155 INTEGER ICAMAX 00156 REAL SCASUM, SLAMCH 00157 EXTERNAL LSAME, ICAMAX, SCASUM, SLAMCH 00158 * .. 00159 * .. External Subroutines .. 00160 EXTERNAL CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA 00161 * .. 00162 * .. Intrinsic Functions .. 00163 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, REAL 00164 * .. 00165 * .. Statement Functions .. 00166 REAL CABS1 00167 * .. 00168 * .. Statement Function definitions .. 00169 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00170 * .. 00171 * .. Executable Statements .. 00172 * 00173 * Decode and test the input parameters 00174 * 00175 BOTHV = LSAME( SIDE, 'B' ) 00176 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 00177 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 00178 * 00179 ALLV = LSAME( HOWMNY, 'A' ) 00180 OVER = LSAME( HOWMNY, 'B' ) 00181 SOMEV = LSAME( HOWMNY, 'S' ) 00182 * 00183 * Set M to the number of columns required to store the selected 00184 * eigenvectors. 00185 * 00186 IF( SOMEV ) THEN 00187 M = 0 00188 DO 10 J = 1, N 00189 IF( SELECT( J ) ) 00190 $ M = M + 1 00191 10 CONTINUE 00192 ELSE 00193 M = N 00194 END IF 00195 * 00196 INFO = 0 00197 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 00198 INFO = -1 00199 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 00200 INFO = -2 00201 ELSE IF( N.LT.0 ) THEN 00202 INFO = -4 00203 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00204 INFO = -6 00205 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 00206 INFO = -8 00207 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 00208 INFO = -10 00209 ELSE IF( MM.LT.M ) THEN 00210 INFO = -11 00211 END IF 00212 IF( INFO.NE.0 ) THEN 00213 CALL XERBLA( 'CTREVC', -INFO ) 00214 RETURN 00215 END IF 00216 * 00217 * Quick return if possible. 00218 * 00219 IF( N.EQ.0 ) 00220 $ RETURN 00221 * 00222 * Set the constants to control overflow. 00223 * 00224 UNFL = SLAMCH( 'Safe minimum' ) 00225 OVFL = ONE / UNFL 00226 CALL SLABAD( UNFL, OVFL ) 00227 ULP = SLAMCH( 'Precision' ) 00228 SMLNUM = UNFL*( N / ULP ) 00229 * 00230 * Store the diagonal elements of T in working array WORK. 00231 * 00232 DO 20 I = 1, N 00233 WORK( I+N ) = T( I, I ) 00234 20 CONTINUE 00235 * 00236 * Compute 1-norm of each column of strictly upper triangular 00237 * part of T to control overflow in triangular solver. 00238 * 00239 RWORK( 1 ) = ZERO 00240 DO 30 J = 2, N 00241 RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 ) 00242 30 CONTINUE 00243 * 00244 IF( RIGHTV ) THEN 00245 * 00246 * Compute right eigenvectors. 00247 * 00248 IS = M 00249 DO 80 KI = N, 1, -1 00250 * 00251 IF( SOMEV ) THEN 00252 IF( .NOT.SELECT( KI ) ) 00253 $ GO TO 80 00254 END IF 00255 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 00256 * 00257 WORK( 1 ) = CMONE 00258 * 00259 * Form right-hand side. 00260 * 00261 DO 40 K = 1, KI - 1 00262 WORK( K ) = -T( K, KI ) 00263 40 CONTINUE 00264 * 00265 * Solve the triangular system: 00266 * (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. 00267 * 00268 DO 50 K = 1, KI - 1 00269 T( K, K ) = T( K, K ) - T( KI, KI ) 00270 IF( CABS1( T( K, K ) ).LT.SMIN ) 00271 $ T( K, K ) = SMIN 00272 50 CONTINUE 00273 * 00274 IF( KI.GT.1 ) THEN 00275 CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', 00276 $ KI-1, T, LDT, WORK( 1 ), SCALE, RWORK, 00277 $ INFO ) 00278 WORK( KI ) = SCALE 00279 END IF 00280 * 00281 * Copy the vector x or Q*x to VR and normalize. 00282 * 00283 IF( .NOT.OVER ) THEN 00284 CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 ) 00285 * 00286 II = ICAMAX( KI, VR( 1, IS ), 1 ) 00287 REMAX = ONE / CABS1( VR( II, IS ) ) 00288 CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00289 * 00290 DO 60 K = KI + 1, N 00291 VR( K, IS ) = CMZERO 00292 60 CONTINUE 00293 ELSE 00294 IF( KI.GT.1 ) 00295 $ CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ), 00296 $ 1, CMPLX( SCALE ), VR( 1, KI ), 1 ) 00297 * 00298 II = ICAMAX( N, VR( 1, KI ), 1 ) 00299 REMAX = ONE / CABS1( VR( II, KI ) ) 00300 CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 ) 00301 END IF 00302 * 00303 * Set back the original diagonal elements of T. 00304 * 00305 DO 70 K = 1, KI - 1 00306 T( K, K ) = WORK( K+N ) 00307 70 CONTINUE 00308 * 00309 IS = IS - 1 00310 80 CONTINUE 00311 END IF 00312 * 00313 IF( LEFTV ) THEN 00314 * 00315 * Compute left eigenvectors. 00316 * 00317 IS = 1 00318 DO 130 KI = 1, N 00319 * 00320 IF( SOMEV ) THEN 00321 IF( .NOT.SELECT( KI ) ) 00322 $ GO TO 130 00323 END IF 00324 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 00325 * 00326 WORK( N ) = CMONE 00327 * 00328 * Form right-hand side. 00329 * 00330 DO 90 K = KI + 1, N 00331 WORK( K ) = -CONJG( T( KI, K ) ) 00332 90 CONTINUE 00333 * 00334 * Solve the triangular system: 00335 * (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK. 00336 * 00337 DO 100 K = KI + 1, N 00338 T( K, K ) = T( K, K ) - T( KI, KI ) 00339 IF( CABS1( T( K, K ) ).LT.SMIN ) 00340 $ T( K, K ) = SMIN 00341 100 CONTINUE 00342 * 00343 IF( KI.LT.N ) THEN 00344 CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 00345 $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, 00346 $ WORK( KI+1 ), SCALE, RWORK, INFO ) 00347 WORK( KI ) = SCALE 00348 END IF 00349 * 00350 * Copy the vector x or Q*x to VL and normalize. 00351 * 00352 IF( .NOT.OVER ) THEN 00353 CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 ) 00354 * 00355 II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 00356 REMAX = ONE / CABS1( VL( II, IS ) ) 00357 CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00358 * 00359 DO 110 K = 1, KI - 1 00360 VL( K, IS ) = CMZERO 00361 110 CONTINUE 00362 ELSE 00363 IF( KI.LT.N ) 00364 $ CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL, 00365 $ WORK( KI+1 ), 1, CMPLX( SCALE ), 00366 $ VL( 1, KI ), 1 ) 00367 * 00368 II = ICAMAX( N, VL( 1, KI ), 1 ) 00369 REMAX = ONE / CABS1( VL( II, KI ) ) 00370 CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 ) 00371 END IF 00372 * 00373 * Set back the original diagonal elements of T. 00374 * 00375 DO 120 K = KI + 1, N 00376 T( K, K ) = WORK( K+N ) 00377 120 CONTINUE 00378 * 00379 IS = IS + 1 00380 130 CONTINUE 00381 END IF 00382 * 00383 RETURN 00384 * 00385 * End of CTREVC 00386 * 00387 END
1.7.3 
