 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00002 $ LDVR, MM, M, WORK, 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 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00016 $ WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * STREVC computes some or all of the right and/or left eigenvectors of 00023 * a real upper quasi-triangular matrix T. 00024 * Matrices of this type are produced by the Schur factorization of 00025 * a real general matrix: A = Q*T*Q**T, as computed by SHSEQR. 00026 * 00027 * The right eigenvector x and the left eigenvector y of T corresponding 00028 * to an eigenvalue w are defined by: 00029 * 00030 * T*x = w*x, (y**T)*T = w*(y**T) 00031 * 00032 * where y**T denotes the transpose of y. 00033 * The eigenvalues are not input to this routine, but are read directly 00034 * from the diagonal blocks of T. 00035 * 00036 * This routine returns the matrices X and/or Y of right and left 00037 * eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 00038 * input matrix. If Q is the orthogonal factor that reduces a matrix 00039 * A to Schur form T, then Q*X and Q*Y are the matrices of right and 00040 * left eigenvectors of A. 00041 * 00042 * Arguments 00043 * ========= 00044 * 00045 * SIDE (input) CHARACTER*1 00046 * = 'R': compute right eigenvectors only; 00047 * = 'L': compute left eigenvectors only; 00048 * = 'B': compute both right and left eigenvectors. 00049 * 00050 * HOWMNY (input) CHARACTER*1 00051 * = 'A': compute all right and/or left eigenvectors; 00052 * = 'B': compute all right and/or left eigenvectors, 00053 * backtransformed by the matrices in VR and/or VL; 00054 * = 'S': compute selected right and/or left eigenvectors, 00055 * as indicated by the logical array SELECT. 00056 * 00057 * SELECT (input/output) LOGICAL array, dimension (N) 00058 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be 00059 * computed. 00060 * If w(j) is a real eigenvalue, the corresponding real 00061 * eigenvector is computed if SELECT(j) is .TRUE.. 00062 * If w(j) and w(j+1) are the real and imaginary parts of a 00063 * complex eigenvalue, the corresponding complex eigenvector is 00064 * computed if either SELECT(j) or SELECT(j+1) is .TRUE., and 00065 * on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to 00066 * .FALSE.. 00067 * Not referenced if HOWMNY = 'A' or 'B'. 00068 * 00069 * N (input) INTEGER 00070 * The order of the matrix T. N >= 0. 00071 * 00072 * T (input) REAL array, dimension (LDT,N) 00073 * The upper quasi-triangular matrix T in Schur canonical form. 00074 * 00075 * LDT (input) INTEGER 00076 * The leading dimension of the array T. LDT >= max(1,N). 00077 * 00078 * VL (input/output) REAL array, dimension (LDVL,MM) 00079 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00080 * contain an N-by-N matrix Q (usually the orthogonal matrix Q 00081 * of Schur vectors returned by SHSEQR). 00082 * On exit, if SIDE = 'L' or 'B', VL contains: 00083 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 00084 * if HOWMNY = 'B', the matrix Q*Y; 00085 * if HOWMNY = 'S', the left eigenvectors of T specified by 00086 * SELECT, stored consecutively in the columns 00087 * of VL, in the same order as their 00088 * eigenvalues. 00089 * A complex eigenvector corresponding to a complex eigenvalue 00090 * is stored in two consecutive columns, the first holding the 00091 * real part, and the second the imaginary part. 00092 * Not referenced if SIDE = 'R'. 00093 * 00094 * LDVL (input) INTEGER 00095 * The leading dimension of the array VL. LDVL >= 1, and if 00096 * SIDE = 'L' or 'B', LDVL >= N. 00097 * 00098 * VR (input/output) REAL array, dimension (LDVR,MM) 00099 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00100 * contain an N-by-N matrix Q (usually the orthogonal matrix Q 00101 * of Schur vectors returned by SHSEQR). 00102 * On exit, if SIDE = 'R' or 'B', VR contains: 00103 * if HOWMNY = 'A', the matrix X of right eigenvectors of T; 00104 * if HOWMNY = 'B', the matrix Q*X; 00105 * if HOWMNY = 'S', the right eigenvectors of T specified by 00106 * SELECT, stored consecutively in the columns 00107 * of VR, in the same order as their 00108 * eigenvalues. 00109 * A complex eigenvector corresponding to a complex eigenvalue 00110 * is stored in two consecutive columns, the first holding the 00111 * real part and the second the imaginary part. 00112 * Not referenced if SIDE = 'L'. 00113 * 00114 * LDVR (input) INTEGER 00115 * The leading dimension of the array VR. LDVR >= 1, and if 00116 * SIDE = 'R' or 'B', LDVR >= N. 00117 * 00118 * MM (input) INTEGER 00119 * The number of columns in the arrays VL and/or VR. MM >= M. 00120 * 00121 * M (output) INTEGER 00122 * The number of columns in the arrays VL and/or VR actually 00123 * used to store the eigenvectors. 00124 * If HOWMNY = 'A' or 'B', M is set to N. 00125 * Each selected real eigenvector occupies one column and each 00126 * selected complex eigenvector occupies two columns. 00127 * 00128 * WORK (workspace) REAL array, dimension (3*N) 00129 * 00130 * INFO (output) INTEGER 00131 * = 0: successful exit 00132 * < 0: if INFO = -i, the i-th argument had an illegal value 00133 * 00134 * Further Details 00135 * =============== 00136 * 00137 * The algorithm used in this program is basically backward (forward) 00138 * substitution, with scaling to make the the code robust against 00139 * possible overflow. 00140 * 00141 * Each eigenvector is normalized so that the element of largest 00142 * magnitude has magnitude 1; here the magnitude of a complex number 00143 * (x,y) is taken to be |x| + |y|. 00144 * 00145 * ===================================================================== 00146 * 00147 * .. Parameters .. 00148 REAL ZERO, ONE 00149 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00150 * .. 00151 * .. Local Scalars .. 00152 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV 00153 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2 00154 REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE, 00155 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR, 00156 $ XNORM 00157 * .. 00158 * .. External Functions .. 00159 LOGICAL LSAME 00160 INTEGER ISAMAX 00161 REAL SDOT, SLAMCH 00162 EXTERNAL LSAME, ISAMAX, SDOT, SLAMCH 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL, 00166 $ XERBLA 00167 * .. 00168 * .. Intrinsic Functions .. 00169 INTRINSIC ABS, MAX, SQRT 00170 * .. 00171 * .. Local Arrays .. 00172 REAL X( 2, 2 ) 00173 * .. 00174 * .. Executable Statements .. 00175 * 00176 * Decode and test the input parameters 00177 * 00178 BOTHV = LSAME( SIDE, 'B' ) 00179 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 00180 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 00181 * 00182 ALLV = LSAME( HOWMNY, 'A' ) 00183 OVER = LSAME( HOWMNY, 'B' ) 00184 SOMEV = LSAME( HOWMNY, 'S' ) 00185 * 00186 INFO = 0 00187 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 00188 INFO = -1 00189 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 00190 INFO = -2 00191 ELSE IF( N.LT.0 ) THEN 00192 INFO = -4 00193 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00194 INFO = -6 00195 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 00196 INFO = -8 00197 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 00198 INFO = -10 00199 ELSE 00200 * 00201 * Set M to the number of columns required to store the selected 00202 * eigenvectors, standardize the array SELECT if necessary, and 00203 * test MM. 00204 * 00205 IF( SOMEV ) THEN 00206 M = 0 00207 PAIR = .FALSE. 00208 DO 10 J = 1, N 00209 IF( PAIR ) THEN 00210 PAIR = .FALSE. 00211 SELECT( J ) = .FALSE. 00212 ELSE 00213 IF( J.LT.N ) THEN 00214 IF( T( J+1, J ).EQ.ZERO ) THEN 00215 IF( SELECT( J ) ) 00216 $ M = M + 1 00217 ELSE 00218 PAIR = .TRUE. 00219 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN 00220 SELECT( J ) = .TRUE. 00221 M = M + 2 00222 END IF 00223 END IF 00224 ELSE 00225 IF( SELECT( N ) ) 00226 $ M = M + 1 00227 END IF 00228 END IF 00229 10 CONTINUE 00230 ELSE 00231 M = N 00232 END IF 00233 * 00234 IF( MM.LT.M ) THEN 00235 INFO = -11 00236 END IF 00237 END IF 00238 IF( INFO.NE.0 ) THEN 00239 CALL XERBLA( 'STREVC', -INFO ) 00240 RETURN 00241 END IF 00242 * 00243 * Quick return if possible. 00244 * 00245 IF( N.EQ.0 ) 00246 $ RETURN 00247 * 00248 * Set the constants to control overflow. 00249 * 00250 UNFL = SLAMCH( 'Safe minimum' ) 00251 OVFL = ONE / UNFL 00252 CALL SLABAD( UNFL, OVFL ) 00253 ULP = SLAMCH( 'Precision' ) 00254 SMLNUM = UNFL*( N / ULP ) 00255 BIGNUM = ( ONE-ULP ) / SMLNUM 00256 * 00257 * Compute 1-norm of each column of strictly upper triangular 00258 * part of T to control overflow in triangular solver. 00259 * 00260 WORK( 1 ) = ZERO 00261 DO 30 J = 2, N 00262 WORK( J ) = ZERO 00263 DO 20 I = 1, J - 1 00264 WORK( J ) = WORK( J ) + ABS( T( I, J ) ) 00265 20 CONTINUE 00266 30 CONTINUE 00267 * 00268 * Index IP is used to specify the real or complex eigenvalue: 00269 * IP = 0, real eigenvalue, 00270 * 1, first of conjugate complex pair: (wr,wi) 00271 * -1, second of conjugate complex pair: (wr,wi) 00272 * 00273 N2 = 2*N 00274 * 00275 IF( RIGHTV ) THEN 00276 * 00277 * Compute right eigenvectors. 00278 * 00279 IP = 0 00280 IS = M 00281 DO 140 KI = N, 1, -1 00282 * 00283 IF( IP.EQ.1 ) 00284 $ GO TO 130 00285 IF( KI.EQ.1 ) 00286 $ GO TO 40 00287 IF( T( KI, KI-1 ).EQ.ZERO ) 00288 $ GO TO 40 00289 IP = -1 00290 * 00291 40 CONTINUE 00292 IF( SOMEV ) THEN 00293 IF( IP.EQ.0 ) THEN 00294 IF( .NOT.SELECT( KI ) ) 00295 $ GO TO 130 00296 ELSE 00297 IF( .NOT.SELECT( KI-1 ) ) 00298 $ GO TO 130 00299 END IF 00300 END IF 00301 * 00302 * Compute the KI-th eigenvalue (WR,WI). 00303 * 00304 WR = T( KI, KI ) 00305 WI = ZERO 00306 IF( IP.NE.0 ) 00307 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* 00308 $ SQRT( ABS( T( KI-1, KI ) ) ) 00309 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00310 * 00311 IF( IP.EQ.0 ) THEN 00312 * 00313 * Real right eigenvector 00314 * 00315 WORK( KI+N ) = ONE 00316 * 00317 * Form right-hand side 00318 * 00319 DO 50 K = 1, KI - 1 00320 WORK( K+N ) = -T( K, KI ) 00321 50 CONTINUE 00322 * 00323 * Solve the upper quasi-triangular system: 00324 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. 00325 * 00326 JNXT = KI - 1 00327 DO 60 J = KI - 1, 1, -1 00328 IF( J.GT.JNXT ) 00329 $ GO TO 60 00330 J1 = J 00331 J2 = J 00332 JNXT = J - 1 00333 IF( J.GT.1 ) THEN 00334 IF( T( J, J-1 ).NE.ZERO ) THEN 00335 J1 = J - 1 00336 JNXT = J - 2 00337 END IF 00338 END IF 00339 * 00340 IF( J1.EQ.J2 ) THEN 00341 * 00342 * 1-by-1 diagonal block 00343 * 00344 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00345 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00346 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00347 * 00348 * Scale X(1,1) to avoid overflow when updating 00349 * the right-hand side. 00350 * 00351 IF( XNORM.GT.ONE ) THEN 00352 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00353 X( 1, 1 ) = X( 1, 1 ) / XNORM 00354 SCALE = SCALE / XNORM 00355 END IF 00356 END IF 00357 * 00358 * Scale if necessary 00359 * 00360 IF( SCALE.NE.ONE ) 00361 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00362 WORK( J+N ) = X( 1, 1 ) 00363 * 00364 * Update right-hand side 00365 * 00366 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00367 $ WORK( 1+N ), 1 ) 00368 * 00369 ELSE 00370 * 00371 * 2-by-2 diagonal block 00372 * 00373 CALL SLALN2( .FALSE., 2, 1, SMIN, ONE, 00374 $ T( J-1, J-1 ), LDT, ONE, ONE, 00375 $ WORK( J-1+N ), N, WR, ZERO, X, 2, 00376 $ SCALE, XNORM, IERR ) 00377 * 00378 * Scale X(1,1) and X(2,1) to avoid overflow when 00379 * updating the right-hand side. 00380 * 00381 IF( XNORM.GT.ONE ) THEN 00382 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00383 IF( BETA.GT.BIGNUM / XNORM ) THEN 00384 X( 1, 1 ) = X( 1, 1 ) / XNORM 00385 X( 2, 1 ) = X( 2, 1 ) / XNORM 00386 SCALE = SCALE / XNORM 00387 END IF 00388 END IF 00389 * 00390 * Scale if necessary 00391 * 00392 IF( SCALE.NE.ONE ) 00393 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00394 WORK( J-1+N ) = X( 1, 1 ) 00395 WORK( J+N ) = X( 2, 1 ) 00396 * 00397 * Update right-hand side 00398 * 00399 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00400 $ WORK( 1+N ), 1 ) 00401 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00402 $ WORK( 1+N ), 1 ) 00403 END IF 00404 60 CONTINUE 00405 * 00406 * Copy the vector x or Q*x to VR and normalize. 00407 * 00408 IF( .NOT.OVER ) THEN 00409 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) 00410 * 00411 II = ISAMAX( KI, VR( 1, IS ), 1 ) 00412 REMAX = ONE / ABS( VR( II, IS ) ) 00413 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00414 * 00415 DO 70 K = KI + 1, N 00416 VR( K, IS ) = ZERO 00417 70 CONTINUE 00418 ELSE 00419 IF( KI.GT.1 ) 00420 $ CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR, 00421 $ WORK( 1+N ), 1, WORK( KI+N ), 00422 $ VR( 1, KI ), 1 ) 00423 * 00424 II = ISAMAX( N, VR( 1, KI ), 1 ) 00425 REMAX = ONE / ABS( VR( II, KI ) ) 00426 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 ) 00427 END IF 00428 * 00429 ELSE 00430 * 00431 * Complex right eigenvector. 00432 * 00433 * Initial solve 00434 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. 00435 * [ (T(KI,KI-1) T(KI,KI) ) ] 00436 * 00437 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN 00438 WORK( KI-1+N ) = ONE 00439 WORK( KI+N2 ) = WI / T( KI-1, KI ) 00440 ELSE 00441 WORK( KI-1+N ) = -WI / T( KI, KI-1 ) 00442 WORK( KI+N2 ) = ONE 00443 END IF 00444 WORK( KI+N ) = ZERO 00445 WORK( KI-1+N2 ) = ZERO 00446 * 00447 * Form right-hand side 00448 * 00449 DO 80 K = 1, KI - 2 00450 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) 00451 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) 00452 80 CONTINUE 00453 * 00454 * Solve upper quasi-triangular system: 00455 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) 00456 * 00457 JNXT = KI - 2 00458 DO 90 J = KI - 2, 1, -1 00459 IF( J.GT.JNXT ) 00460 $ GO TO 90 00461 J1 = J 00462 J2 = J 00463 JNXT = J - 1 00464 IF( J.GT.1 ) THEN 00465 IF( T( J, J-1 ).NE.ZERO ) THEN 00466 J1 = J - 1 00467 JNXT = J - 2 00468 END IF 00469 END IF 00470 * 00471 IF( J1.EQ.J2 ) THEN 00472 * 00473 * 1-by-1 diagonal block 00474 * 00475 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00476 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, 00477 $ X, 2, SCALE, XNORM, IERR ) 00478 * 00479 * Scale X(1,1) and X(1,2) to avoid overflow when 00480 * updating the right-hand side. 00481 * 00482 IF( XNORM.GT.ONE ) THEN 00483 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00484 X( 1, 1 ) = X( 1, 1 ) / XNORM 00485 X( 1, 2 ) = X( 1, 2 ) / XNORM 00486 SCALE = SCALE / XNORM 00487 END IF 00488 END IF 00489 * 00490 * Scale if necessary 00491 * 00492 IF( SCALE.NE.ONE ) THEN 00493 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00494 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00495 END IF 00496 WORK( J+N ) = X( 1, 1 ) 00497 WORK( J+N2 ) = X( 1, 2 ) 00498 * 00499 * Update the right-hand side 00500 * 00501 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00502 $ WORK( 1+N ), 1 ) 00503 CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, 00504 $ WORK( 1+N2 ), 1 ) 00505 * 00506 ELSE 00507 * 00508 * 2-by-2 diagonal block 00509 * 00510 CALL SLALN2( .FALSE., 2, 2, SMIN, ONE, 00511 $ T( J-1, J-1 ), LDT, ONE, ONE, 00512 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, 00513 $ XNORM, IERR ) 00514 * 00515 * Scale X to avoid overflow when updating 00516 * the right-hand side. 00517 * 00518 IF( XNORM.GT.ONE ) THEN 00519 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00520 IF( BETA.GT.BIGNUM / XNORM ) THEN 00521 REC = ONE / XNORM 00522 X( 1, 1 ) = X( 1, 1 )*REC 00523 X( 1, 2 ) = X( 1, 2 )*REC 00524 X( 2, 1 ) = X( 2, 1 )*REC 00525 X( 2, 2 ) = X( 2, 2 )*REC 00526 SCALE = SCALE*REC 00527 END IF 00528 END IF 00529 * 00530 * Scale if necessary 00531 * 00532 IF( SCALE.NE.ONE ) THEN 00533 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00534 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00535 END IF 00536 WORK( J-1+N ) = X( 1, 1 ) 00537 WORK( J+N ) = X( 2, 1 ) 00538 WORK( J-1+N2 ) = X( 1, 2 ) 00539 WORK( J+N2 ) = X( 2, 2 ) 00540 * 00541 * Update the right-hand side 00542 * 00543 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00544 $ WORK( 1+N ), 1 ) 00545 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00546 $ WORK( 1+N ), 1 ) 00547 CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, 00548 $ WORK( 1+N2 ), 1 ) 00549 CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, 00550 $ WORK( 1+N2 ), 1 ) 00551 END IF 00552 90 CONTINUE 00553 * 00554 * Copy the vector x or Q*x to VR and normalize. 00555 * 00556 IF( .NOT.OVER ) THEN 00557 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) 00558 CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) 00559 * 00560 EMAX = ZERO 00561 DO 100 K = 1, KI 00562 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ 00563 $ ABS( VR( K, IS ) ) ) 00564 100 CONTINUE 00565 * 00566 REMAX = ONE / EMAX 00567 CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) 00568 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00569 * 00570 DO 110 K = KI + 1, N 00571 VR( K, IS-1 ) = ZERO 00572 VR( K, IS ) = ZERO 00573 110 CONTINUE 00574 * 00575 ELSE 00576 * 00577 IF( KI.GT.2 ) THEN 00578 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00579 $ WORK( 1+N ), 1, WORK( KI-1+N ), 00580 $ VR( 1, KI-1 ), 1 ) 00581 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00582 $ WORK( 1+N2 ), 1, WORK( KI+N2 ), 00583 $ VR( 1, KI ), 1 ) 00584 ELSE 00585 CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) 00586 CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) 00587 END IF 00588 * 00589 EMAX = ZERO 00590 DO 120 K = 1, N 00591 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ 00592 $ ABS( VR( K, KI ) ) ) 00593 120 CONTINUE 00594 REMAX = ONE / EMAX 00595 CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) 00596 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 ) 00597 END IF 00598 END IF 00599 * 00600 IS = IS - 1 00601 IF( IP.NE.0 ) 00602 $ IS = IS - 1 00603 130 CONTINUE 00604 IF( IP.EQ.1 ) 00605 $ IP = 0 00606 IF( IP.EQ.-1 ) 00607 $ IP = 1 00608 140 CONTINUE 00609 END IF 00610 * 00611 IF( LEFTV ) THEN 00612 * 00613 * Compute left eigenvectors. 00614 * 00615 IP = 0 00616 IS = 1 00617 DO 260 KI = 1, N 00618 * 00619 IF( IP.EQ.-1 ) 00620 $ GO TO 250 00621 IF( KI.EQ.N ) 00622 $ GO TO 150 00623 IF( T( KI+1, KI ).EQ.ZERO ) 00624 $ GO TO 150 00625 IP = 1 00626 * 00627 150 CONTINUE 00628 IF( SOMEV ) THEN 00629 IF( .NOT.SELECT( KI ) ) 00630 $ GO TO 250 00631 END IF 00632 * 00633 * Compute the KI-th eigenvalue (WR,WI). 00634 * 00635 WR = T( KI, KI ) 00636 WI = ZERO 00637 IF( IP.NE.0 ) 00638 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* 00639 $ SQRT( ABS( T( KI+1, KI ) ) ) 00640 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00641 * 00642 IF( IP.EQ.0 ) THEN 00643 * 00644 * Real left eigenvector. 00645 * 00646 WORK( KI+N ) = ONE 00647 * 00648 * Form right-hand side 00649 * 00650 DO 160 K = KI + 1, N 00651 WORK( K+N ) = -T( KI, K ) 00652 160 CONTINUE 00653 * 00654 * Solve the quasi-triangular system: 00655 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK 00656 * 00657 VMAX = ONE 00658 VCRIT = BIGNUM 00659 * 00660 JNXT = KI + 1 00661 DO 170 J = KI + 1, N 00662 IF( J.LT.JNXT ) 00663 $ GO TO 170 00664 J1 = J 00665 J2 = J 00666 JNXT = J + 1 00667 IF( J.LT.N ) THEN 00668 IF( T( J+1, J ).NE.ZERO ) THEN 00669 J2 = J + 1 00670 JNXT = J + 2 00671 END IF 00672 END IF 00673 * 00674 IF( J1.EQ.J2 ) THEN 00675 * 00676 * 1-by-1 diagonal block 00677 * 00678 * Scale if necessary to avoid overflow when forming 00679 * the right-hand side. 00680 * 00681 IF( WORK( J ).GT.VCRIT ) THEN 00682 REC = ONE / VMAX 00683 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00684 VMAX = ONE 00685 VCRIT = BIGNUM 00686 END IF 00687 * 00688 WORK( J+N ) = WORK( J+N ) - 00689 $ SDOT( J-KI-1, T( KI+1, J ), 1, 00690 $ WORK( KI+1+N ), 1 ) 00691 * 00692 * Solve (T(J,J)-WR)**T*X = WORK 00693 * 00694 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00695 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00696 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00697 * 00698 * Scale if necessary 00699 * 00700 IF( SCALE.NE.ONE ) 00701 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00702 WORK( J+N ) = X( 1, 1 ) 00703 VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) 00704 VCRIT = BIGNUM / VMAX 00705 * 00706 ELSE 00707 * 00708 * 2-by-2 diagonal block 00709 * 00710 * Scale if necessary to avoid overflow when forming 00711 * the right-hand side. 00712 * 00713 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00714 IF( BETA.GT.VCRIT ) THEN 00715 REC = ONE / VMAX 00716 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00717 VMAX = ONE 00718 VCRIT = BIGNUM 00719 END IF 00720 * 00721 WORK( J+N ) = WORK( J+N ) - 00722 $ SDOT( J-KI-1, T( KI+1, J ), 1, 00723 $ WORK( KI+1+N ), 1 ) 00724 * 00725 WORK( J+1+N ) = WORK( J+1+N ) - 00726 $ SDOT( J-KI-1, T( KI+1, J+1 ), 1, 00727 $ WORK( KI+1+N ), 1 ) 00728 * 00729 * Solve 00730 * [T(J,J)-WR T(J,J+1) ]**T* X = SCALE*( WORK1 ) 00731 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) 00732 * 00733 CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), 00734 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00735 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00736 * 00737 * Scale if necessary 00738 * 00739 IF( SCALE.NE.ONE ) 00740 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00741 WORK( J+N ) = X( 1, 1 ) 00742 WORK( J+1+N ) = X( 2, 1 ) 00743 * 00744 VMAX = MAX( ABS( WORK( J+N ) ), 00745 $ ABS( WORK( J+1+N ) ), VMAX ) 00746 VCRIT = BIGNUM / VMAX 00747 * 00748 END IF 00749 170 CONTINUE 00750 * 00751 * Copy the vector x or Q*x to VL and normalize. 00752 * 00753 IF( .NOT.OVER ) THEN 00754 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 00755 * 00756 II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 00757 REMAX = ONE / ABS( VL( II, IS ) ) 00758 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00759 * 00760 DO 180 K = 1, KI - 1 00761 VL( K, IS ) = ZERO 00762 180 CONTINUE 00763 * 00764 ELSE 00765 * 00766 IF( KI.LT.N ) 00767 $ CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, 00768 $ WORK( KI+1+N ), 1, WORK( KI+N ), 00769 $ VL( 1, KI ), 1 ) 00770 * 00771 II = ISAMAX( N, VL( 1, KI ), 1 ) 00772 REMAX = ONE / ABS( VL( II, KI ) ) 00773 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 ) 00774 * 00775 END IF 00776 * 00777 ELSE 00778 * 00779 * Complex left eigenvector. 00780 * 00781 * Initial solve: 00782 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0. 00783 * ((T(KI+1,KI) T(KI+1,KI+1)) ) 00784 * 00785 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN 00786 WORK( KI+N ) = WI / T( KI, KI+1 ) 00787 WORK( KI+1+N2 ) = ONE 00788 ELSE 00789 WORK( KI+N ) = ONE 00790 WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) 00791 END IF 00792 WORK( KI+1+N ) = ZERO 00793 WORK( KI+N2 ) = ZERO 00794 * 00795 * Form right-hand side 00796 * 00797 DO 190 K = KI + 2, N 00798 WORK( K+N ) = -WORK( KI+N )*T( KI, K ) 00799 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) 00800 190 CONTINUE 00801 * 00802 * Solve complex quasi-triangular system: 00803 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 00804 * 00805 VMAX = ONE 00806 VCRIT = BIGNUM 00807 * 00808 JNXT = KI + 2 00809 DO 200 J = KI + 2, N 00810 IF( J.LT.JNXT ) 00811 $ GO TO 200 00812 J1 = J 00813 J2 = J 00814 JNXT = J + 1 00815 IF( J.LT.N ) THEN 00816 IF( T( J+1, J ).NE.ZERO ) THEN 00817 J2 = J + 1 00818 JNXT = J + 2 00819 END IF 00820 END IF 00821 * 00822 IF( J1.EQ.J2 ) THEN 00823 * 00824 * 1-by-1 diagonal block 00825 * 00826 * Scale if necessary to avoid overflow when 00827 * forming the right-hand side elements. 00828 * 00829 IF( WORK( J ).GT.VCRIT ) THEN 00830 REC = ONE / VMAX 00831 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00832 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00833 VMAX = ONE 00834 VCRIT = BIGNUM 00835 END IF 00836 * 00837 WORK( J+N ) = WORK( J+N ) - 00838 $ SDOT( J-KI-2, T( KI+2, J ), 1, 00839 $ WORK( KI+2+N ), 1 ) 00840 WORK( J+N2 ) = WORK( J+N2 ) - 00841 $ SDOT( J-KI-2, T( KI+2, J ), 1, 00842 $ WORK( KI+2+N2 ), 1 ) 00843 * 00844 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 00845 * 00846 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00847 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00848 $ -WI, X, 2, SCALE, XNORM, IERR ) 00849 * 00850 * Scale if necessary 00851 * 00852 IF( SCALE.NE.ONE ) THEN 00853 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00854 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 00855 END IF 00856 WORK( J+N ) = X( 1, 1 ) 00857 WORK( J+N2 ) = X( 1, 2 ) 00858 VMAX = MAX( ABS( WORK( J+N ) ), 00859 $ ABS( WORK( J+N2 ) ), VMAX ) 00860 VCRIT = BIGNUM / VMAX 00861 * 00862 ELSE 00863 * 00864 * 2-by-2 diagonal block 00865 * 00866 * Scale if necessary to avoid overflow when forming 00867 * the right-hand side elements. 00868 * 00869 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00870 IF( BETA.GT.VCRIT ) THEN 00871 REC = ONE / VMAX 00872 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00873 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00874 VMAX = ONE 00875 VCRIT = BIGNUM 00876 END IF 00877 * 00878 WORK( J+N ) = WORK( J+N ) - 00879 $ SDOT( J-KI-2, T( KI+2, J ), 1, 00880 $ WORK( KI+2+N ), 1 ) 00881 * 00882 WORK( J+N2 ) = WORK( J+N2 ) - 00883 $ SDOT( J-KI-2, T( KI+2, J ), 1, 00884 $ WORK( KI+2+N2 ), 1 ) 00885 * 00886 WORK( J+1+N ) = WORK( J+1+N ) - 00887 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1, 00888 $ WORK( KI+2+N ), 1 ) 00889 * 00890 WORK( J+1+N2 ) = WORK( J+1+N2 ) - 00891 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1, 00892 $ WORK( KI+2+N2 ), 1 ) 00893 * 00894 * Solve 2-by-2 complex linear equation 00895 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B 00896 * ([T(j+1,j) T(j+1,j+1)] ) 00897 * 00898 CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), 00899 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00900 $ -WI, X, 2, SCALE, XNORM, IERR ) 00901 * 00902 * Scale if necessary 00903 * 00904 IF( SCALE.NE.ONE ) THEN 00905 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00906 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 00907 END IF 00908 WORK( J+N ) = X( 1, 1 ) 00909 WORK( J+N2 ) = X( 1, 2 ) 00910 WORK( J+1+N ) = X( 2, 1 ) 00911 WORK( J+1+N2 ) = X( 2, 2 ) 00912 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), 00913 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) 00914 VCRIT = BIGNUM / VMAX 00915 * 00916 END IF 00917 200 CONTINUE 00918 * 00919 * Copy the vector x or Q*x to VL and normalize. 00920 * 00921 IF( .NOT.OVER ) THEN 00922 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 00923 CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), 00924 $ 1 ) 00925 * 00926 EMAX = ZERO 00927 DO 220 K = KI, N 00928 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ 00929 $ ABS( VL( K, IS+1 ) ) ) 00930 220 CONTINUE 00931 REMAX = ONE / EMAX 00932 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00933 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) 00934 * 00935 DO 230 K = 1, KI - 1 00936 VL( K, IS ) = ZERO 00937 VL( K, IS+1 ) = ZERO 00938 230 CONTINUE 00939 ELSE 00940 IF( KI.LT.N-1 ) THEN 00941 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 00942 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), 00943 $ VL( 1, KI ), 1 ) 00944 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 00945 $ LDVL, WORK( KI+2+N2 ), 1, 00946 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 00947 ELSE 00948 CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) 00949 CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 00950 END IF 00951 * 00952 EMAX = ZERO 00953 DO 240 K = 1, N 00954 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ 00955 $ ABS( VL( K, KI+1 ) ) ) 00956 240 CONTINUE 00957 REMAX = ONE / EMAX 00958 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 ) 00959 CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) 00960 * 00961 END IF 00962 * 00963 END IF 00964 * 00965 IS = IS + 1 00966 IF( IP.NE.0 ) 00967 $ IS = IS + 1 00968 250 CONTINUE 00969 IF( IP.EQ.-1 ) 00970 $ IP = 0 00971 IF( IP.EQ.1 ) 00972 $ IP = -1 00973 * 00974 260 CONTINUE 00975 * 00976 END IF 00977 * 00978 RETURN 00979 * 00980 * End of STREVC 00981 * 00982 END
1.7.3 
