 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 00002 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, 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, LDP, LDS, LDVL, LDVR, M, MM, N 00012 * .. 00013 * .. Array Arguments .. 00014 LOGICAL SELECT( * ) 00015 REAL RWORK( * ) 00016 COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 00017 $ VR( LDVR, * ), WORK( * ) 00018 * .. 00019 * 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CTGEVC computes some or all of the right and/or left eigenvectors of 00025 * a pair of complex matrices (S,P), where S and P are upper triangular. 00026 * Matrix pairs of this type are produced by the generalized Schur 00027 * factorization of a complex matrix pair (A,B): 00028 * 00029 * A = Q*S*Z**H, B = Q*P*Z**H 00030 * 00031 * as computed by CGGHRD + CHGEQZ. 00032 * 00033 * The right eigenvector x and the left eigenvector y of (S,P) 00034 * corresponding to an eigenvalue w are defined by: 00035 * 00036 * S*x = w*P*x, (y**H)*S = w*(y**H)*P, 00037 * 00038 * where y**H denotes the conjugate tranpose of y. 00039 * The eigenvalues are not input to this routine, but are computed 00040 * directly from the diagonal elements of S and P. 00041 * 00042 * This routine returns the matrices X and/or Y of right and left 00043 * eigenvectors of (S,P), or the products Z*X and/or Q*Y, 00044 * where Z and Q are input matrices. 00045 * If Q and Z are the unitary factors from the generalized Schur 00046 * factorization of a matrix pair (A,B), then Z*X and Q*Y 00047 * are the matrices of right and left eigenvectors of (A,B). 00048 * 00049 * Arguments 00050 * ========= 00051 * 00052 * SIDE (input) CHARACTER*1 00053 * = 'R': compute right eigenvectors only; 00054 * = 'L': compute left eigenvectors only; 00055 * = 'B': compute both right and left eigenvectors. 00056 * 00057 * HOWMNY (input) CHARACTER*1 00058 * = 'A': compute all right and/or left eigenvectors; 00059 * = 'B': compute all right and/or left eigenvectors, 00060 * backtransformed by the matrices in VR and/or VL; 00061 * = 'S': compute selected right and/or left eigenvectors, 00062 * specified by the logical array SELECT. 00063 * 00064 * SELECT (input) LOGICAL array, dimension (N) 00065 * If HOWMNY='S', SELECT specifies the eigenvectors to be 00066 * computed. The eigenvector corresponding to the j-th 00067 * eigenvalue is computed if SELECT(j) = .TRUE.. 00068 * Not referenced if HOWMNY = 'A' or 'B'. 00069 * 00070 * N (input) INTEGER 00071 * The order of the matrices S and P. N >= 0. 00072 * 00073 * S (input) COMPLEX array, dimension (LDS,N) 00074 * The upper triangular matrix S from a generalized Schur 00075 * factorization, as computed by CHGEQZ. 00076 * 00077 * LDS (input) INTEGER 00078 * The leading dimension of array S. LDS >= max(1,N). 00079 * 00080 * P (input) COMPLEX array, dimension (LDP,N) 00081 * The upper triangular matrix P from a generalized Schur 00082 * factorization, as computed by CHGEQZ. P must have real 00083 * diagonal elements. 00084 * 00085 * LDP (input) INTEGER 00086 * The leading dimension of array P. LDP >= max(1,N). 00087 * 00088 * VL (input/output) COMPLEX array, dimension (LDVL,MM) 00089 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00090 * contain an N-by-N matrix Q (usually the unitary matrix Q 00091 * of left Schur vectors returned by CHGEQZ). 00092 * On exit, if SIDE = 'L' or 'B', VL contains: 00093 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); 00094 * if HOWMNY = 'B', the matrix Q*Y; 00095 * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by 00096 * SELECT, stored consecutively in the columns of 00097 * VL, in the same order as their eigenvalues. 00098 * Not referenced if SIDE = 'R'. 00099 * 00100 * LDVL (input) INTEGER 00101 * The leading dimension of array VL. LDVL >= 1, and if 00102 * SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N. 00103 * 00104 * VR (input/output) COMPLEX array, dimension (LDVR,MM) 00105 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00106 * contain an N-by-N matrix Q (usually the unitary matrix Z 00107 * of right Schur vectors returned by CHGEQZ). 00108 * On exit, if SIDE = 'R' or 'B', VR contains: 00109 * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); 00110 * if HOWMNY = 'B', the matrix Z*X; 00111 * if HOWMNY = 'S', the right eigenvectors of (S,P) specified by 00112 * SELECT, stored consecutively in the columns of 00113 * VR, in the same order as their eigenvalues. 00114 * Not referenced if SIDE = 'L'. 00115 * 00116 * LDVR (input) INTEGER 00117 * The leading dimension of the array VR. LDVR >= 1, and if 00118 * SIDE = 'R' or 'B', LDVR >= N. 00119 * 00120 * MM (input) INTEGER 00121 * The number of columns in the arrays VL and/or VR. MM >= M. 00122 * 00123 * M (output) INTEGER 00124 * The number of columns in the arrays VL and/or VR actually 00125 * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 00126 * is set to N. Each selected eigenvector occupies one column. 00127 * 00128 * WORK (workspace) COMPLEX array, dimension (2*N) 00129 * 00130 * RWORK (workspace) REAL array, dimension (2*N) 00131 * 00132 * INFO (output) INTEGER 00133 * = 0: successful exit. 00134 * < 0: if INFO = -i, the i-th argument had an illegal value. 00135 * 00136 * ===================================================================== 00137 * 00138 * .. Parameters .. 00139 REAL ZERO, ONE 00140 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00141 COMPLEX CZERO, CONE 00142 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00143 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00144 * .. 00145 * .. Local Scalars .. 00146 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP, 00147 $ LSA, LSB 00148 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC, 00149 $ J, JE, JR 00150 REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG, 00151 $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA, 00152 $ SCALE, SMALL, TEMP, ULP, XMAX 00153 COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X 00154 * .. 00155 * .. External Functions .. 00156 LOGICAL LSAME 00157 REAL SLAMCH 00158 COMPLEX CLADIV 00159 EXTERNAL LSAME, SLAMCH, CLADIV 00160 * .. 00161 * .. External Subroutines .. 00162 EXTERNAL CGEMV, SLABAD, XERBLA 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL 00166 * .. 00167 * .. Statement Functions .. 00168 REAL ABS1 00169 * .. 00170 * .. Statement Function definitions .. 00171 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Decode and Test the input parameters 00176 * 00177 IF( LSAME( HOWMNY, 'A' ) ) THEN 00178 IHWMNY = 1 00179 ILALL = .TRUE. 00180 ILBACK = .FALSE. 00181 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 00182 IHWMNY = 2 00183 ILALL = .FALSE. 00184 ILBACK = .FALSE. 00185 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN 00186 IHWMNY = 3 00187 ILALL = .TRUE. 00188 ILBACK = .TRUE. 00189 ELSE 00190 IHWMNY = -1 00191 END IF 00192 * 00193 IF( LSAME( SIDE, 'R' ) ) THEN 00194 ISIDE = 1 00195 COMPL = .FALSE. 00196 COMPR = .TRUE. 00197 ELSE IF( LSAME( SIDE, 'L' ) ) THEN 00198 ISIDE = 2 00199 COMPL = .TRUE. 00200 COMPR = .FALSE. 00201 ELSE IF( LSAME( SIDE, 'B' ) ) THEN 00202 ISIDE = 3 00203 COMPL = .TRUE. 00204 COMPR = .TRUE. 00205 ELSE 00206 ISIDE = -1 00207 END IF 00208 * 00209 INFO = 0 00210 IF( ISIDE.LT.0 ) THEN 00211 INFO = -1 00212 ELSE IF( IHWMNY.LT.0 ) THEN 00213 INFO = -2 00214 ELSE IF( N.LT.0 ) THEN 00215 INFO = -4 00216 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN 00217 INFO = -6 00218 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN 00219 INFO = -8 00220 END IF 00221 IF( INFO.NE.0 ) THEN 00222 CALL XERBLA( 'CTGEVC', -INFO ) 00223 RETURN 00224 END IF 00225 * 00226 * Count the number of eigenvectors 00227 * 00228 IF( .NOT.ILALL ) THEN 00229 IM = 0 00230 DO 10 J = 1, N 00231 IF( SELECT( J ) ) 00232 $ IM = IM + 1 00233 10 CONTINUE 00234 ELSE 00235 IM = N 00236 END IF 00237 * 00238 * Check diagonal of B 00239 * 00240 ILBBAD = .FALSE. 00241 DO 20 J = 1, N 00242 IF( AIMAG( P( J, J ) ).NE.ZERO ) 00243 $ ILBBAD = .TRUE. 00244 20 CONTINUE 00245 * 00246 IF( ILBBAD ) THEN 00247 INFO = -7 00248 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 00249 INFO = -10 00250 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 00251 INFO = -12 00252 ELSE IF( MM.LT.IM ) THEN 00253 INFO = -13 00254 END IF 00255 IF( INFO.NE.0 ) THEN 00256 CALL XERBLA( 'CTGEVC', -INFO ) 00257 RETURN 00258 END IF 00259 * 00260 * Quick return if possible 00261 * 00262 M = IM 00263 IF( N.EQ.0 ) 00264 $ RETURN 00265 * 00266 * Machine Constants 00267 * 00268 SAFMIN = SLAMCH( 'Safe minimum' ) 00269 BIG = ONE / SAFMIN 00270 CALL SLABAD( SAFMIN, BIG ) 00271 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00272 SMALL = SAFMIN*N / ULP 00273 BIG = ONE / SMALL 00274 BIGNUM = ONE / ( SAFMIN*N ) 00275 * 00276 * Compute the 1-norm of each column of the strictly upper triangular 00277 * part of A and B to check for possible overflow in the triangular 00278 * solver. 00279 * 00280 ANORM = ABS1( S( 1, 1 ) ) 00281 BNORM = ABS1( P( 1, 1 ) ) 00282 RWORK( 1 ) = ZERO 00283 RWORK( N+1 ) = ZERO 00284 DO 40 J = 2, N 00285 RWORK( J ) = ZERO 00286 RWORK( N+J ) = ZERO 00287 DO 30 I = 1, J - 1 00288 RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) ) 00289 RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) ) 00290 30 CONTINUE 00291 ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) ) 00292 BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) ) 00293 40 CONTINUE 00294 * 00295 ASCALE = ONE / MAX( ANORM, SAFMIN ) 00296 BSCALE = ONE / MAX( BNORM, SAFMIN ) 00297 * 00298 * Left eigenvectors 00299 * 00300 IF( COMPL ) THEN 00301 IEIG = 0 00302 * 00303 * Main loop over eigenvalues 00304 * 00305 DO 140 JE = 1, N 00306 IF( ILALL ) THEN 00307 ILCOMP = .TRUE. 00308 ELSE 00309 ILCOMP = SELECT( JE ) 00310 END IF 00311 IF( ILCOMP ) THEN 00312 IEIG = IEIG + 1 00313 * 00314 IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. 00315 $ ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN 00316 * 00317 * Singular matrix pencil -- return unit eigenvector 00318 * 00319 DO 50 JR = 1, N 00320 VL( JR, IEIG ) = CZERO 00321 50 CONTINUE 00322 VL( IEIG, IEIG ) = CONE 00323 GO TO 140 00324 END IF 00325 * 00326 * Non-singular eigenvalue: 00327 * Compute coefficients a and b in 00328 * H 00329 * y ( a A - b B ) = 0 00330 * 00331 TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, 00332 $ ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN ) 00333 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE 00334 SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE 00335 ACOEFF = SBETA*ASCALE 00336 BCOEFF = SALPHA*BSCALE 00337 * 00338 * Scale to avoid underflow 00339 * 00340 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL 00341 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. 00342 $ SMALL 00343 * 00344 SCALE = ONE 00345 IF( LSA ) 00346 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00347 IF( LSB ) 00348 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* 00349 $ MIN( BNORM, BIG ) ) 00350 IF( LSA .OR. LSB ) THEN 00351 SCALE = MIN( SCALE, ONE / 00352 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), 00353 $ ABS1( BCOEFF ) ) ) ) 00354 IF( LSA ) THEN 00355 ACOEFF = ASCALE*( SCALE*SBETA ) 00356 ELSE 00357 ACOEFF = SCALE*ACOEFF 00358 END IF 00359 IF( LSB ) THEN 00360 BCOEFF = BSCALE*( SCALE*SALPHA ) 00361 ELSE 00362 BCOEFF = SCALE*BCOEFF 00363 END IF 00364 END IF 00365 * 00366 ACOEFA = ABS( ACOEFF ) 00367 BCOEFA = ABS1( BCOEFF ) 00368 XMAX = ONE 00369 DO 60 JR = 1, N 00370 WORK( JR ) = CZERO 00371 60 CONTINUE 00372 WORK( JE ) = CONE 00373 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00374 * 00375 * H 00376 * Triangular solve of (a A - b B) y = 0 00377 * 00378 * H 00379 * (rowwise in (a A - b B) , or columnwise in a A - b B) 00380 * 00381 DO 100 J = JE + 1, N 00382 * 00383 * Compute 00384 * j-1 00385 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) 00386 * k=je 00387 * (Scale if necessary) 00388 * 00389 TEMP = ONE / XMAX 00390 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM* 00391 $ TEMP ) THEN 00392 DO 70 JR = JE, J - 1 00393 WORK( JR ) = TEMP*WORK( JR ) 00394 70 CONTINUE 00395 XMAX = ONE 00396 END IF 00397 SUMA = CZERO 00398 SUMB = CZERO 00399 * 00400 DO 80 JR = JE, J - 1 00401 SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR ) 00402 SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR ) 00403 80 CONTINUE 00404 SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB 00405 * 00406 * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) ) 00407 * 00408 * with scaling and perturbation of the denominator 00409 * 00410 D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) ) 00411 IF( ABS1( D ).LE.DMIN ) 00412 $ D = CMPLX( DMIN ) 00413 * 00414 IF( ABS1( D ).LT.ONE ) THEN 00415 IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN 00416 TEMP = ONE / ABS1( SUM ) 00417 DO 90 JR = JE, J - 1 00418 WORK( JR ) = TEMP*WORK( JR ) 00419 90 CONTINUE 00420 XMAX = TEMP*XMAX 00421 SUM = TEMP*SUM 00422 END IF 00423 END IF 00424 WORK( J ) = CLADIV( -SUM, D ) 00425 XMAX = MAX( XMAX, ABS1( WORK( J ) ) ) 00426 100 CONTINUE 00427 * 00428 * Back transform eigenvector if HOWMNY='B'. 00429 * 00430 IF( ILBACK ) THEN 00431 CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL, 00432 $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 ) 00433 ISRC = 2 00434 IBEG = 1 00435 ELSE 00436 ISRC = 1 00437 IBEG = JE 00438 END IF 00439 * 00440 * Copy and scale eigenvector into column of VL 00441 * 00442 XMAX = ZERO 00443 DO 110 JR = IBEG, N 00444 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) 00445 110 CONTINUE 00446 * 00447 IF( XMAX.GT.SAFMIN ) THEN 00448 TEMP = ONE / XMAX 00449 DO 120 JR = IBEG, N 00450 VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) 00451 120 CONTINUE 00452 ELSE 00453 IBEG = N + 1 00454 END IF 00455 * 00456 DO 130 JR = 1, IBEG - 1 00457 VL( JR, IEIG ) = CZERO 00458 130 CONTINUE 00459 * 00460 END IF 00461 140 CONTINUE 00462 END IF 00463 * 00464 * Right eigenvectors 00465 * 00466 IF( COMPR ) THEN 00467 IEIG = IM + 1 00468 * 00469 * Main loop over eigenvalues 00470 * 00471 DO 250 JE = N, 1, -1 00472 IF( ILALL ) THEN 00473 ILCOMP = .TRUE. 00474 ELSE 00475 ILCOMP = SELECT( JE ) 00476 END IF 00477 IF( ILCOMP ) THEN 00478 IEIG = IEIG - 1 00479 * 00480 IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND. 00481 $ ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN 00482 * 00483 * Singular matrix pencil -- return unit eigenvector 00484 * 00485 DO 150 JR = 1, N 00486 VR( JR, IEIG ) = CZERO 00487 150 CONTINUE 00488 VR( IEIG, IEIG ) = CONE 00489 GO TO 250 00490 END IF 00491 * 00492 * Non-singular eigenvalue: 00493 * Compute coefficients a and b in 00494 * 00495 * ( a A - b B ) x = 0 00496 * 00497 TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE, 00498 $ ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN ) 00499 SALPHA = ( TEMP*S( JE, JE ) )*ASCALE 00500 SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE 00501 ACOEFF = SBETA*ASCALE 00502 BCOEFF = SALPHA*BSCALE 00503 * 00504 * Scale to avoid underflow 00505 * 00506 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL 00507 LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT. 00508 $ SMALL 00509 * 00510 SCALE = ONE 00511 IF( LSA ) 00512 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00513 IF( LSB ) 00514 $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )* 00515 $ MIN( BNORM, BIG ) ) 00516 IF( LSA .OR. LSB ) THEN 00517 SCALE = MIN( SCALE, ONE / 00518 $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ), 00519 $ ABS1( BCOEFF ) ) ) ) 00520 IF( LSA ) THEN 00521 ACOEFF = ASCALE*( SCALE*SBETA ) 00522 ELSE 00523 ACOEFF = SCALE*ACOEFF 00524 END IF 00525 IF( LSB ) THEN 00526 BCOEFF = BSCALE*( SCALE*SALPHA ) 00527 ELSE 00528 BCOEFF = SCALE*BCOEFF 00529 END IF 00530 END IF 00531 * 00532 ACOEFA = ABS( ACOEFF ) 00533 BCOEFA = ABS1( BCOEFF ) 00534 XMAX = ONE 00535 DO 160 JR = 1, N 00536 WORK( JR ) = CZERO 00537 160 CONTINUE 00538 WORK( JE ) = CONE 00539 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00540 * 00541 * Triangular solve of (a A - b B) x = 0 (columnwise) 00542 * 00543 * WORK(1:j-1) contains sums w, 00544 * WORK(j+1:JE) contains x 00545 * 00546 DO 170 JR = 1, JE - 1 00547 WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE ) 00548 170 CONTINUE 00549 WORK( JE ) = CONE 00550 * 00551 DO 210 J = JE - 1, 1, -1 00552 * 00553 * Form x(j) := - w(j) / d 00554 * with scaling and perturbation of the denominator 00555 * 00556 D = ACOEFF*S( J, J ) - BCOEFF*P( J, J ) 00557 IF( ABS1( D ).LE.DMIN ) 00558 $ D = CMPLX( DMIN ) 00559 * 00560 IF( ABS1( D ).LT.ONE ) THEN 00561 IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN 00562 TEMP = ONE / ABS1( WORK( J ) ) 00563 DO 180 JR = 1, JE 00564 WORK( JR ) = TEMP*WORK( JR ) 00565 180 CONTINUE 00566 END IF 00567 END IF 00568 * 00569 WORK( J ) = CLADIV( -WORK( J ), D ) 00570 * 00571 IF( J.GT.1 ) THEN 00572 * 00573 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 00574 * 00575 IF( ABS1( WORK( J ) ).GT.ONE ) THEN 00576 TEMP = ONE / ABS1( WORK( J ) ) 00577 IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE. 00578 $ BIGNUM*TEMP ) THEN 00579 DO 190 JR = 1, JE 00580 WORK( JR ) = TEMP*WORK( JR ) 00581 190 CONTINUE 00582 END IF 00583 END IF 00584 * 00585 CA = ACOEFF*WORK( J ) 00586 CB = BCOEFF*WORK( J ) 00587 DO 200 JR = 1, J - 1 00588 WORK( JR ) = WORK( JR ) + CA*S( JR, J ) - 00589 $ CB*P( JR, J ) 00590 200 CONTINUE 00591 END IF 00592 210 CONTINUE 00593 * 00594 * Back transform eigenvector if HOWMNY='B'. 00595 * 00596 IF( ILBACK ) THEN 00597 CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1, 00598 $ CZERO, WORK( N+1 ), 1 ) 00599 ISRC = 2 00600 IEND = N 00601 ELSE 00602 ISRC = 1 00603 IEND = JE 00604 END IF 00605 * 00606 * Copy and scale eigenvector into column of VR 00607 * 00608 XMAX = ZERO 00609 DO 220 JR = 1, IEND 00610 XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) ) 00611 220 CONTINUE 00612 * 00613 IF( XMAX.GT.SAFMIN ) THEN 00614 TEMP = ONE / XMAX 00615 DO 230 JR = 1, IEND 00616 VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR ) 00617 230 CONTINUE 00618 ELSE 00619 IEND = 0 00620 END IF 00621 * 00622 DO 240 JR = IEND + 1, N 00623 VR( JR, IEIG ) = CZERO 00624 240 CONTINUE 00625 * 00626 END IF 00627 250 CONTINUE 00628 END IF 00629 * 00630 RETURN 00631 * 00632 * End of CTGEVC 00633 * 00634 END
1.7.3 
