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