![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGEGV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgegv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgegv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgegv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, 00022 * VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBVL, JOBVR 00026 * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL RWORK( * ) 00030 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00031 * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), 00032 * $ WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> This routine is deprecated and has been replaced by routine CGGEV. 00042 *> 00043 *> CGEGV computes the eigenvalues and, optionally, the left and/or right 00044 *> eigenvectors of a complex matrix pair (A,B). 00045 *> Given two square matrices A and B, 00046 *> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the 00047 *> eigenvalues lambda and corresponding (non-zero) eigenvectors x such 00048 *> that 00049 *> A*x = lambda*B*x. 00050 *> 00051 *> An alternate form is to find the eigenvalues mu and corresponding 00052 *> eigenvectors y such that 00053 *> mu*A*y = B*y. 00054 *> 00055 *> These two forms are equivalent with mu = 1/lambda and x = y if 00056 *> neither lambda nor mu is zero. In order to deal with the case that 00057 *> lambda or mu is zero or small, two values alpha and beta are returned 00058 *> for each eigenvalue, such that lambda = alpha/beta and 00059 *> mu = beta/alpha. 00060 *> 00061 *> The vectors x and y in the above equations are right eigenvectors of 00062 *> the matrix pair (A,B). Vectors u and v satisfying 00063 *> u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B 00064 *> are left eigenvectors of (A,B). 00065 *> 00066 *> Note: this routine performs "full balancing" on A and B 00067 *> \endverbatim 00068 * 00069 * Arguments: 00070 * ========== 00071 * 00072 *> \param[in] JOBVL 00073 *> \verbatim 00074 *> JOBVL is CHARACTER*1 00075 *> = 'N': do not compute the left generalized eigenvectors; 00076 *> = 'V': compute the left generalized eigenvectors (returned 00077 *> in VL). 00078 *> \endverbatim 00079 *> 00080 *> \param[in] JOBVR 00081 *> \verbatim 00082 *> JOBVR is CHARACTER*1 00083 *> = 'N': do not compute the right generalized eigenvectors; 00084 *> = 'V': compute the right generalized eigenvectors (returned 00085 *> in VR). 00086 *> \endverbatim 00087 *> 00088 *> \param[in] N 00089 *> \verbatim 00090 *> N is INTEGER 00091 *> The order of the matrices A, B, VL, and VR. N >= 0. 00092 *> \endverbatim 00093 *> 00094 *> \param[in,out] A 00095 *> \verbatim 00096 *> A is COMPLEX array, dimension (LDA, N) 00097 *> On entry, the matrix A. 00098 *> If JOBVL = 'V' or JOBVR = 'V', then on exit A 00099 *> contains the Schur form of A from the generalized Schur 00100 *> factorization of the pair (A,B) after balancing. If no 00101 *> eigenvectors were computed, then only the diagonal elements 00102 *> of the Schur form will be correct. See CGGHRD and CHGEQZ 00103 *> for details. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDA 00107 *> \verbatim 00108 *> LDA is INTEGER 00109 *> The leading dimension of A. LDA >= max(1,N). 00110 *> \endverbatim 00111 *> 00112 *> \param[in,out] B 00113 *> \verbatim 00114 *> B is COMPLEX array, dimension (LDB, N) 00115 *> On entry, the matrix B. 00116 *> If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the 00117 *> upper triangular matrix obtained from B in the generalized 00118 *> Schur factorization of the pair (A,B) after balancing. 00119 *> If no eigenvectors were computed, then only the diagonal 00120 *> elements of B will be correct. See CGGHRD and CHGEQZ for 00121 *> details. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] LDB 00125 *> \verbatim 00126 *> LDB is INTEGER 00127 *> The leading dimension of B. LDB >= max(1,N). 00128 *> \endverbatim 00129 *> 00130 *> \param[out] ALPHA 00131 *> \verbatim 00132 *> ALPHA is COMPLEX array, dimension (N) 00133 *> The complex scalars alpha that define the eigenvalues of 00134 *> GNEP. 00135 *> \endverbatim 00136 *> 00137 *> \param[out] BETA 00138 *> \verbatim 00139 *> BETA is COMPLEX array, dimension (N) 00140 *> The complex scalars beta that define the eigenvalues of GNEP. 00141 *> 00142 *> Together, the quantities alpha = ALPHA(j) and beta = BETA(j) 00143 *> represent the j-th eigenvalue of the matrix pair (A,B), in 00144 *> one of the forms lambda = alpha/beta or mu = beta/alpha. 00145 *> Since either lambda or mu may overflow, they should not, 00146 *> in general, be computed. 00147 *> \endverbatim 00148 *> 00149 *> \param[out] VL 00150 *> \verbatim 00151 *> VL is COMPLEX array, dimension (LDVL,N) 00152 *> If JOBVL = 'V', the left eigenvectors u(j) are stored 00153 *> in the columns of VL, in the same order as their eigenvalues. 00154 *> Each eigenvector is scaled so that its largest component has 00155 *> abs(real part) + abs(imag. part) = 1, except for eigenvectors 00156 *> corresponding to an eigenvalue with alpha = beta = 0, which 00157 *> are set to zero. 00158 *> Not referenced if JOBVL = 'N'. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LDVL 00162 *> \verbatim 00163 *> LDVL is INTEGER 00164 *> The leading dimension of the matrix VL. LDVL >= 1, and 00165 *> if JOBVL = 'V', LDVL >= N. 00166 *> \endverbatim 00167 *> 00168 *> \param[out] VR 00169 *> \verbatim 00170 *> VR is COMPLEX array, dimension (LDVR,N) 00171 *> If JOBVR = 'V', the right eigenvectors x(j) are stored 00172 *> in the columns of VR, in the same order as their eigenvalues. 00173 *> Each eigenvector is scaled so that its largest component has 00174 *> abs(real part) + abs(imag. part) = 1, except for eigenvectors 00175 *> corresponding to an eigenvalue with alpha = beta = 0, which 00176 *> are set to zero. 00177 *> Not referenced if JOBVR = 'N'. 00178 *> \endverbatim 00179 *> 00180 *> \param[in] LDVR 00181 *> \verbatim 00182 *> LDVR is INTEGER 00183 *> The leading dimension of the matrix VR. LDVR >= 1, and 00184 *> if JOBVR = 'V', LDVR >= N. 00185 *> \endverbatim 00186 *> 00187 *> \param[out] WORK 00188 *> \verbatim 00189 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00190 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00191 *> \endverbatim 00192 *> 00193 *> \param[in] LWORK 00194 *> \verbatim 00195 *> LWORK is INTEGER 00196 *> The dimension of the array WORK. LWORK >= max(1,2*N). 00197 *> For good performance, LWORK must generally be larger. 00198 *> To compute the optimal value of LWORK, call ILAENV to get 00199 *> blocksizes (for CGEQRF, CUNMQR, and CUNGQR.) Then compute: 00200 *> NB -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR; 00201 *> The optimal LWORK is MAX( 2*N, N*(NB+1) ). 00202 *> 00203 *> If LWORK = -1, then a workspace query is assumed; the routine 00204 *> only calculates the optimal size of the WORK array, returns 00205 *> this value as the first entry of the WORK array, and no error 00206 *> message related to LWORK is issued by XERBLA. 00207 *> \endverbatim 00208 *> 00209 *> \param[out] RWORK 00210 *> \verbatim 00211 *> RWORK is REAL array, dimension (8*N) 00212 *> \endverbatim 00213 *> 00214 *> \param[out] INFO 00215 *> \verbatim 00216 *> INFO is INTEGER 00217 *> = 0: successful exit 00218 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00219 *> =1,...,N: 00220 *> The QZ iteration failed. No eigenvectors have been 00221 *> calculated, but ALPHA(j) and BETA(j) should be 00222 *> correct for j=INFO+1,...,N. 00223 *> > N: errors that usually indicate LAPACK problems: 00224 *> =N+1: error return from CGGBAL 00225 *> =N+2: error return from CGEQRF 00226 *> =N+3: error return from CUNMQR 00227 *> =N+4: error return from CUNGQR 00228 *> =N+5: error return from CGGHRD 00229 *> =N+6: error return from CHGEQZ (other than failed 00230 *> iteration) 00231 *> =N+7: error return from CTGEVC 00232 *> =N+8: error return from CGGBAK (computing VL) 00233 *> =N+9: error return from CGGBAK (computing VR) 00234 *> =N+10: error return from CLASCL (various calls) 00235 *> \endverbatim 00236 * 00237 * Authors: 00238 * ======== 00239 * 00240 *> \author Univ. of Tennessee 00241 *> \author Univ. of California Berkeley 00242 *> \author Univ. of Colorado Denver 00243 *> \author NAG Ltd. 00244 * 00245 *> \date November 2011 00246 * 00247 *> \ingroup complexGEeigen 00248 * 00249 *> \par Further Details: 00250 * ===================== 00251 *> 00252 *> \verbatim 00253 *> 00254 *> Balancing 00255 *> --------- 00256 *> 00257 *> This driver calls CGGBAL to both permute and scale rows and columns 00258 *> of A and B. The permutations PL and PR are chosen so that PL*A*PR 00259 *> and PL*B*R will be upper triangular except for the diagonal blocks 00260 *> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as 00261 *> possible. The diagonal scaling matrices DL and DR are chosen so 00262 *> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to 00263 *> one (except for the elements that start out zero.) 00264 *> 00265 *> After the eigenvalues and eigenvectors of the balanced matrices 00266 *> have been computed, CGGBAK transforms the eigenvectors back to what 00267 *> they would have been (in perfect arithmetic) if they had not been 00268 *> balanced. 00269 *> 00270 *> Contents of A and B on Exit 00271 *> -------- -- - --- - -- ---- 00272 *> 00273 *> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or 00274 *> both), then on exit the arrays A and B will contain the complex Schur 00275 *> form[*] of the "balanced" versions of A and B. If no eigenvectors 00276 *> are computed, then only the diagonal blocks will be correct. 00277 *> 00278 *> [*] In other words, upper triangular form. 00279 *> \endverbatim 00280 *> 00281 * ===================================================================== 00282 SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, 00283 $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) 00284 * 00285 * -- LAPACK driver routine (version 3.4.0) -- 00286 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00287 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00288 * November 2011 00289 * 00290 * .. Scalar Arguments .. 00291 CHARACTER JOBVL, JOBVR 00292 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00293 * .. 00294 * .. Array Arguments .. 00295 REAL RWORK( * ) 00296 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00297 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), 00298 $ WORK( * ) 00299 * .. 00300 * 00301 * ===================================================================== 00302 * 00303 * .. Parameters .. 00304 REAL ZERO, ONE 00305 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00306 COMPLEX CZERO, CONE 00307 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 00308 $ CONE = ( 1.0E0, 0.0E0 ) ) 00309 * .. 00310 * .. Local Scalars .. 00311 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY 00312 CHARACTER CHTEMP 00313 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO, 00314 $ IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR, 00315 $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3 00316 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM, 00317 $ BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI, 00318 $ SALFAR, SBETA, SCALE, TEMP 00319 COMPLEX X 00320 * .. 00321 * .. Local Arrays .. 00322 LOGICAL LDUMMA( 1 ) 00323 * .. 00324 * .. External Subroutines .. 00325 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY, 00326 $ CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, XERBLA 00327 * .. 00328 * .. External Functions .. 00329 LOGICAL LSAME 00330 INTEGER ILAENV 00331 REAL CLANGE, SLAMCH 00332 EXTERNAL ILAENV, LSAME, CLANGE, SLAMCH 00333 * .. 00334 * .. Intrinsic Functions .. 00335 INTRINSIC ABS, AIMAG, CMPLX, INT, MAX, REAL 00336 * .. 00337 * .. Statement Functions .. 00338 REAL ABS1 00339 * .. 00340 * .. Statement Function definitions .. 00341 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00342 * .. 00343 * .. Executable Statements .. 00344 * 00345 * Decode the input arguments 00346 * 00347 IF( LSAME( JOBVL, 'N' ) ) THEN 00348 IJOBVL = 1 00349 ILVL = .FALSE. 00350 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN 00351 IJOBVL = 2 00352 ILVL = .TRUE. 00353 ELSE 00354 IJOBVL = -1 00355 ILVL = .FALSE. 00356 END IF 00357 * 00358 IF( LSAME( JOBVR, 'N' ) ) THEN 00359 IJOBVR = 1 00360 ILVR = .FALSE. 00361 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN 00362 IJOBVR = 2 00363 ILVR = .TRUE. 00364 ELSE 00365 IJOBVR = -1 00366 ILVR = .FALSE. 00367 END IF 00368 ILV = ILVL .OR. ILVR 00369 * 00370 * Test the input arguments 00371 * 00372 LWKMIN = MAX( 2*N, 1 ) 00373 LWKOPT = LWKMIN 00374 WORK( 1 ) = LWKOPT 00375 LQUERY = ( LWORK.EQ.-1 ) 00376 INFO = 0 00377 IF( IJOBVL.LE.0 ) THEN 00378 INFO = -1 00379 ELSE IF( IJOBVR.LE.0 ) THEN 00380 INFO = -2 00381 ELSE IF( N.LT.0 ) THEN 00382 INFO = -3 00383 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00384 INFO = -5 00385 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00386 INFO = -7 00387 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN 00388 INFO = -11 00389 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN 00390 INFO = -13 00391 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00392 INFO = -15 00393 END IF 00394 * 00395 IF( INFO.EQ.0 ) THEN 00396 NB1 = ILAENV( 1, 'CGEQRF', ' ', N, N, -1, -1 ) 00397 NB2 = ILAENV( 1, 'CUNMQR', ' ', N, N, N, -1 ) 00398 NB3 = ILAENV( 1, 'CUNGQR', ' ', N, N, N, -1 ) 00399 NB = MAX( NB1, NB2, NB3 ) 00400 LOPT = MAX( 2*N, N*(NB+1) ) 00401 WORK( 1 ) = LOPT 00402 END IF 00403 * 00404 IF( INFO.NE.0 ) THEN 00405 CALL XERBLA( 'CGEGV ', -INFO ) 00406 RETURN 00407 ELSE IF( LQUERY ) THEN 00408 RETURN 00409 END IF 00410 * 00411 * Quick return if possible 00412 * 00413 IF( N.EQ.0 ) 00414 $ RETURN 00415 * 00416 * Get machine constants 00417 * 00418 EPS = SLAMCH( 'E' )*SLAMCH( 'B' ) 00419 SAFMIN = SLAMCH( 'S' ) 00420 SAFMIN = SAFMIN + SAFMIN 00421 SAFMAX = ONE / SAFMIN 00422 * 00423 * Scale A 00424 * 00425 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK ) 00426 ANRM1 = ANRM 00427 ANRM2 = ONE 00428 IF( ANRM.LT.ONE ) THEN 00429 IF( SAFMAX*ANRM.LT.ONE ) THEN 00430 ANRM1 = SAFMIN 00431 ANRM2 = SAFMAX*ANRM 00432 END IF 00433 END IF 00434 * 00435 IF( ANRM.GT.ZERO ) THEN 00436 CALL CLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO ) 00437 IF( IINFO.NE.0 ) THEN 00438 INFO = N + 10 00439 RETURN 00440 END IF 00441 END IF 00442 * 00443 * Scale B 00444 * 00445 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK ) 00446 BNRM1 = BNRM 00447 BNRM2 = ONE 00448 IF( BNRM.LT.ONE ) THEN 00449 IF( SAFMAX*BNRM.LT.ONE ) THEN 00450 BNRM1 = SAFMIN 00451 BNRM2 = SAFMAX*BNRM 00452 END IF 00453 END IF 00454 * 00455 IF( BNRM.GT.ZERO ) THEN 00456 CALL CLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO ) 00457 IF( IINFO.NE.0 ) THEN 00458 INFO = N + 10 00459 RETURN 00460 END IF 00461 END IF 00462 * 00463 * Permute the matrix to make it more nearly triangular 00464 * Also "balance" the matrix. 00465 * 00466 ILEFT = 1 00467 IRIGHT = N + 1 00468 IRWORK = IRIGHT + N 00469 CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), 00470 $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO ) 00471 IF( IINFO.NE.0 ) THEN 00472 INFO = N + 1 00473 GO TO 80 00474 END IF 00475 * 00476 * Reduce B to triangular form, and initialize VL and/or VR 00477 * 00478 IROWS = IHI + 1 - ILO 00479 IF( ILV ) THEN 00480 ICOLS = N + 1 - ILO 00481 ELSE 00482 ICOLS = IROWS 00483 END IF 00484 ITAU = 1 00485 IWORK = ITAU + IROWS 00486 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 00487 $ WORK( IWORK ), LWORK+1-IWORK, IINFO ) 00488 IF( IINFO.GE.0 ) 00489 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00490 IF( IINFO.NE.0 ) THEN 00491 INFO = N + 2 00492 GO TO 80 00493 END IF 00494 * 00495 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 00496 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ), 00497 $ LWORK+1-IWORK, IINFO ) 00498 IF( IINFO.GE.0 ) 00499 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00500 IF( IINFO.NE.0 ) THEN 00501 INFO = N + 3 00502 GO TO 80 00503 END IF 00504 * 00505 IF( ILVL ) THEN 00506 CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL ) 00507 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 00508 $ VL( ILO+1, ILO ), LDVL ) 00509 CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, 00510 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK, 00511 $ IINFO ) 00512 IF( IINFO.GE.0 ) 00513 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00514 IF( IINFO.NE.0 ) THEN 00515 INFO = N + 4 00516 GO TO 80 00517 END IF 00518 END IF 00519 * 00520 IF( ILVR ) 00521 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR ) 00522 * 00523 * Reduce to generalized Hessenberg form 00524 * 00525 IF( ILV ) THEN 00526 * 00527 * Eigenvectors requested -- work on whole matrix. 00528 * 00529 CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, 00530 $ LDVL, VR, LDVR, IINFO ) 00531 ELSE 00532 CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, 00533 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO ) 00534 END IF 00535 IF( IINFO.NE.0 ) THEN 00536 INFO = N + 5 00537 GO TO 80 00538 END IF 00539 * 00540 * Perform QZ algorithm 00541 * 00542 IWORK = ITAU 00543 IF( ILV ) THEN 00544 CHTEMP = 'S' 00545 ELSE 00546 CHTEMP = 'E' 00547 END IF 00548 CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, 00549 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ), 00550 $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO ) 00551 IF( IINFO.GE.0 ) 00552 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00553 IF( IINFO.NE.0 ) THEN 00554 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN 00555 INFO = IINFO 00556 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN 00557 INFO = IINFO - N 00558 ELSE 00559 INFO = N + 6 00560 END IF 00561 GO TO 80 00562 END IF 00563 * 00564 IF( ILV ) THEN 00565 * 00566 * Compute Eigenvectors 00567 * 00568 IF( ILVL ) THEN 00569 IF( ILVR ) THEN 00570 CHTEMP = 'B' 00571 ELSE 00572 CHTEMP = 'L' 00573 END IF 00574 ELSE 00575 CHTEMP = 'R' 00576 END IF 00577 * 00578 CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, 00579 $ VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ), 00580 $ IINFO ) 00581 IF( IINFO.NE.0 ) THEN 00582 INFO = N + 7 00583 GO TO 80 00584 END IF 00585 * 00586 * Undo balancing on VL and VR, rescale 00587 * 00588 IF( ILVL ) THEN 00589 CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), 00590 $ RWORK( IRIGHT ), N, VL, LDVL, IINFO ) 00591 IF( IINFO.NE.0 ) THEN 00592 INFO = N + 8 00593 GO TO 80 00594 END IF 00595 DO 30 JC = 1, N 00596 TEMP = ZERO 00597 DO 10 JR = 1, N 00598 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) ) 00599 10 CONTINUE 00600 IF( TEMP.LT.SAFMIN ) 00601 $ GO TO 30 00602 TEMP = ONE / TEMP 00603 DO 20 JR = 1, N 00604 VL( JR, JC ) = VL( JR, JC )*TEMP 00605 20 CONTINUE 00606 30 CONTINUE 00607 END IF 00608 IF( ILVR ) THEN 00609 CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), 00610 $ RWORK( IRIGHT ), N, VR, LDVR, IINFO ) 00611 IF( IINFO.NE.0 ) THEN 00612 INFO = N + 9 00613 GO TO 80 00614 END IF 00615 DO 60 JC = 1, N 00616 TEMP = ZERO 00617 DO 40 JR = 1, N 00618 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) ) 00619 40 CONTINUE 00620 IF( TEMP.LT.SAFMIN ) 00621 $ GO TO 60 00622 TEMP = ONE / TEMP 00623 DO 50 JR = 1, N 00624 VR( JR, JC ) = VR( JR, JC )*TEMP 00625 50 CONTINUE 00626 60 CONTINUE 00627 END IF 00628 * 00629 * End of eigenvector calculation 00630 * 00631 END IF 00632 * 00633 * Undo scaling in alpha, beta 00634 * 00635 * Note: this does not give the alpha and beta for the unscaled 00636 * problem. 00637 * 00638 * Un-scaling is limited to avoid underflow in alpha and beta 00639 * if they are significant. 00640 * 00641 DO 70 JC = 1, N 00642 ABSAR = ABS( REAL( ALPHA( JC ) ) ) 00643 ABSAI = ABS( AIMAG( ALPHA( JC ) ) ) 00644 ABSB = ABS( REAL( BETA( JC ) ) ) 00645 SALFAR = ANRM*REAL( ALPHA( JC ) ) 00646 SALFAI = ANRM*AIMAG( ALPHA( JC ) ) 00647 SBETA = BNRM*REAL( BETA( JC ) ) 00648 ILIMIT = .FALSE. 00649 SCALE = ONE 00650 * 00651 * Check for significant underflow in imaginary part of ALPHA 00652 * 00653 IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE. 00654 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN 00655 ILIMIT = .TRUE. 00656 SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI ) 00657 END IF 00658 * 00659 * Check for significant underflow in real part of ALPHA 00660 * 00661 IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE. 00662 $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN 00663 ILIMIT = .TRUE. 00664 SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) / 00665 $ MAX( SAFMIN, ANRM2*ABSAR ) ) 00666 END IF 00667 * 00668 * Check for significant underflow in BETA 00669 * 00670 IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE. 00671 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN 00672 ILIMIT = .TRUE. 00673 SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) / 00674 $ MAX( SAFMIN, BNRM2*ABSB ) ) 00675 END IF 00676 * 00677 * Check for possible overflow when limiting scaling 00678 * 00679 IF( ILIMIT ) THEN 00680 TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ), 00681 $ ABS( SBETA ) ) 00682 IF( TEMP.GT.ONE ) 00683 $ SCALE = SCALE / TEMP 00684 IF( SCALE.LT.ONE ) 00685 $ ILIMIT = .FALSE. 00686 END IF 00687 * 00688 * Recompute un-scaled ALPHA, BETA if necessary. 00689 * 00690 IF( ILIMIT ) THEN 00691 SALFAR = ( SCALE*REAL( ALPHA( JC ) ) )*ANRM 00692 SALFAI = ( SCALE*AIMAG( ALPHA( JC ) ) )*ANRM 00693 SBETA = ( SCALE*BETA( JC ) )*BNRM 00694 END IF 00695 ALPHA( JC ) = CMPLX( SALFAR, SALFAI ) 00696 BETA( JC ) = SBETA 00697 70 CONTINUE 00698 * 00699 80 CONTINUE 00700 WORK( 1 ) = LWKOPT 00701 * 00702 RETURN 00703 * 00704 * End of CGEGV 00705 * 00706 END