![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTREVC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTREVC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTREVC( 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 * DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00031 * $ WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DTREVC 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 DHSEQR. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DHSEQR). 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 DOUBLE PRECISION 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 DHSEQR). 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 DOUBLE PRECISION 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 doubleOTHERcomputational 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 DTREVC( 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 DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00237 $ WORK( * ) 00238 * .. 00239 * 00240 * ===================================================================== 00241 * 00242 * .. Parameters .. 00243 DOUBLE PRECISION ZERO, ONE 00244 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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 DOUBLE PRECISION 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 IDAMAX 00256 DOUBLE PRECISION DDOT, DLAMCH 00257 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 00258 * .. 00259 * .. External Subroutines .. 00260 EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA 00261 * .. 00262 * .. Intrinsic Functions .. 00263 INTRINSIC ABS, MAX, SQRT 00264 * .. 00265 * .. Local Arrays .. 00266 DOUBLE PRECISION X( 2, 2 ) 00267 * .. 00268 * .. Executable Statements .. 00269 * 00270 * Decode and test the input parameters 00271 * 00272 BOTHV = LSAME( SIDE, 'B' ) 00273 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 00274 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 00275 * 00276 ALLV = LSAME( HOWMNY, 'A' ) 00277 OVER = LSAME( HOWMNY, 'B' ) 00278 SOMEV = LSAME( HOWMNY, 'S' ) 00279 * 00280 INFO = 0 00281 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 00282 INFO = -1 00283 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 00284 INFO = -2 00285 ELSE IF( N.LT.0 ) THEN 00286 INFO = -4 00287 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00288 INFO = -6 00289 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 00290 INFO = -8 00291 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 00292 INFO = -10 00293 ELSE 00294 * 00295 * Set M to the number of columns required to store the selected 00296 * eigenvectors, standardize the array SELECT if necessary, and 00297 * test MM. 00298 * 00299 IF( SOMEV ) THEN 00300 M = 0 00301 PAIR = .FALSE. 00302 DO 10 J = 1, N 00303 IF( PAIR ) THEN 00304 PAIR = .FALSE. 00305 SELECT( J ) = .FALSE. 00306 ELSE 00307 IF( J.LT.N ) THEN 00308 IF( T( J+1, J ).EQ.ZERO ) THEN 00309 IF( SELECT( J ) ) 00310 $ M = M + 1 00311 ELSE 00312 PAIR = .TRUE. 00313 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN 00314 SELECT( J ) = .TRUE. 00315 M = M + 2 00316 END IF 00317 END IF 00318 ELSE 00319 IF( SELECT( N ) ) 00320 $ M = M + 1 00321 END IF 00322 END IF 00323 10 CONTINUE 00324 ELSE 00325 M = N 00326 END IF 00327 * 00328 IF( MM.LT.M ) THEN 00329 INFO = -11 00330 END IF 00331 END IF 00332 IF( INFO.NE.0 ) THEN 00333 CALL XERBLA( 'DTREVC', -INFO ) 00334 RETURN 00335 END IF 00336 * 00337 * Quick return if possible. 00338 * 00339 IF( N.EQ.0 ) 00340 $ RETURN 00341 * 00342 * Set the constants to control overflow. 00343 * 00344 UNFL = DLAMCH( 'Safe minimum' ) 00345 OVFL = ONE / UNFL 00346 CALL DLABAD( UNFL, OVFL ) 00347 ULP = DLAMCH( 'Precision' ) 00348 SMLNUM = UNFL*( N / ULP ) 00349 BIGNUM = ( ONE-ULP ) / SMLNUM 00350 * 00351 * Compute 1-norm of each column of strictly upper triangular 00352 * part of T to control overflow in triangular solver. 00353 * 00354 WORK( 1 ) = ZERO 00355 DO 30 J = 2, N 00356 WORK( J ) = ZERO 00357 DO 20 I = 1, J - 1 00358 WORK( J ) = WORK( J ) + ABS( T( I, J ) ) 00359 20 CONTINUE 00360 30 CONTINUE 00361 * 00362 * Index IP is used to specify the real or complex eigenvalue: 00363 * IP = 0, real eigenvalue, 00364 * 1, first of conjugate complex pair: (wr,wi) 00365 * -1, second of conjugate complex pair: (wr,wi) 00366 * 00367 N2 = 2*N 00368 * 00369 IF( RIGHTV ) THEN 00370 * 00371 * Compute right eigenvectors. 00372 * 00373 IP = 0 00374 IS = M 00375 DO 140 KI = N, 1, -1 00376 * 00377 IF( IP.EQ.1 ) 00378 $ GO TO 130 00379 IF( KI.EQ.1 ) 00380 $ GO TO 40 00381 IF( T( KI, KI-1 ).EQ.ZERO ) 00382 $ GO TO 40 00383 IP = -1 00384 * 00385 40 CONTINUE 00386 IF( SOMEV ) THEN 00387 IF( IP.EQ.0 ) THEN 00388 IF( .NOT.SELECT( KI ) ) 00389 $ GO TO 130 00390 ELSE 00391 IF( .NOT.SELECT( KI-1 ) ) 00392 $ GO TO 130 00393 END IF 00394 END IF 00395 * 00396 * Compute the KI-th eigenvalue (WR,WI). 00397 * 00398 WR = T( KI, KI ) 00399 WI = ZERO 00400 IF( IP.NE.0 ) 00401 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* 00402 $ SQRT( ABS( T( KI-1, KI ) ) ) 00403 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00404 * 00405 IF( IP.EQ.0 ) THEN 00406 * 00407 * Real right eigenvector 00408 * 00409 WORK( KI+N ) = ONE 00410 * 00411 * Form right-hand side 00412 * 00413 DO 50 K = 1, KI - 1 00414 WORK( K+N ) = -T( K, KI ) 00415 50 CONTINUE 00416 * 00417 * Solve the upper quasi-triangular system: 00418 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. 00419 * 00420 JNXT = KI - 1 00421 DO 60 J = KI - 1, 1, -1 00422 IF( J.GT.JNXT ) 00423 $ GO TO 60 00424 J1 = J 00425 J2 = J 00426 JNXT = J - 1 00427 IF( J.GT.1 ) THEN 00428 IF( T( J, J-1 ).NE.ZERO ) THEN 00429 J1 = J - 1 00430 JNXT = J - 2 00431 END IF 00432 END IF 00433 * 00434 IF( J1.EQ.J2 ) THEN 00435 * 00436 * 1-by-1 diagonal block 00437 * 00438 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00439 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00440 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00441 * 00442 * Scale X(1,1) to avoid overflow when updating 00443 * the right-hand side. 00444 * 00445 IF( XNORM.GT.ONE ) THEN 00446 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00447 X( 1, 1 ) = X( 1, 1 ) / XNORM 00448 SCALE = SCALE / XNORM 00449 END IF 00450 END IF 00451 * 00452 * Scale if necessary 00453 * 00454 IF( SCALE.NE.ONE ) 00455 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00456 WORK( J+N ) = X( 1, 1 ) 00457 * 00458 * Update right-hand side 00459 * 00460 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00461 $ WORK( 1+N ), 1 ) 00462 * 00463 ELSE 00464 * 00465 * 2-by-2 diagonal block 00466 * 00467 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, 00468 $ T( J-1, J-1 ), LDT, ONE, ONE, 00469 $ WORK( J-1+N ), N, WR, ZERO, X, 2, 00470 $ SCALE, XNORM, IERR ) 00471 * 00472 * Scale X(1,1) and X(2,1) to avoid overflow when 00473 * updating the right-hand side. 00474 * 00475 IF( XNORM.GT.ONE ) THEN 00476 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00477 IF( BETA.GT.BIGNUM / XNORM ) THEN 00478 X( 1, 1 ) = X( 1, 1 ) / XNORM 00479 X( 2, 1 ) = X( 2, 1 ) / XNORM 00480 SCALE = SCALE / XNORM 00481 END IF 00482 END IF 00483 * 00484 * Scale if necessary 00485 * 00486 IF( SCALE.NE.ONE ) 00487 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00488 WORK( J-1+N ) = X( 1, 1 ) 00489 WORK( J+N ) = X( 2, 1 ) 00490 * 00491 * Update right-hand side 00492 * 00493 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00494 $ WORK( 1+N ), 1 ) 00495 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00496 $ WORK( 1+N ), 1 ) 00497 END IF 00498 60 CONTINUE 00499 * 00500 * Copy the vector x or Q*x to VR and normalize. 00501 * 00502 IF( .NOT.OVER ) THEN 00503 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) 00504 * 00505 II = IDAMAX( KI, VR( 1, IS ), 1 ) 00506 REMAX = ONE / ABS( VR( II, IS ) ) 00507 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00508 * 00509 DO 70 K = KI + 1, N 00510 VR( K, IS ) = ZERO 00511 70 CONTINUE 00512 ELSE 00513 IF( KI.GT.1 ) 00514 $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR, 00515 $ WORK( 1+N ), 1, WORK( KI+N ), 00516 $ VR( 1, KI ), 1 ) 00517 * 00518 II = IDAMAX( N, VR( 1, KI ), 1 ) 00519 REMAX = ONE / ABS( VR( II, KI ) ) 00520 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 00521 END IF 00522 * 00523 ELSE 00524 * 00525 * Complex right eigenvector. 00526 * 00527 * Initial solve 00528 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. 00529 * [ (T(KI,KI-1) T(KI,KI) ) ] 00530 * 00531 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN 00532 WORK( KI-1+N ) = ONE 00533 WORK( KI+N2 ) = WI / T( KI-1, KI ) 00534 ELSE 00535 WORK( KI-1+N ) = -WI / T( KI, KI-1 ) 00536 WORK( KI+N2 ) = ONE 00537 END IF 00538 WORK( KI+N ) = ZERO 00539 WORK( KI-1+N2 ) = ZERO 00540 * 00541 * Form right-hand side 00542 * 00543 DO 80 K = 1, KI - 2 00544 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) 00545 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) 00546 80 CONTINUE 00547 * 00548 * Solve upper quasi-triangular system: 00549 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) 00550 * 00551 JNXT = KI - 2 00552 DO 90 J = KI - 2, 1, -1 00553 IF( J.GT.JNXT ) 00554 $ GO TO 90 00555 J1 = J 00556 J2 = J 00557 JNXT = J - 1 00558 IF( J.GT.1 ) THEN 00559 IF( T( J, J-1 ).NE.ZERO ) THEN 00560 J1 = J - 1 00561 JNXT = J - 2 00562 END IF 00563 END IF 00564 * 00565 IF( J1.EQ.J2 ) THEN 00566 * 00567 * 1-by-1 diagonal block 00568 * 00569 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00570 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, 00571 $ X, 2, SCALE, XNORM, IERR ) 00572 * 00573 * Scale X(1,1) and X(1,2) to avoid overflow when 00574 * updating the right-hand side. 00575 * 00576 IF( XNORM.GT.ONE ) THEN 00577 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00578 X( 1, 1 ) = X( 1, 1 ) / XNORM 00579 X( 1, 2 ) = X( 1, 2 ) / XNORM 00580 SCALE = SCALE / XNORM 00581 END IF 00582 END IF 00583 * 00584 * Scale if necessary 00585 * 00586 IF( SCALE.NE.ONE ) THEN 00587 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00588 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00589 END IF 00590 WORK( J+N ) = X( 1, 1 ) 00591 WORK( J+N2 ) = X( 1, 2 ) 00592 * 00593 * Update the right-hand side 00594 * 00595 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00596 $ WORK( 1+N ), 1 ) 00597 CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, 00598 $ WORK( 1+N2 ), 1 ) 00599 * 00600 ELSE 00601 * 00602 * 2-by-2 diagonal block 00603 * 00604 CALL DLALN2( .FALSE., 2, 2, SMIN, ONE, 00605 $ T( J-1, J-1 ), LDT, ONE, ONE, 00606 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, 00607 $ XNORM, IERR ) 00608 * 00609 * Scale X to avoid overflow when updating 00610 * the right-hand side. 00611 * 00612 IF( XNORM.GT.ONE ) THEN 00613 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00614 IF( BETA.GT.BIGNUM / XNORM ) THEN 00615 REC = ONE / XNORM 00616 X( 1, 1 ) = X( 1, 1 )*REC 00617 X( 1, 2 ) = X( 1, 2 )*REC 00618 X( 2, 1 ) = X( 2, 1 )*REC 00619 X( 2, 2 ) = X( 2, 2 )*REC 00620 SCALE = SCALE*REC 00621 END IF 00622 END IF 00623 * 00624 * Scale if necessary 00625 * 00626 IF( SCALE.NE.ONE ) THEN 00627 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00628 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00629 END IF 00630 WORK( J-1+N ) = X( 1, 1 ) 00631 WORK( J+N ) = X( 2, 1 ) 00632 WORK( J-1+N2 ) = X( 1, 2 ) 00633 WORK( J+N2 ) = X( 2, 2 ) 00634 * 00635 * Update the right-hand side 00636 * 00637 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00638 $ WORK( 1+N ), 1 ) 00639 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00640 $ WORK( 1+N ), 1 ) 00641 CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, 00642 $ WORK( 1+N2 ), 1 ) 00643 CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, 00644 $ WORK( 1+N2 ), 1 ) 00645 END IF 00646 90 CONTINUE 00647 * 00648 * Copy the vector x or Q*x to VR and normalize. 00649 * 00650 IF( .NOT.OVER ) THEN 00651 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) 00652 CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) 00653 * 00654 EMAX = ZERO 00655 DO 100 K = 1, KI 00656 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ 00657 $ ABS( VR( K, IS ) ) ) 00658 100 CONTINUE 00659 * 00660 REMAX = ONE / EMAX 00661 CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) 00662 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00663 * 00664 DO 110 K = KI + 1, N 00665 VR( K, IS-1 ) = ZERO 00666 VR( K, IS ) = ZERO 00667 110 CONTINUE 00668 * 00669 ELSE 00670 * 00671 IF( KI.GT.2 ) THEN 00672 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00673 $ WORK( 1+N ), 1, WORK( KI-1+N ), 00674 $ VR( 1, KI-1 ), 1 ) 00675 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00676 $ WORK( 1+N2 ), 1, WORK( KI+N2 ), 00677 $ VR( 1, KI ), 1 ) 00678 ELSE 00679 CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) 00680 CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) 00681 END IF 00682 * 00683 EMAX = ZERO 00684 DO 120 K = 1, N 00685 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ 00686 $ ABS( VR( K, KI ) ) ) 00687 120 CONTINUE 00688 REMAX = ONE / EMAX 00689 CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) 00690 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 00691 END IF 00692 END IF 00693 * 00694 IS = IS - 1 00695 IF( IP.NE.0 ) 00696 $ IS = IS - 1 00697 130 CONTINUE 00698 IF( IP.EQ.1 ) 00699 $ IP = 0 00700 IF( IP.EQ.-1 ) 00701 $ IP = 1 00702 140 CONTINUE 00703 END IF 00704 * 00705 IF( LEFTV ) THEN 00706 * 00707 * Compute left eigenvectors. 00708 * 00709 IP = 0 00710 IS = 1 00711 DO 260 KI = 1, N 00712 * 00713 IF( IP.EQ.-1 ) 00714 $ GO TO 250 00715 IF( KI.EQ.N ) 00716 $ GO TO 150 00717 IF( T( KI+1, KI ).EQ.ZERO ) 00718 $ GO TO 150 00719 IP = 1 00720 * 00721 150 CONTINUE 00722 IF( SOMEV ) THEN 00723 IF( .NOT.SELECT( KI ) ) 00724 $ GO TO 250 00725 END IF 00726 * 00727 * Compute the KI-th eigenvalue (WR,WI). 00728 * 00729 WR = T( KI, KI ) 00730 WI = ZERO 00731 IF( IP.NE.0 ) 00732 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* 00733 $ SQRT( ABS( T( KI+1, KI ) ) ) 00734 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00735 * 00736 IF( IP.EQ.0 ) THEN 00737 * 00738 * Real left eigenvector. 00739 * 00740 WORK( KI+N ) = ONE 00741 * 00742 * Form right-hand side 00743 * 00744 DO 160 K = KI + 1, N 00745 WORK( K+N ) = -T( KI, K ) 00746 160 CONTINUE 00747 * 00748 * Solve the quasi-triangular system: 00749 * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK 00750 * 00751 VMAX = ONE 00752 VCRIT = BIGNUM 00753 * 00754 JNXT = KI + 1 00755 DO 170 J = KI + 1, N 00756 IF( J.LT.JNXT ) 00757 $ GO TO 170 00758 J1 = J 00759 J2 = J 00760 JNXT = J + 1 00761 IF( J.LT.N ) THEN 00762 IF( T( J+1, J ).NE.ZERO ) THEN 00763 J2 = J + 1 00764 JNXT = J + 2 00765 END IF 00766 END IF 00767 * 00768 IF( J1.EQ.J2 ) THEN 00769 * 00770 * 1-by-1 diagonal block 00771 * 00772 * Scale if necessary to avoid overflow when forming 00773 * the right-hand side. 00774 * 00775 IF( WORK( J ).GT.VCRIT ) THEN 00776 REC = ONE / VMAX 00777 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00778 VMAX = ONE 00779 VCRIT = BIGNUM 00780 END IF 00781 * 00782 WORK( J+N ) = WORK( J+N ) - 00783 $ DDOT( J-KI-1, T( KI+1, J ), 1, 00784 $ WORK( KI+1+N ), 1 ) 00785 * 00786 * Solve (T(J,J)-WR)**T*X = WORK 00787 * 00788 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00789 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00790 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00791 * 00792 * Scale if necessary 00793 * 00794 IF( SCALE.NE.ONE ) 00795 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00796 WORK( J+N ) = X( 1, 1 ) 00797 VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) 00798 VCRIT = BIGNUM / VMAX 00799 * 00800 ELSE 00801 * 00802 * 2-by-2 diagonal block 00803 * 00804 * Scale if necessary to avoid overflow when forming 00805 * the right-hand side. 00806 * 00807 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00808 IF( BETA.GT.VCRIT ) THEN 00809 REC = ONE / VMAX 00810 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00811 VMAX = ONE 00812 VCRIT = BIGNUM 00813 END IF 00814 * 00815 WORK( J+N ) = WORK( J+N ) - 00816 $ DDOT( J-KI-1, T( KI+1, J ), 1, 00817 $ WORK( KI+1+N ), 1 ) 00818 * 00819 WORK( J+1+N ) = WORK( J+1+N ) - 00820 $ DDOT( J-KI-1, T( KI+1, J+1 ), 1, 00821 $ WORK( KI+1+N ), 1 ) 00822 * 00823 * Solve 00824 * [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 ) 00825 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) 00826 * 00827 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), 00828 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00829 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00830 * 00831 * Scale if necessary 00832 * 00833 IF( SCALE.NE.ONE ) 00834 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00835 WORK( J+N ) = X( 1, 1 ) 00836 WORK( J+1+N ) = X( 2, 1 ) 00837 * 00838 VMAX = MAX( ABS( WORK( J+N ) ), 00839 $ ABS( WORK( J+1+N ) ), VMAX ) 00840 VCRIT = BIGNUM / VMAX 00841 * 00842 END IF 00843 170 CONTINUE 00844 * 00845 * Copy the vector x or Q*x to VL and normalize. 00846 * 00847 IF( .NOT.OVER ) THEN 00848 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 00849 * 00850 II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 00851 REMAX = ONE / ABS( VL( II, IS ) ) 00852 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00853 * 00854 DO 180 K = 1, KI - 1 00855 VL( K, IS ) = ZERO 00856 180 CONTINUE 00857 * 00858 ELSE 00859 * 00860 IF( KI.LT.N ) 00861 $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, 00862 $ WORK( KI+1+N ), 1, WORK( KI+N ), 00863 $ VL( 1, KI ), 1 ) 00864 * 00865 II = IDAMAX( N, VL( 1, KI ), 1 ) 00866 REMAX = ONE / ABS( VL( II, KI ) ) 00867 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 00868 * 00869 END IF 00870 * 00871 ELSE 00872 * 00873 * Complex left eigenvector. 00874 * 00875 * Initial solve: 00876 * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0. 00877 * ((T(KI+1,KI) T(KI+1,KI+1)) ) 00878 * 00879 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN 00880 WORK( KI+N ) = WI / T( KI, KI+1 ) 00881 WORK( KI+1+N2 ) = ONE 00882 ELSE 00883 WORK( KI+N ) = ONE 00884 WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) 00885 END IF 00886 WORK( KI+1+N ) = ZERO 00887 WORK( KI+N2 ) = ZERO 00888 * 00889 * Form right-hand side 00890 * 00891 DO 190 K = KI + 2, N 00892 WORK( K+N ) = -WORK( KI+N )*T( KI, K ) 00893 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) 00894 190 CONTINUE 00895 * 00896 * Solve complex quasi-triangular system: 00897 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 00898 * 00899 VMAX = ONE 00900 VCRIT = BIGNUM 00901 * 00902 JNXT = KI + 2 00903 DO 200 J = KI + 2, N 00904 IF( J.LT.JNXT ) 00905 $ GO TO 200 00906 J1 = J 00907 J2 = J 00908 JNXT = J + 1 00909 IF( J.LT.N ) THEN 00910 IF( T( J+1, J ).NE.ZERO ) THEN 00911 J2 = J + 1 00912 JNXT = J + 2 00913 END IF 00914 END IF 00915 * 00916 IF( J1.EQ.J2 ) THEN 00917 * 00918 * 1-by-1 diagonal block 00919 * 00920 * Scale if necessary to avoid overflow when 00921 * forming the right-hand side elements. 00922 * 00923 IF( WORK( J ).GT.VCRIT ) THEN 00924 REC = ONE / VMAX 00925 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00926 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00927 VMAX = ONE 00928 VCRIT = BIGNUM 00929 END IF 00930 * 00931 WORK( J+N ) = WORK( J+N ) - 00932 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00933 $ WORK( KI+2+N ), 1 ) 00934 WORK( J+N2 ) = WORK( J+N2 ) - 00935 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00936 $ WORK( KI+2+N2 ), 1 ) 00937 * 00938 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 00939 * 00940 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00941 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00942 $ -WI, X, 2, SCALE, XNORM, IERR ) 00943 * 00944 * Scale if necessary 00945 * 00946 IF( SCALE.NE.ONE ) THEN 00947 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00948 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 00949 END IF 00950 WORK( J+N ) = X( 1, 1 ) 00951 WORK( J+N2 ) = X( 1, 2 ) 00952 VMAX = MAX( ABS( WORK( J+N ) ), 00953 $ ABS( WORK( J+N2 ) ), VMAX ) 00954 VCRIT = BIGNUM / VMAX 00955 * 00956 ELSE 00957 * 00958 * 2-by-2 diagonal block 00959 * 00960 * Scale if necessary to avoid overflow when forming 00961 * the right-hand side elements. 00962 * 00963 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00964 IF( BETA.GT.VCRIT ) THEN 00965 REC = ONE / VMAX 00966 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00967 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00968 VMAX = ONE 00969 VCRIT = BIGNUM 00970 END IF 00971 * 00972 WORK( J+N ) = WORK( J+N ) - 00973 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00974 $ WORK( KI+2+N ), 1 ) 00975 * 00976 WORK( J+N2 ) = WORK( J+N2 ) - 00977 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00978 $ WORK( KI+2+N2 ), 1 ) 00979 * 00980 WORK( J+1+N ) = WORK( J+1+N ) - 00981 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 00982 $ WORK( KI+2+N ), 1 ) 00983 * 00984 WORK( J+1+N2 ) = WORK( J+1+N2 ) - 00985 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 00986 $ WORK( KI+2+N2 ), 1 ) 00987 * 00988 * Solve 2-by-2 complex linear equation 00989 * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B 00990 * ([T(j+1,j) T(j+1,j+1)] ) 00991 * 00992 CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), 00993 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00994 $ -WI, X, 2, SCALE, XNORM, IERR ) 00995 * 00996 * Scale if necessary 00997 * 00998 IF( SCALE.NE.ONE ) THEN 00999 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 01000 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 01001 END IF 01002 WORK( J+N ) = X( 1, 1 ) 01003 WORK( J+N2 ) = X( 1, 2 ) 01004 WORK( J+1+N ) = X( 2, 1 ) 01005 WORK( J+1+N2 ) = X( 2, 2 ) 01006 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), 01007 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) 01008 VCRIT = BIGNUM / VMAX 01009 * 01010 END IF 01011 200 CONTINUE 01012 * 01013 * Copy the vector x or Q*x to VL and normalize. 01014 * 01015 IF( .NOT.OVER ) THEN 01016 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 01017 CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), 01018 $ 1 ) 01019 * 01020 EMAX = ZERO 01021 DO 220 K = KI, N 01022 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ 01023 $ ABS( VL( K, IS+1 ) ) ) 01024 220 CONTINUE 01025 REMAX = ONE / EMAX 01026 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 01027 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) 01028 * 01029 DO 230 K = 1, KI - 1 01030 VL( K, IS ) = ZERO 01031 VL( K, IS+1 ) = ZERO 01032 230 CONTINUE 01033 ELSE 01034 IF( KI.LT.N-1 ) THEN 01035 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 01036 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), 01037 $ VL( 1, KI ), 1 ) 01038 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 01039 $ LDVL, WORK( KI+2+N2 ), 1, 01040 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 01041 ELSE 01042 CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) 01043 CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 01044 END IF 01045 * 01046 EMAX = ZERO 01047 DO 240 K = 1, N 01048 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ 01049 $ ABS( VL( K, KI+1 ) ) ) 01050 240 CONTINUE 01051 REMAX = ONE / EMAX 01052 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 01053 CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) 01054 * 01055 END IF 01056 * 01057 END IF 01058 * 01059 IS = IS + 1 01060 IF( IP.NE.0 ) 01061 $ IS = IS + 1 01062 250 CONTINUE 01063 IF( IP.EQ.-1 ) 01064 $ IP = 0 01065 IF( IP.EQ.1 ) 01066 $ IP = -1 01067 * 01068 260 CONTINUE 01069 * 01070 END IF 01071 * 01072 RETURN 01073 * 01074 * End of DTREVC 01075 * 01076 END