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