|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00002 $ LDVR, MM, M, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 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 DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00016 $ WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DTREVC 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 DHSEQR. 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**H)*T = w*(y**H) 00031 * 00032 * where y**H denotes the conjugate 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DHSEQR). 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) DOUBLE PRECISION 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 DHSEQR). 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) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE 00149 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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 DOUBLE PRECISION 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 IDAMAX 00161 DOUBLE PRECISION DDOT, DLAMCH 00162 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA 00166 * .. 00167 * .. Intrinsic Functions .. 00168 INTRINSIC ABS, MAX, SQRT 00169 * .. 00170 * .. Local Arrays .. 00171 DOUBLE PRECISION X( 2, 2 ) 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Decode and test the input parameters 00176 * 00177 BOTHV = LSAME( SIDE, 'B' ) 00178 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 00179 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 00180 * 00181 ALLV = LSAME( HOWMNY, 'A' ) 00182 OVER = LSAME( HOWMNY, 'B' ) 00183 SOMEV = LSAME( HOWMNY, 'S' ) 00184 * 00185 INFO = 0 00186 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 00187 INFO = -1 00188 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 00189 INFO = -2 00190 ELSE IF( N.LT.0 ) THEN 00191 INFO = -4 00192 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00193 INFO = -6 00194 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 00195 INFO = -8 00196 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 00197 INFO = -10 00198 ELSE 00199 * 00200 * Set M to the number of columns required to store the selected 00201 * eigenvectors, standardize the array SELECT if necessary, and 00202 * test MM. 00203 * 00204 IF( SOMEV ) THEN 00205 M = 0 00206 PAIR = .FALSE. 00207 DO 10 J = 1, N 00208 IF( PAIR ) THEN 00209 PAIR = .FALSE. 00210 SELECT( J ) = .FALSE. 00211 ELSE 00212 IF( J.LT.N ) THEN 00213 IF( T( J+1, J ).EQ.ZERO ) THEN 00214 IF( SELECT( J ) ) 00215 $ M = M + 1 00216 ELSE 00217 PAIR = .TRUE. 00218 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN 00219 SELECT( J ) = .TRUE. 00220 M = M + 2 00221 END IF 00222 END IF 00223 ELSE 00224 IF( SELECT( N ) ) 00225 $ M = M + 1 00226 END IF 00227 END IF 00228 10 CONTINUE 00229 ELSE 00230 M = N 00231 END IF 00232 * 00233 IF( MM.LT.M ) THEN 00234 INFO = -11 00235 END IF 00236 END IF 00237 IF( INFO.NE.0 ) THEN 00238 CALL XERBLA( 'DTREVC', -INFO ) 00239 RETURN 00240 END IF 00241 * 00242 * Quick return if possible. 00243 * 00244 IF( N.EQ.0 ) 00245 $ RETURN 00246 * 00247 * Set the constants to control overflow. 00248 * 00249 UNFL = DLAMCH( 'Safe minimum' ) 00250 OVFL = ONE / UNFL 00251 CALL DLABAD( UNFL, OVFL ) 00252 ULP = DLAMCH( 'Precision' ) 00253 SMLNUM = UNFL*( N / ULP ) 00254 BIGNUM = ( ONE-ULP ) / SMLNUM 00255 * 00256 * Compute 1-norm of each column of strictly upper triangular 00257 * part of T to control overflow in triangular solver. 00258 * 00259 WORK( 1 ) = ZERO 00260 DO 30 J = 2, N 00261 WORK( J ) = ZERO 00262 DO 20 I = 1, J - 1 00263 WORK( J ) = WORK( J ) + ABS( T( I, J ) ) 00264 20 CONTINUE 00265 30 CONTINUE 00266 * 00267 * Index IP is used to specify the real or complex eigenvalue: 00268 * IP = 0, real eigenvalue, 00269 * 1, first of conjugate complex pair: (wr,wi) 00270 * -1, second of conjugate complex pair: (wr,wi) 00271 * 00272 N2 = 2*N 00273 * 00274 IF( RIGHTV ) THEN 00275 * 00276 * Compute right eigenvectors. 00277 * 00278 IP = 0 00279 IS = M 00280 DO 140 KI = N, 1, -1 00281 * 00282 IF( IP.EQ.1 ) 00283 $ GO TO 130 00284 IF( KI.EQ.1 ) 00285 $ GO TO 40 00286 IF( T( KI, KI-1 ).EQ.ZERO ) 00287 $ GO TO 40 00288 IP = -1 00289 * 00290 40 CONTINUE 00291 IF( SOMEV ) THEN 00292 IF( IP.EQ.0 ) THEN 00293 IF( .NOT.SELECT( KI ) ) 00294 $ GO TO 130 00295 ELSE 00296 IF( .NOT.SELECT( KI-1 ) ) 00297 $ GO TO 130 00298 END IF 00299 END IF 00300 * 00301 * Compute the KI-th eigenvalue (WR,WI). 00302 * 00303 WR = T( KI, KI ) 00304 WI = ZERO 00305 IF( IP.NE.0 ) 00306 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* 00307 $ SQRT( ABS( T( KI-1, KI ) ) ) 00308 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00309 * 00310 IF( IP.EQ.0 ) THEN 00311 * 00312 * Real right eigenvector 00313 * 00314 WORK( KI+N ) = ONE 00315 * 00316 * Form right-hand side 00317 * 00318 DO 50 K = 1, KI - 1 00319 WORK( K+N ) = -T( K, KI ) 00320 50 CONTINUE 00321 * 00322 * Solve the upper quasi-triangular system: 00323 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. 00324 * 00325 JNXT = KI - 1 00326 DO 60 J = KI - 1, 1, -1 00327 IF( J.GT.JNXT ) 00328 $ GO TO 60 00329 J1 = J 00330 J2 = J 00331 JNXT = J - 1 00332 IF( J.GT.1 ) THEN 00333 IF( T( J, J-1 ).NE.ZERO ) THEN 00334 J1 = J - 1 00335 JNXT = J - 2 00336 END IF 00337 END IF 00338 * 00339 IF( J1.EQ.J2 ) THEN 00340 * 00341 * 1-by-1 diagonal block 00342 * 00343 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00344 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00345 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00346 * 00347 * Scale X(1,1) to avoid overflow when updating 00348 * the right-hand side. 00349 * 00350 IF( XNORM.GT.ONE ) THEN 00351 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00352 X( 1, 1 ) = X( 1, 1 ) / XNORM 00353 SCALE = SCALE / XNORM 00354 END IF 00355 END IF 00356 * 00357 * Scale if necessary 00358 * 00359 IF( SCALE.NE.ONE ) 00360 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00361 WORK( J+N ) = X( 1, 1 ) 00362 * 00363 * Update right-hand side 00364 * 00365 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00366 $ WORK( 1+N ), 1 ) 00367 * 00368 ELSE 00369 * 00370 * 2-by-2 diagonal block 00371 * 00372 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, 00373 $ T( J-1, J-1 ), LDT, ONE, ONE, 00374 $ WORK( J-1+N ), N, WR, ZERO, X, 2, 00375 $ SCALE, XNORM, IERR ) 00376 * 00377 * Scale X(1,1) and X(2,1) to avoid overflow when 00378 * updating the right-hand side. 00379 * 00380 IF( XNORM.GT.ONE ) THEN 00381 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00382 IF( BETA.GT.BIGNUM / XNORM ) THEN 00383 X( 1, 1 ) = X( 1, 1 ) / XNORM 00384 X( 2, 1 ) = X( 2, 1 ) / XNORM 00385 SCALE = SCALE / XNORM 00386 END IF 00387 END IF 00388 * 00389 * Scale if necessary 00390 * 00391 IF( SCALE.NE.ONE ) 00392 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00393 WORK( J-1+N ) = X( 1, 1 ) 00394 WORK( J+N ) = X( 2, 1 ) 00395 * 00396 * Update right-hand side 00397 * 00398 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00399 $ WORK( 1+N ), 1 ) 00400 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00401 $ WORK( 1+N ), 1 ) 00402 END IF 00403 60 CONTINUE 00404 * 00405 * Copy the vector x or Q*x to VR and normalize. 00406 * 00407 IF( .NOT.OVER ) THEN 00408 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) 00409 * 00410 II = IDAMAX( KI, VR( 1, IS ), 1 ) 00411 REMAX = ONE / ABS( VR( II, IS ) ) 00412 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00413 * 00414 DO 70 K = KI + 1, N 00415 VR( K, IS ) = ZERO 00416 70 CONTINUE 00417 ELSE 00418 IF( KI.GT.1 ) 00419 $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR, 00420 $ WORK( 1+N ), 1, WORK( KI+N ), 00421 $ VR( 1, KI ), 1 ) 00422 * 00423 II = IDAMAX( N, VR( 1, KI ), 1 ) 00424 REMAX = ONE / ABS( VR( II, KI ) ) 00425 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 00426 END IF 00427 * 00428 ELSE 00429 * 00430 * Complex right eigenvector. 00431 * 00432 * Initial solve 00433 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. 00434 * [ (T(KI,KI-1) T(KI,KI) ) ] 00435 * 00436 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN 00437 WORK( KI-1+N ) = ONE 00438 WORK( KI+N2 ) = WI / T( KI-1, KI ) 00439 ELSE 00440 WORK( KI-1+N ) = -WI / T( KI, KI-1 ) 00441 WORK( KI+N2 ) = ONE 00442 END IF 00443 WORK( KI+N ) = ZERO 00444 WORK( KI-1+N2 ) = ZERO 00445 * 00446 * Form right-hand side 00447 * 00448 DO 80 K = 1, KI - 2 00449 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) 00450 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) 00451 80 CONTINUE 00452 * 00453 * Solve upper quasi-triangular system: 00454 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) 00455 * 00456 JNXT = KI - 2 00457 DO 90 J = KI - 2, 1, -1 00458 IF( J.GT.JNXT ) 00459 $ GO TO 90 00460 J1 = J 00461 J2 = J 00462 JNXT = J - 1 00463 IF( J.GT.1 ) THEN 00464 IF( T( J, J-1 ).NE.ZERO ) THEN 00465 J1 = J - 1 00466 JNXT = J - 2 00467 END IF 00468 END IF 00469 * 00470 IF( J1.EQ.J2 ) THEN 00471 * 00472 * 1-by-1 diagonal block 00473 * 00474 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00475 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, 00476 $ X, 2, SCALE, XNORM, IERR ) 00477 * 00478 * Scale X(1,1) and X(1,2) to avoid overflow when 00479 * updating the right-hand side. 00480 * 00481 IF( XNORM.GT.ONE ) THEN 00482 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 00483 X( 1, 1 ) = X( 1, 1 ) / XNORM 00484 X( 1, 2 ) = X( 1, 2 ) / XNORM 00485 SCALE = SCALE / XNORM 00486 END IF 00487 END IF 00488 * 00489 * Scale if necessary 00490 * 00491 IF( SCALE.NE.ONE ) THEN 00492 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00493 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00494 END IF 00495 WORK( J+N ) = X( 1, 1 ) 00496 WORK( J+N2 ) = X( 1, 2 ) 00497 * 00498 * Update the right-hand side 00499 * 00500 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 00501 $ WORK( 1+N ), 1 ) 00502 CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, 00503 $ WORK( 1+N2 ), 1 ) 00504 * 00505 ELSE 00506 * 00507 * 2-by-2 diagonal block 00508 * 00509 CALL DLALN2( .FALSE., 2, 2, SMIN, ONE, 00510 $ T( J-1, J-1 ), LDT, ONE, ONE, 00511 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, 00512 $ XNORM, IERR ) 00513 * 00514 * Scale X to avoid overflow when updating 00515 * the right-hand side. 00516 * 00517 IF( XNORM.GT.ONE ) THEN 00518 BETA = MAX( WORK( J-1 ), WORK( J ) ) 00519 IF( BETA.GT.BIGNUM / XNORM ) THEN 00520 REC = ONE / XNORM 00521 X( 1, 1 ) = X( 1, 1 )*REC 00522 X( 1, 2 ) = X( 1, 2 )*REC 00523 X( 2, 1 ) = X( 2, 1 )*REC 00524 X( 2, 2 ) = X( 2, 2 )*REC 00525 SCALE = SCALE*REC 00526 END IF 00527 END IF 00528 * 00529 * Scale if necessary 00530 * 00531 IF( SCALE.NE.ONE ) THEN 00532 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) 00533 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 00534 END IF 00535 WORK( J-1+N ) = X( 1, 1 ) 00536 WORK( J+N ) = X( 2, 1 ) 00537 WORK( J-1+N2 ) = X( 1, 2 ) 00538 WORK( J+N2 ) = X( 2, 2 ) 00539 * 00540 * Update the right-hand side 00541 * 00542 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 00543 $ WORK( 1+N ), 1 ) 00544 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 00545 $ WORK( 1+N ), 1 ) 00546 CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, 00547 $ WORK( 1+N2 ), 1 ) 00548 CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, 00549 $ WORK( 1+N2 ), 1 ) 00550 END IF 00551 90 CONTINUE 00552 * 00553 * Copy the vector x or Q*x to VR and normalize. 00554 * 00555 IF( .NOT.OVER ) THEN 00556 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) 00557 CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) 00558 * 00559 EMAX = ZERO 00560 DO 100 K = 1, KI 00561 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ 00562 $ ABS( VR( K, IS ) ) ) 00563 100 CONTINUE 00564 * 00565 REMAX = ONE / EMAX 00566 CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) 00567 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) 00568 * 00569 DO 110 K = KI + 1, N 00570 VR( K, IS-1 ) = ZERO 00571 VR( K, IS ) = ZERO 00572 110 CONTINUE 00573 * 00574 ELSE 00575 * 00576 IF( KI.GT.2 ) THEN 00577 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00578 $ WORK( 1+N ), 1, WORK( KI-1+N ), 00579 $ VR( 1, KI-1 ), 1 ) 00580 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, 00581 $ WORK( 1+N2 ), 1, WORK( KI+N2 ), 00582 $ VR( 1, KI ), 1 ) 00583 ELSE 00584 CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) 00585 CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) 00586 END IF 00587 * 00588 EMAX = ZERO 00589 DO 120 K = 1, N 00590 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ 00591 $ ABS( VR( K, KI ) ) ) 00592 120 CONTINUE 00593 REMAX = ONE / EMAX 00594 CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) 00595 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) 00596 END IF 00597 END IF 00598 * 00599 IS = IS - 1 00600 IF( IP.NE.0 ) 00601 $ IS = IS - 1 00602 130 CONTINUE 00603 IF( IP.EQ.1 ) 00604 $ IP = 0 00605 IF( IP.EQ.-1 ) 00606 $ IP = 1 00607 140 CONTINUE 00608 END IF 00609 * 00610 IF( LEFTV ) THEN 00611 * 00612 * Compute left eigenvectors. 00613 * 00614 IP = 0 00615 IS = 1 00616 DO 260 KI = 1, N 00617 * 00618 IF( IP.EQ.-1 ) 00619 $ GO TO 250 00620 IF( KI.EQ.N ) 00621 $ GO TO 150 00622 IF( T( KI+1, KI ).EQ.ZERO ) 00623 $ GO TO 150 00624 IP = 1 00625 * 00626 150 CONTINUE 00627 IF( SOMEV ) THEN 00628 IF( .NOT.SELECT( KI ) ) 00629 $ GO TO 250 00630 END IF 00631 * 00632 * Compute the KI-th eigenvalue (WR,WI). 00633 * 00634 WR = T( KI, KI ) 00635 WI = ZERO 00636 IF( IP.NE.0 ) 00637 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* 00638 $ SQRT( ABS( T( KI+1, KI ) ) ) 00639 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 00640 * 00641 IF( IP.EQ.0 ) THEN 00642 * 00643 * Real left eigenvector. 00644 * 00645 WORK( KI+N ) = ONE 00646 * 00647 * Form right-hand side 00648 * 00649 DO 160 K = KI + 1, N 00650 WORK( K+N ) = -T( KI, K ) 00651 160 CONTINUE 00652 * 00653 * Solve the quasi-triangular system: 00654 * (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK 00655 * 00656 VMAX = ONE 00657 VCRIT = BIGNUM 00658 * 00659 JNXT = KI + 1 00660 DO 170 J = KI + 1, N 00661 IF( J.LT.JNXT ) 00662 $ GO TO 170 00663 J1 = J 00664 J2 = J 00665 JNXT = J + 1 00666 IF( J.LT.N ) THEN 00667 IF( T( J+1, J ).NE.ZERO ) THEN 00668 J2 = J + 1 00669 JNXT = J + 2 00670 END IF 00671 END IF 00672 * 00673 IF( J1.EQ.J2 ) THEN 00674 * 00675 * 1-by-1 diagonal block 00676 * 00677 * Scale if necessary to avoid overflow when forming 00678 * the right-hand side. 00679 * 00680 IF( WORK( J ).GT.VCRIT ) THEN 00681 REC = ONE / VMAX 00682 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00683 VMAX = ONE 00684 VCRIT = BIGNUM 00685 END IF 00686 * 00687 WORK( J+N ) = WORK( J+N ) - 00688 $ DDOT( J-KI-1, T( KI+1, J ), 1, 00689 $ WORK( KI+1+N ), 1 ) 00690 * 00691 * Solve (T(J,J)-WR)'*X = WORK 00692 * 00693 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 00694 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00695 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00696 * 00697 * Scale if necessary 00698 * 00699 IF( SCALE.NE.ONE ) 00700 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00701 WORK( J+N ) = X( 1, 1 ) 00702 VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) 00703 VCRIT = BIGNUM / VMAX 00704 * 00705 ELSE 00706 * 00707 * 2-by-2 diagonal block 00708 * 00709 * Scale if necessary to avoid overflow when forming 00710 * the right-hand side. 00711 * 00712 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00713 IF( BETA.GT.VCRIT ) THEN 00714 REC = ONE / VMAX 00715 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00716 VMAX = ONE 00717 VCRIT = BIGNUM 00718 END IF 00719 * 00720 WORK( J+N ) = WORK( J+N ) - 00721 $ DDOT( J-KI-1, T( KI+1, J ), 1, 00722 $ WORK( KI+1+N ), 1 ) 00723 * 00724 WORK( J+1+N ) = WORK( J+1+N ) - 00725 $ DDOT( J-KI-1, T( KI+1, J+1 ), 1, 00726 $ WORK( KI+1+N ), 1 ) 00727 * 00728 * Solve 00729 * [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 ) 00730 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) 00731 * 00732 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), 00733 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00734 $ ZERO, X, 2, SCALE, XNORM, IERR ) 00735 * 00736 * Scale if necessary 00737 * 00738 IF( SCALE.NE.ONE ) 00739 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00740 WORK( J+N ) = X( 1, 1 ) 00741 WORK( J+1+N ) = X( 2, 1 ) 00742 * 00743 VMAX = MAX( ABS( WORK( J+N ) ), 00744 $ ABS( WORK( J+1+N ) ), VMAX ) 00745 VCRIT = BIGNUM / VMAX 00746 * 00747 END IF 00748 170 CONTINUE 00749 * 00750 * Copy the vector x or Q*x to VL and normalize. 00751 * 00752 IF( .NOT.OVER ) THEN 00753 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 00754 * 00755 II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 00756 REMAX = ONE / ABS( VL( II, IS ) ) 00757 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00758 * 00759 DO 180 K = 1, KI - 1 00760 VL( K, IS ) = ZERO 00761 180 CONTINUE 00762 * 00763 ELSE 00764 * 00765 IF( KI.LT.N ) 00766 $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, 00767 $ WORK( KI+1+N ), 1, WORK( KI+N ), 00768 $ VL( 1, KI ), 1 ) 00769 * 00770 II = IDAMAX( N, VL( 1, KI ), 1 ) 00771 REMAX = ONE / ABS( VL( II, KI ) ) 00772 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 00773 * 00774 END IF 00775 * 00776 ELSE 00777 * 00778 * Complex left eigenvector. 00779 * 00780 * Initial solve: 00781 * ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0. 00782 * ((T(KI+1,KI) T(KI+1,KI+1)) ) 00783 * 00784 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN 00785 WORK( KI+N ) = WI / T( KI, KI+1 ) 00786 WORK( KI+1+N2 ) = ONE 00787 ELSE 00788 WORK( KI+N ) = ONE 00789 WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) 00790 END IF 00791 WORK( KI+1+N ) = ZERO 00792 WORK( KI+N2 ) = ZERO 00793 * 00794 * Form right-hand side 00795 * 00796 DO 190 K = KI + 2, N 00797 WORK( K+N ) = -WORK( KI+N )*T( KI, K ) 00798 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) 00799 190 CONTINUE 00800 * 00801 * Solve complex quasi-triangular system: 00802 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 00803 * 00804 VMAX = ONE 00805 VCRIT = BIGNUM 00806 * 00807 JNXT = KI + 2 00808 DO 200 J = KI + 2, N 00809 IF( J.LT.JNXT ) 00810 $ GO TO 200 00811 J1 = J 00812 J2 = J 00813 JNXT = J + 1 00814 IF( J.LT.N ) THEN 00815 IF( T( J+1, J ).NE.ZERO ) THEN 00816 J2 = J + 1 00817 JNXT = J + 2 00818 END IF 00819 END IF 00820 * 00821 IF( J1.EQ.J2 ) THEN 00822 * 00823 * 1-by-1 diagonal block 00824 * 00825 * Scale if necessary to avoid overflow when 00826 * forming the right-hand side elements. 00827 * 00828 IF( WORK( J ).GT.VCRIT ) THEN 00829 REC = ONE / VMAX 00830 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00831 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00832 VMAX = ONE 00833 VCRIT = BIGNUM 00834 END IF 00835 * 00836 WORK( J+N ) = WORK( J+N ) - 00837 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00838 $ WORK( KI+2+N ), 1 ) 00839 WORK( J+N2 ) = WORK( J+N2 ) - 00840 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00841 $ WORK( KI+2+N2 ), 1 ) 00842 * 00843 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 00844 * 00845 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 00846 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00847 $ -WI, X, 2, SCALE, XNORM, IERR ) 00848 * 00849 * Scale if necessary 00850 * 00851 IF( SCALE.NE.ONE ) THEN 00852 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00853 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 00854 END IF 00855 WORK( J+N ) = X( 1, 1 ) 00856 WORK( J+N2 ) = X( 1, 2 ) 00857 VMAX = MAX( ABS( WORK( J+N ) ), 00858 $ ABS( WORK( J+N2 ) ), VMAX ) 00859 VCRIT = BIGNUM / VMAX 00860 * 00861 ELSE 00862 * 00863 * 2-by-2 diagonal block 00864 * 00865 * Scale if necessary to avoid overflow when forming 00866 * the right-hand side elements. 00867 * 00868 BETA = MAX( WORK( J ), WORK( J+1 ) ) 00869 IF( BETA.GT.VCRIT ) THEN 00870 REC = ONE / VMAX 00871 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 00872 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 00873 VMAX = ONE 00874 VCRIT = BIGNUM 00875 END IF 00876 * 00877 WORK( J+N ) = WORK( J+N ) - 00878 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00879 $ WORK( KI+2+N ), 1 ) 00880 * 00881 WORK( J+N2 ) = WORK( J+N2 ) - 00882 $ DDOT( J-KI-2, T( KI+2, J ), 1, 00883 $ WORK( KI+2+N2 ), 1 ) 00884 * 00885 WORK( J+1+N ) = WORK( J+1+N ) - 00886 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 00887 $ WORK( KI+2+N ), 1 ) 00888 * 00889 WORK( J+1+N2 ) = WORK( J+1+N2 ) - 00890 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, 00891 $ WORK( KI+2+N2 ), 1 ) 00892 * 00893 * Solve 2-by-2 complex linear equation 00894 * ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B 00895 * ([T(j+1,j) T(j+1,j+1)] ) 00896 * 00897 CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), 00898 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 00899 $ -WI, X, 2, SCALE, XNORM, IERR ) 00900 * 00901 * Scale if necessary 00902 * 00903 IF( SCALE.NE.ONE ) THEN 00904 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 00905 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 00906 END IF 00907 WORK( J+N ) = X( 1, 1 ) 00908 WORK( J+N2 ) = X( 1, 2 ) 00909 WORK( J+1+N ) = X( 2, 1 ) 00910 WORK( J+1+N2 ) = X( 2, 2 ) 00911 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), 00912 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) 00913 VCRIT = BIGNUM / VMAX 00914 * 00915 END IF 00916 200 CONTINUE 00917 * 00918 * Copy the vector x or Q*x to VL and normalize. 00919 * 00920 IF( .NOT.OVER ) THEN 00921 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 00922 CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), 00923 $ 1 ) 00924 * 00925 EMAX = ZERO 00926 DO 220 K = KI, N 00927 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ 00928 $ ABS( VL( K, IS+1 ) ) ) 00929 220 CONTINUE 00930 REMAX = ONE / EMAX 00931 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 00932 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) 00933 * 00934 DO 230 K = 1, KI - 1 00935 VL( K, IS ) = ZERO 00936 VL( K, IS+1 ) = ZERO 00937 230 CONTINUE 00938 ELSE 00939 IF( KI.LT.N-1 ) THEN 00940 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 00941 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), 00942 $ VL( 1, KI ), 1 ) 00943 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 00944 $ LDVL, WORK( KI+2+N2 ), 1, 00945 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 00946 ELSE 00947 CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) 00948 CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 00949 END IF 00950 * 00951 EMAX = ZERO 00952 DO 240 K = 1, N 00953 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ 00954 $ ABS( VL( K, KI+1 ) ) ) 00955 240 CONTINUE 00956 REMAX = ONE / EMAX 00957 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) 00958 CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) 00959 * 00960 END IF 00961 * 00962 END IF 00963 * 00964 IS = IS + 1 00965 IF( IP.NE.0 ) 00966 $ IS = IS + 1 00967 250 CONTINUE 00968 IF( IP.EQ.-1 ) 00969 $ IP = 0 00970 IF( IP.EQ.1 ) 00971 $ IP = -1 00972 * 00973 260 CONTINUE 00974 * 00975 END IF 00976 * 00977 RETURN 00978 * 00979 * End of DTREVC 00980 * 00981 END
1.7.3 
