1*> \brief \b CTREVC 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CTREVC + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrevc.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrevc.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrevc.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 22* LDVR, MM, M, WORK, RWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER HOWMNY, SIDE 26* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 27* .. 28* .. Array Arguments .. 29* LOGICAL SELECT( * ) 30* REAL RWORK( * ) 31* COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 32* $ WORK( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> CTREVC computes some or all of the right and/or left eigenvectors of 42*> a complex upper triangular matrix T. 43*> Matrices of this type are produced by the Schur factorization of 44*> a complex general matrix: A = Q*T*Q**H, as computed by CHSEQR. 45*> 46*> The right eigenvector x and the left eigenvector y of T corresponding 47*> to an eigenvalue w are defined by: 48*> 49*> T*x = w*x, (y**H)*T = w*(y**H) 50*> 51*> where y**H denotes the conjugate transpose of the vector y. 52*> The eigenvalues are not input to this routine, but are read directly 53*> from the diagonal of T. 54*> 55*> This routine returns the matrices X and/or Y of right and left 56*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 57*> input matrix. If Q is the unitary factor that reduces a matrix A to 58*> Schur form T, then Q*X and Q*Y are the matrices of right and left 59*> eigenvectors of A. 60*> \endverbatim 61* 62* Arguments: 63* ========== 64* 65*> \param[in] SIDE 66*> \verbatim 67*> SIDE is CHARACTER*1 68*> = 'R': compute right eigenvectors only; 69*> = 'L': compute left eigenvectors only; 70*> = 'B': compute both right and left eigenvectors. 71*> \endverbatim 72*> 73*> \param[in] HOWMNY 74*> \verbatim 75*> HOWMNY is CHARACTER*1 76*> = 'A': compute all right and/or left eigenvectors; 77*> = 'B': compute all right and/or left eigenvectors, 78*> backtransformed using the matrices supplied in 79*> VR and/or VL; 80*> = 'S': compute selected right and/or left eigenvectors, 81*> as indicated by the logical array SELECT. 82*> \endverbatim 83*> 84*> \param[in] SELECT 85*> \verbatim 86*> SELECT is LOGICAL array, dimension (N) 87*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be 88*> computed. 89*> The eigenvector corresponding to the j-th eigenvalue is 90*> computed if SELECT(j) = .TRUE.. 91*> Not referenced if HOWMNY = 'A' or 'B'. 92*> \endverbatim 93*> 94*> \param[in] N 95*> \verbatim 96*> N is INTEGER 97*> The order of the matrix T. N >= 0. 98*> \endverbatim 99*> 100*> \param[in,out] T 101*> \verbatim 102*> T is COMPLEX array, dimension (LDT,N) 103*> The upper triangular matrix T. T is modified, but restored 104*> on exit. 105*> \endverbatim 106*> 107*> \param[in] LDT 108*> \verbatim 109*> LDT is INTEGER 110*> The leading dimension of the array T. LDT >= max(1,N). 111*> \endverbatim 112*> 113*> \param[in,out] VL 114*> \verbatim 115*> VL is COMPLEX array, dimension (LDVL,MM) 116*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 117*> contain an N-by-N matrix Q (usually the unitary matrix Q of 118*> Schur vectors returned by CHSEQR). 119*> On exit, if SIDE = 'L' or 'B', VL contains: 120*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 121*> if HOWMNY = 'B', the matrix Q*Y; 122*> if HOWMNY = 'S', the left eigenvectors of T specified by 123*> SELECT, stored consecutively in the columns 124*> of VL, in the same order as their 125*> eigenvalues. 126*> Not referenced if SIDE = 'R'. 127*> \endverbatim 128*> 129*> \param[in] LDVL 130*> \verbatim 131*> LDVL is INTEGER 132*> The leading dimension of the array VL. LDVL >= 1, and if 133*> SIDE = 'L' or 'B', LDVL >= N. 134*> \endverbatim 135*> 136*> \param[in,out] VR 137*> \verbatim 138*> VR is COMPLEX array, dimension (LDVR,MM) 139*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 140*> contain an N-by-N matrix Q (usually the unitary matrix Q of 141*> Schur vectors returned by CHSEQR). 142*> On exit, if SIDE = 'R' or 'B', VR contains: 143*> if HOWMNY = 'A', the matrix X of right eigenvectors of T; 144*> if HOWMNY = 'B', the matrix Q*X; 145*> if HOWMNY = 'S', the right eigenvectors of T specified by 146*> SELECT, stored consecutively in the columns 147*> of VR, in the same order as their 148*> eigenvalues. 149*> Not referenced if SIDE = 'L'. 150*> \endverbatim 151*> 152*> \param[in] LDVR 153*> \verbatim 154*> LDVR is INTEGER 155*> The leading dimension of the array VR. LDVR >= 1, and if 156*> SIDE = 'R' or 'B'; LDVR >= N. 157*> \endverbatim 158*> 159*> \param[in] MM 160*> \verbatim 161*> MM is INTEGER 162*> The number of columns in the arrays VL and/or VR. MM >= M. 163*> \endverbatim 164*> 165*> \param[out] M 166*> \verbatim 167*> M is INTEGER 168*> The number of columns in the arrays VL and/or VR actually 169*> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 170*> is set to N. Each selected eigenvector occupies one 171*> column. 172*> \endverbatim 173*> 174*> \param[out] WORK 175*> \verbatim 176*> WORK is COMPLEX array, dimension (2*N) 177*> \endverbatim 178*> 179*> \param[out] RWORK 180*> \verbatim 181*> RWORK is REAL array, dimension (N) 182*> \endverbatim 183*> 184*> \param[out] INFO 185*> \verbatim 186*> INFO is INTEGER 187*> = 0: successful exit 188*> < 0: if INFO = -i, the i-th argument had an illegal value 189*> \endverbatim 190* 191* Authors: 192* ======== 193* 194*> \author Univ. of Tennessee 195*> \author Univ. of California Berkeley 196*> \author Univ. of Colorado Denver 197*> \author NAG Ltd. 198* 199*> \date November 2011 200* 201*> \ingroup complexOTHERcomputational 202* 203*> \par Further Details: 204* ===================== 205*> 206*> \verbatim 207*> 208*> The algorithm used in this program is basically backward (forward) 209*> substitution, with scaling to make the the code robust against 210*> possible overflow. 211*> 212*> Each eigenvector is normalized so that the element of largest 213*> magnitude has magnitude 1; here the magnitude of a complex number 214*> (x,y) is taken to be |x| + |y|. 215*> \endverbatim 216*> 217* ===================================================================== 218 SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 219 $ LDVR, MM, M, WORK, RWORK, INFO ) 220* 221* -- LAPACK computational routine (version 3.4.0) -- 222* -- LAPACK is a software package provided by Univ. of Tennessee, -- 223* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 224* November 2011 225* 226* .. Scalar Arguments .. 227 CHARACTER HOWMNY, SIDE 228 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 229* .. 230* .. Array Arguments .. 231 LOGICAL SELECT( * ) 232 REAL RWORK( * ) 233 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 234 $ WORK( * ) 235* .. 236* 237* ===================================================================== 238* 239* .. Parameters .. 240 REAL ZERO, ONE 241 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 242 COMPLEX CMZERO, CMONE 243 PARAMETER ( CMZERO = ( 0.0E+0, 0.0E+0 ), 244 $ CMONE = ( 1.0E+0, 0.0E+0 ) ) 245* .. 246* .. Local Scalars .. 247 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV 248 INTEGER I, II, IS, J, K, KI 249 REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL 250 COMPLEX CDUM 251* .. 252* .. External Functions .. 253 LOGICAL LSAME 254 INTEGER ICAMAX 255 REAL SCASUM, SLAMCH 256 EXTERNAL LSAME, ICAMAX, SCASUM, SLAMCH 257* .. 258* .. External Subroutines .. 259 EXTERNAL CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA 260* .. 261* .. Intrinsic Functions .. 262 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, REAL 263* .. 264* .. Statement Functions .. 265 REAL CABS1 266* .. 267* .. Statement Function definitions .. 268 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 269* .. 270* .. Executable Statements .. 271* 272* Decode and test the input parameters 273* 274 BOTHV = LSAME( SIDE, 'B' ) 275 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 276 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 277* 278 ALLV = LSAME( HOWMNY, 'A' ) 279 OVER = LSAME( HOWMNY, 'B' ) 280 SOMEV = LSAME( HOWMNY, 'S' ) 281* 282* Set M to the number of columns required to store the selected 283* eigenvectors. 284* 285 IF( SOMEV ) THEN 286 M = 0 287 DO 10 J = 1, N 288 IF( SELECT( J ) ) 289 $ M = M + 1 290 10 CONTINUE 291 ELSE 292 M = N 293 END IF 294* 295 INFO = 0 296 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 297 INFO = -1 298 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 299 INFO = -2 300 ELSE IF( N.LT.0 ) THEN 301 INFO = -4 302 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 303 INFO = -6 304 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 305 INFO = -8 306 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 307 INFO = -10 308 ELSE IF( MM.LT.M ) THEN 309 INFO = -11 310 END IF 311 IF( INFO.NE.0 ) THEN 312 CALL XERBLA( 'CTREVC', -INFO ) 313 RETURN 314 END IF 315* 316* Quick return if possible. 317* 318 IF( N.EQ.0 ) 319 $ RETURN 320* 321* Set the constants to control overflow. 322* 323 UNFL = SLAMCH( 'Safe minimum' ) 324 OVFL = ONE / UNFL 325 CALL SLABAD( UNFL, OVFL ) 326 ULP = SLAMCH( 'Precision' ) 327 SMLNUM = UNFL*( N / ULP ) 328* 329* Store the diagonal elements of T in working array WORK. 330* 331 DO 20 I = 1, N 332 WORK( I+N ) = T( I, I ) 333 20 CONTINUE 334* 335* Compute 1-norm of each column of strictly upper triangular 336* part of T to control overflow in triangular solver. 337* 338 RWORK( 1 ) = ZERO 339 DO 30 J = 2, N 340 RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 ) 341 30 CONTINUE 342* 343 IF( RIGHTV ) THEN 344* 345* Compute right eigenvectors. 346* 347 IS = M 348 DO 80 KI = N, 1, -1 349* 350 IF( SOMEV ) THEN 351 IF( .NOT.SELECT( KI ) ) 352 $ GO TO 80 353 END IF 354 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 355* 356 WORK( 1 ) = CMONE 357* 358* Form right-hand side. 359* 360 DO 40 K = 1, KI - 1 361 WORK( K ) = -T( K, KI ) 362 40 CONTINUE 363* 364* Solve the triangular system: 365* (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. 366* 367 DO 50 K = 1, KI - 1 368 T( K, K ) = T( K, K ) - T( KI, KI ) 369 IF( CABS1( T( K, K ) ).LT.SMIN ) 370 $ T( K, K ) = SMIN 371 50 CONTINUE 372* 373 IF( KI.GT.1 ) THEN 374 CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y', 375 $ KI-1, T, LDT, WORK( 1 ), SCALE, RWORK, 376 $ INFO ) 377 WORK( KI ) = SCALE 378 END IF 379* 380* Copy the vector x or Q*x to VR and normalize. 381* 382 IF( .NOT.OVER ) THEN 383 CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 ) 384* 385 II = ICAMAX( KI, VR( 1, IS ), 1 ) 386 REMAX = ONE / CABS1( VR( II, IS ) ) 387 CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 ) 388* 389 DO 60 K = KI + 1, N 390 VR( K, IS ) = CMZERO 391 60 CONTINUE 392 ELSE 393 IF( KI.GT.1 ) 394 $ CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ), 395 $ 1, CMPLX( SCALE ), VR( 1, KI ), 1 ) 396* 397 II = ICAMAX( N, VR( 1, KI ), 1 ) 398 REMAX = ONE / CABS1( VR( II, KI ) ) 399 CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 ) 400 END IF 401* 402* Set back the original diagonal elements of T. 403* 404 DO 70 K = 1, KI - 1 405 T( K, K ) = WORK( K+N ) 406 70 CONTINUE 407* 408 IS = IS - 1 409 80 CONTINUE 410 END IF 411* 412 IF( LEFTV ) THEN 413* 414* Compute left eigenvectors. 415* 416 IS = 1 417 DO 130 KI = 1, N 418* 419 IF( SOMEV ) THEN 420 IF( .NOT.SELECT( KI ) ) 421 $ GO TO 130 422 END IF 423 SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM ) 424* 425 WORK( N ) = CMONE 426* 427* Form right-hand side. 428* 429 DO 90 K = KI + 1, N 430 WORK( K ) = -CONJG( T( KI, K ) ) 431 90 CONTINUE 432* 433* Solve the triangular system: 434* (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK. 435* 436 DO 100 K = KI + 1, N 437 T( K, K ) = T( K, K ) - T( KI, KI ) 438 IF( CABS1( T( K, K ) ).LT.SMIN ) 439 $ T( K, K ) = SMIN 440 100 CONTINUE 441* 442 IF( KI.LT.N ) THEN 443 CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 444 $ 'Y', N-KI, T( KI+1, KI+1 ), LDT, 445 $ WORK( KI+1 ), SCALE, RWORK, INFO ) 446 WORK( KI ) = SCALE 447 END IF 448* 449* Copy the vector x or Q*x to VL and normalize. 450* 451 IF( .NOT.OVER ) THEN 452 CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 ) 453* 454 II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 455 REMAX = ONE / CABS1( VL( II, IS ) ) 456 CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 457* 458 DO 110 K = 1, KI - 1 459 VL( K, IS ) = CMZERO 460 110 CONTINUE 461 ELSE 462 IF( KI.LT.N ) 463 $ CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL, 464 $ WORK( KI+1 ), 1, CMPLX( SCALE ), 465 $ VL( 1, KI ), 1 ) 466* 467 II = ICAMAX( N, VL( 1, KI ), 1 ) 468 REMAX = ONE / CABS1( VL( II, KI ) ) 469 CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 ) 470 END IF 471* 472* Set back the original diagonal elements of T. 473* 474 DO 120 K = KI + 1, N 475 T( K, K ) = WORK( K+N ) 476 120 CONTINUE 477* 478 IS = IS + 1 479 130 CONTINUE 480 END IF 481* 482 RETURN 483* 484* End of CTREVC 485* 486 END 487c $Id$ 488