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