![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CGGEVX 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 CGGEVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggevx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggevx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggevx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, 00022 * ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI, 00023 * LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, 00024 * WORK, LWORK, RWORK, IWORK, BWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * CHARACTER BALANC, JOBVL, JOBVR, SENSE 00028 * INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00029 * REAL ABNRM, BBNRM 00030 * .. 00031 * .. Array Arguments .. 00032 * LOGICAL BWORK( * ) 00033 * INTEGER IWORK( * ) 00034 * REAL LSCALE( * ), RCONDE( * ), RCONDV( * ), 00035 * $ RSCALE( * ), RWORK( * ) 00036 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00037 * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), 00038 * $ WORK( * ) 00039 * .. 00040 * 00041 * 00042 *> \par Purpose: 00043 * ============= 00044 *> 00045 *> \verbatim 00046 *> 00047 *> CGGEVX computes for a pair of N-by-N complex nonsymmetric matrices 00048 *> (A,B) the generalized eigenvalues, and optionally, the left and/or 00049 *> right generalized eigenvectors. 00050 *> 00051 *> Optionally, it also computes a balancing transformation to improve 00052 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI, 00053 *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for 00054 *> the eigenvalues (RCONDE), and reciprocal condition numbers for the 00055 *> right eigenvectors (RCONDV). 00056 *> 00057 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar 00058 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is 00059 *> singular. It is usually represented as the pair (alpha,beta), as 00060 *> there is a reasonable interpretation for beta=0, and even for both 00061 *> being zero. 00062 *> 00063 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j) 00064 *> of (A,B) satisfies 00065 *> A * v(j) = lambda(j) * B * v(j) . 00066 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j) 00067 *> of (A,B) satisfies 00068 *> u(j)**H * A = lambda(j) * u(j)**H * B. 00069 *> where u(j)**H is the conjugate-transpose of u(j). 00070 *> 00071 *> \endverbatim 00072 * 00073 * Arguments: 00074 * ========== 00075 * 00076 *> \param[in] BALANC 00077 *> \verbatim 00078 *> BALANC is CHARACTER*1 00079 *> Specifies the balance option to be performed: 00080 *> = 'N': do not diagonally scale or permute; 00081 *> = 'P': permute only; 00082 *> = 'S': scale only; 00083 *> = 'B': both permute and scale. 00084 *> Computed reciprocal condition numbers will be for the 00085 *> matrices after permuting and/or balancing. Permuting does 00086 *> not change condition numbers (in exact arithmetic), but 00087 *> balancing does. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] JOBVL 00091 *> \verbatim 00092 *> JOBVL is CHARACTER*1 00093 *> = 'N': do not compute the left generalized eigenvectors; 00094 *> = 'V': compute the left generalized eigenvectors. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] JOBVR 00098 *> \verbatim 00099 *> JOBVR is CHARACTER*1 00100 *> = 'N': do not compute the right generalized eigenvectors; 00101 *> = 'V': compute the right generalized eigenvectors. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] SENSE 00105 *> \verbatim 00106 *> SENSE is CHARACTER*1 00107 *> Determines which reciprocal condition numbers are computed. 00108 *> = 'N': none are computed; 00109 *> = 'E': computed for eigenvalues only; 00110 *> = 'V': computed for eigenvectors only; 00111 *> = 'B': computed for eigenvalues and eigenvectors. 00112 *> \endverbatim 00113 *> 00114 *> \param[in] N 00115 *> \verbatim 00116 *> N is INTEGER 00117 *> The order of the matrices A, B, VL, and VR. N >= 0. 00118 *> \endverbatim 00119 *> 00120 *> \param[in,out] A 00121 *> \verbatim 00122 *> A is COMPLEX array, dimension (LDA, N) 00123 *> On entry, the matrix A in the pair (A,B). 00124 *> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V' 00125 *> or both, then A contains the first part of the complex Schur 00126 *> form of the "balanced" versions of the input A and B. 00127 *> \endverbatim 00128 *> 00129 *> \param[in] LDA 00130 *> \verbatim 00131 *> LDA is INTEGER 00132 *> The leading dimension of A. LDA >= max(1,N). 00133 *> \endverbatim 00134 *> 00135 *> \param[in,out] B 00136 *> \verbatim 00137 *> B is COMPLEX array, dimension (LDB, N) 00138 *> On entry, the matrix B in the pair (A,B). 00139 *> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V' 00140 *> or both, then B contains the second part of the complex 00141 *> Schur form of the "balanced" versions of the input A and B. 00142 *> \endverbatim 00143 *> 00144 *> \param[in] LDB 00145 *> \verbatim 00146 *> LDB is INTEGER 00147 *> The leading dimension of B. LDB >= max(1,N). 00148 *> \endverbatim 00149 *> 00150 *> \param[out] ALPHA 00151 *> \verbatim 00152 *> ALPHA is COMPLEX array, dimension (N) 00153 *> \endverbatim 00154 *> 00155 *> \param[out] BETA 00156 *> \verbatim 00157 *> BETA is COMPLEX array, dimension (N) 00158 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized 00159 *> eigenvalues. 00160 *> 00161 *> Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or 00162 *> underflow, and BETA(j) may even be zero. Thus, the user 00163 *> should avoid naively computing the ratio ALPHA/BETA. 00164 *> However, ALPHA will be always less than and usually 00165 *> comparable with norm(A) in magnitude, and BETA always less 00166 *> than and usually comparable with norm(B). 00167 *> \endverbatim 00168 *> 00169 *> \param[out] VL 00170 *> \verbatim 00171 *> VL is COMPLEX array, dimension (LDVL,N) 00172 *> If JOBVL = 'V', the left generalized eigenvectors u(j) are 00173 *> stored one after another in the columns of VL, in the same 00174 *> order as their eigenvalues. 00175 *> Each eigenvector will be scaled so the largest component 00176 *> will have abs(real part) + abs(imag. part) = 1. 00177 *> Not referenced if JOBVL = 'N'. 00178 *> \endverbatim 00179 *> 00180 *> \param[in] LDVL 00181 *> \verbatim 00182 *> LDVL is INTEGER 00183 *> The leading dimension of the matrix VL. LDVL >= 1, and 00184 *> if JOBVL = 'V', LDVL >= N. 00185 *> \endverbatim 00186 *> 00187 *> \param[out] VR 00188 *> \verbatim 00189 *> VR is COMPLEX array, dimension (LDVR,N) 00190 *> If JOBVR = 'V', the right generalized eigenvectors v(j) are 00191 *> stored one after another in the columns of VR, in the same 00192 *> order as their eigenvalues. 00193 *> Each eigenvector will be scaled so the largest component 00194 *> will have abs(real part) + abs(imag. part) = 1. 00195 *> Not referenced if JOBVR = 'N'. 00196 *> \endverbatim 00197 *> 00198 *> \param[in] LDVR 00199 *> \verbatim 00200 *> LDVR is INTEGER 00201 *> The leading dimension of the matrix VR. LDVR >= 1, and 00202 *> if JOBVR = 'V', LDVR >= N. 00203 *> \endverbatim 00204 *> 00205 *> \param[out] ILO 00206 *> \verbatim 00207 *> ILO is INTEGER 00208 *> \endverbatim 00209 *> 00210 *> \param[out] IHI 00211 *> \verbatim 00212 *> IHI is INTEGER 00213 *> ILO and IHI are integer values such that on exit 00214 *> A(i,j) = 0 and B(i,j) = 0 if i > j and 00215 *> j = 1,...,ILO-1 or i = IHI+1,...,N. 00216 *> If BALANC = 'N' or 'S', ILO = 1 and IHI = N. 00217 *> \endverbatim 00218 *> 00219 *> \param[out] LSCALE 00220 *> \verbatim 00221 *> LSCALE is REAL array, dimension (N) 00222 *> Details of the permutations and scaling factors applied 00223 *> to the left side of A and B. If PL(j) is the index of the 00224 *> row interchanged with row j, and DL(j) is the scaling 00225 *> factor applied to row j, then 00226 *> LSCALE(j) = PL(j) for j = 1,...,ILO-1 00227 *> = DL(j) for j = ILO,...,IHI 00228 *> = PL(j) for j = IHI+1,...,N. 00229 *> The order in which the interchanges are made is N to IHI+1, 00230 *> then 1 to ILO-1. 00231 *> \endverbatim 00232 *> 00233 *> \param[out] RSCALE 00234 *> \verbatim 00235 *> RSCALE is REAL array, dimension (N) 00236 *> Details of the permutations and scaling factors applied 00237 *> to the right side of A and B. If PR(j) is the index of the 00238 *> column interchanged with column j, and DR(j) is the scaling 00239 *> factor applied to column j, then 00240 *> RSCALE(j) = PR(j) for j = 1,...,ILO-1 00241 *> = DR(j) for j = ILO,...,IHI 00242 *> = PR(j) for j = IHI+1,...,N 00243 *> The order in which the interchanges are made is N to IHI+1, 00244 *> then 1 to ILO-1. 00245 *> \endverbatim 00246 *> 00247 *> \param[out] ABNRM 00248 *> \verbatim 00249 *> ABNRM is REAL 00250 *> The one-norm of the balanced matrix A. 00251 *> \endverbatim 00252 *> 00253 *> \param[out] BBNRM 00254 *> \verbatim 00255 *> BBNRM is REAL 00256 *> The one-norm of the balanced matrix B. 00257 *> \endverbatim 00258 *> 00259 *> \param[out] RCONDE 00260 *> \verbatim 00261 *> RCONDE is REAL array, dimension (N) 00262 *> If SENSE = 'E' or 'B', the reciprocal condition numbers of 00263 *> the eigenvalues, stored in consecutive elements of the array. 00264 *> If SENSE = 'N' or 'V', RCONDE is not referenced. 00265 *> \endverbatim 00266 *> 00267 *> \param[out] RCONDV 00268 *> \verbatim 00269 *> RCONDV is REAL array, dimension (N) 00270 *> If SENSE = 'V' or 'B', the estimated reciprocal condition 00271 *> numbers of the eigenvectors, stored in consecutive elements 00272 *> of the array. If the eigenvalues cannot be reordered to 00273 *> compute RCONDV(j), RCONDV(j) is set to 0; this can only occur 00274 *> when the true value would be very small anyway. 00275 *> If SENSE = 'N' or 'E', RCONDV is not referenced. 00276 *> \endverbatim 00277 *> 00278 *> \param[out] WORK 00279 *> \verbatim 00280 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00281 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00282 *> \endverbatim 00283 *> 00284 *> \param[in] LWORK 00285 *> \verbatim 00286 *> LWORK is INTEGER 00287 *> The dimension of the array WORK. LWORK >= max(1,2*N). 00288 *> If SENSE = 'E', LWORK >= max(1,4*N). 00289 *> If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N). 00290 *> 00291 *> If LWORK = -1, then a workspace query is assumed; the routine 00292 *> only calculates the optimal size of the WORK array, returns 00293 *> this value as the first entry of the WORK array, and no error 00294 *> message related to LWORK is issued by XERBLA. 00295 *> \endverbatim 00296 *> 00297 *> \param[out] RWORK 00298 *> \verbatim 00299 *> RWORK is REAL array, dimension (lrwork) 00300 *> lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B', 00301 *> and at least max(1,2*N) otherwise. 00302 *> Real workspace. 00303 *> \endverbatim 00304 *> 00305 *> \param[out] IWORK 00306 *> \verbatim 00307 *> IWORK is INTEGER array, dimension (N+2) 00308 *> If SENSE = 'E', IWORK is not referenced. 00309 *> \endverbatim 00310 *> 00311 *> \param[out] BWORK 00312 *> \verbatim 00313 *> BWORK is LOGICAL array, dimension (N) 00314 *> If SENSE = 'N', BWORK is not referenced. 00315 *> \endverbatim 00316 *> 00317 *> \param[out] INFO 00318 *> \verbatim 00319 *> INFO is INTEGER 00320 *> = 0: successful exit 00321 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00322 *> = 1,...,N: 00323 *> The QZ iteration failed. No eigenvectors have been 00324 *> calculated, but ALPHA(j) and BETA(j) should be correct 00325 *> for j=INFO+1,...,N. 00326 *> > N: =N+1: other than QZ iteration failed in CHGEQZ. 00327 *> =N+2: error return from CTGEVC. 00328 *> \endverbatim 00329 * 00330 * Authors: 00331 * ======== 00332 * 00333 *> \author Univ. of Tennessee 00334 *> \author Univ. of California Berkeley 00335 *> \author Univ. of Colorado Denver 00336 *> \author NAG Ltd. 00337 * 00338 *> \date April 2012 00339 * 00340 *> \ingroup complexGEeigen 00341 * 00342 *> \par Further Details: 00343 * ===================== 00344 *> 00345 *> \verbatim 00346 *> 00347 *> Balancing a matrix pair (A,B) includes, first, permuting rows and 00348 *> columns to isolate eigenvalues, second, applying diagonal similarity 00349 *> transformation to the rows and columns to make the rows and columns 00350 *> as close in norm as possible. The computed reciprocal condition 00351 *> numbers correspond to the balanced matrix. Permuting rows and columns 00352 *> will not change the condition numbers (in exact arithmetic) but 00353 *> diagonal scaling will. For further explanation of balancing, see 00354 *> section 4.11.1.2 of LAPACK Users' Guide. 00355 *> 00356 *> An approximate error bound on the chordal distance between the i-th 00357 *> computed generalized eigenvalue w and the corresponding exact 00358 *> eigenvalue lambda is 00359 *> 00360 *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I) 00361 *> 00362 *> An approximate error bound for the angle between the i-th computed 00363 *> eigenvector VL(i) or VR(i) is given by 00364 *> 00365 *> EPS * norm(ABNRM, BBNRM) / DIF(i). 00366 *> 00367 *> For further explanation of the reciprocal condition numbers RCONDE 00368 *> and RCONDV, see section 4.11 of LAPACK User's Guide. 00369 *> \endverbatim 00370 *> 00371 * ===================================================================== 00372 SUBROUTINE CGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, 00373 $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI, 00374 $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, 00375 $ WORK, LWORK, RWORK, IWORK, BWORK, INFO ) 00376 * 00377 * -- LAPACK driver routine (version 3.4.1) -- 00378 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00379 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00380 * April 2012 00381 * 00382 * .. Scalar Arguments .. 00383 CHARACTER BALANC, JOBVL, JOBVR, SENSE 00384 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00385 REAL ABNRM, BBNRM 00386 * .. 00387 * .. Array Arguments .. 00388 LOGICAL BWORK( * ) 00389 INTEGER IWORK( * ) 00390 REAL LSCALE( * ), RCONDE( * ), RCONDV( * ), 00391 $ RSCALE( * ), RWORK( * ) 00392 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00393 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), 00394 $ WORK( * ) 00395 * .. 00396 * 00397 * ===================================================================== 00398 * 00399 * .. Parameters .. 00400 REAL ZERO, ONE 00401 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00402 COMPLEX CZERO, CONE 00403 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00404 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00405 * .. 00406 * .. Local Scalars .. 00407 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL, 00408 $ WANTSB, WANTSE, WANTSN, WANTSV 00409 CHARACTER CHTEMP 00410 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS, 00411 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK 00412 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, 00413 $ SMLNUM, TEMP 00414 COMPLEX X 00415 * .. 00416 * .. Local Arrays .. 00417 LOGICAL LDUMMA( 1 ) 00418 * .. 00419 * .. External Subroutines .. 00420 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY, 00421 $ CLASCL, CLASET, CTGEVC, CTGSNA, CUNGQR, CUNMQR, 00422 $ SLABAD, SLASCL, XERBLA 00423 * .. 00424 * .. External Functions .. 00425 LOGICAL LSAME 00426 INTEGER ILAENV 00427 REAL CLANGE, SLAMCH 00428 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH 00429 * .. 00430 * .. Intrinsic Functions .. 00431 INTRINSIC ABS, AIMAG, MAX, REAL, SQRT 00432 * .. 00433 * .. Statement Functions .. 00434 REAL ABS1 00435 * .. 00436 * .. Statement Function definitions .. 00437 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00438 * .. 00439 * .. Executable Statements .. 00440 * 00441 * Decode the input arguments 00442 * 00443 IF( LSAME( JOBVL, 'N' ) ) THEN 00444 IJOBVL = 1 00445 ILVL = .FALSE. 00446 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN 00447 IJOBVL = 2 00448 ILVL = .TRUE. 00449 ELSE 00450 IJOBVL = -1 00451 ILVL = .FALSE. 00452 END IF 00453 * 00454 IF( LSAME( JOBVR, 'N' ) ) THEN 00455 IJOBVR = 1 00456 ILVR = .FALSE. 00457 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN 00458 IJOBVR = 2 00459 ILVR = .TRUE. 00460 ELSE 00461 IJOBVR = -1 00462 ILVR = .FALSE. 00463 END IF 00464 ILV = ILVL .OR. ILVR 00465 * 00466 NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' ) 00467 WANTSN = LSAME( SENSE, 'N' ) 00468 WANTSE = LSAME( SENSE, 'E' ) 00469 WANTSV = LSAME( SENSE, 'V' ) 00470 WANTSB = LSAME( SENSE, 'B' ) 00471 * 00472 * Test the input arguments 00473 * 00474 INFO = 0 00475 LQUERY = ( LWORK.EQ.-1 ) 00476 IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR. 00477 $ LSAME( BALANC, 'B' ) ) ) THEN 00478 INFO = -1 00479 ELSE IF( IJOBVL.LE.0 ) THEN 00480 INFO = -2 00481 ELSE IF( IJOBVR.LE.0 ) THEN 00482 INFO = -3 00483 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) ) 00484 $ THEN 00485 INFO = -4 00486 ELSE IF( N.LT.0 ) THEN 00487 INFO = -5 00488 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00489 INFO = -7 00490 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00491 INFO = -9 00492 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN 00493 INFO = -13 00494 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN 00495 INFO = -15 00496 END IF 00497 * 00498 * Compute workspace 00499 * (Note: Comments in the code beginning "Workspace:" describe the 00500 * minimal amount of workspace needed at that point in the code, 00501 * as well as the preferred amount for good performance. 00502 * NB refers to the optimal block size for the immediately 00503 * following subroutine, as returned by ILAENV. The workspace is 00504 * computed assuming ILO = 1 and IHI = N, the worst case.) 00505 * 00506 IF( INFO.EQ.0 ) THEN 00507 IF( N.EQ.0 ) THEN 00508 MINWRK = 1 00509 MAXWRK = 1 00510 ELSE 00511 MINWRK = 2*N 00512 IF( WANTSE ) THEN 00513 MINWRK = 4*N 00514 ELSE IF( WANTSV .OR. WANTSB ) THEN 00515 MINWRK = 2*N*( N + 1) 00516 END IF 00517 MAXWRK = MINWRK 00518 MAXWRK = MAX( MAXWRK, 00519 $ N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) ) 00520 MAXWRK = MAX( MAXWRK, 00521 $ N + N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) ) 00522 IF( ILVL ) THEN 00523 MAXWRK = MAX( MAXWRK, N + 00524 $ N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, 0 ) ) 00525 END IF 00526 END IF 00527 WORK( 1 ) = MAXWRK 00528 * 00529 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00530 INFO = -25 00531 END IF 00532 END IF 00533 * 00534 IF( INFO.NE.0 ) THEN 00535 CALL XERBLA( 'CGGEVX', -INFO ) 00536 RETURN 00537 ELSE IF( LQUERY ) THEN 00538 RETURN 00539 END IF 00540 * 00541 * Quick return if possible 00542 * 00543 IF( N.EQ.0 ) 00544 $ RETURN 00545 * 00546 * Get machine constants 00547 * 00548 EPS = SLAMCH( 'P' ) 00549 SMLNUM = SLAMCH( 'S' ) 00550 BIGNUM = ONE / SMLNUM 00551 CALL SLABAD( SMLNUM, BIGNUM ) 00552 SMLNUM = SQRT( SMLNUM ) / EPS 00553 BIGNUM = ONE / SMLNUM 00554 * 00555 * Scale A if max element outside range [SMLNUM,BIGNUM] 00556 * 00557 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK ) 00558 ILASCL = .FALSE. 00559 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00560 ANRMTO = SMLNUM 00561 ILASCL = .TRUE. 00562 ELSE IF( ANRM.GT.BIGNUM ) THEN 00563 ANRMTO = BIGNUM 00564 ILASCL = .TRUE. 00565 END IF 00566 IF( ILASCL ) 00567 $ CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 00568 * 00569 * Scale B if max element outside range [SMLNUM,BIGNUM] 00570 * 00571 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK ) 00572 ILBSCL = .FALSE. 00573 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00574 BNRMTO = SMLNUM 00575 ILBSCL = .TRUE. 00576 ELSE IF( BNRM.GT.BIGNUM ) THEN 00577 BNRMTO = BIGNUM 00578 ILBSCL = .TRUE. 00579 END IF 00580 IF( ILBSCL ) 00581 $ CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 00582 * 00583 * Permute and/or balance the matrix pair (A,B) 00584 * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise) 00585 * 00586 CALL CGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, 00587 $ RWORK, IERR ) 00588 * 00589 * Compute ABNRM and BBNRM 00590 * 00591 ABNRM = CLANGE( '1', N, N, A, LDA, RWORK( 1 ) ) 00592 IF( ILASCL ) THEN 00593 RWORK( 1 ) = ABNRM 00594 CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1, 00595 $ IERR ) 00596 ABNRM = RWORK( 1 ) 00597 END IF 00598 * 00599 BBNRM = CLANGE( '1', N, N, B, LDB, RWORK( 1 ) ) 00600 IF( ILBSCL ) THEN 00601 RWORK( 1 ) = BBNRM 00602 CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1, 00603 $ IERR ) 00604 BBNRM = RWORK( 1 ) 00605 END IF 00606 * 00607 * Reduce B to triangular form (QR decomposition of B) 00608 * (Complex Workspace: need N, prefer N*NB ) 00609 * 00610 IROWS = IHI + 1 - ILO 00611 IF( ILV .OR. .NOT.WANTSN ) THEN 00612 ICOLS = N + 1 - ILO 00613 ELSE 00614 ICOLS = IROWS 00615 END IF 00616 ITAU = 1 00617 IWRK = ITAU + IROWS 00618 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 00619 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 00620 * 00621 * Apply the unitary transformation to A 00622 * (Complex Workspace: need N, prefer N*NB) 00623 * 00624 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 00625 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 00626 $ LWORK+1-IWRK, IERR ) 00627 * 00628 * Initialize VL and/or VR 00629 * (Workspace: need N, prefer N*NB) 00630 * 00631 IF( ILVL ) THEN 00632 CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL ) 00633 IF( IROWS.GT.1 ) THEN 00634 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 00635 $ VL( ILO+1, ILO ), LDVL ) 00636 END IF 00637 CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, 00638 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 00639 END IF 00640 * 00641 IF( ILVR ) 00642 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR ) 00643 * 00644 * Reduce to generalized Hessenberg form 00645 * (Workspace: none needed) 00646 * 00647 IF( ILV .OR. .NOT.WANTSN ) THEN 00648 * 00649 * Eigenvectors requested -- work on whole matrix. 00650 * 00651 CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, 00652 $ LDVL, VR, LDVR, IERR ) 00653 ELSE 00654 CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, 00655 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) 00656 END IF 00657 * 00658 * Perform QZ algorithm (Compute eigenvalues, and optionally, the 00659 * Schur forms and Schur vectors) 00660 * (Complex Workspace: need N) 00661 * (Real Workspace: need N) 00662 * 00663 IWRK = ITAU 00664 IF( ILV .OR. .NOT.WANTSN ) THEN 00665 CHTEMP = 'S' 00666 ELSE 00667 CHTEMP = 'E' 00668 END IF 00669 * 00670 CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, 00671 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ), 00672 $ LWORK+1-IWRK, RWORK, IERR ) 00673 IF( IERR.NE.0 ) THEN 00674 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 00675 INFO = IERR 00676 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 00677 INFO = IERR - N 00678 ELSE 00679 INFO = N + 1 00680 END IF 00681 GO TO 90 00682 END IF 00683 * 00684 * Compute Eigenvectors and estimate condition numbers if desired 00685 * CTGEVC: (Complex Workspace: need 2*N ) 00686 * (Real Workspace: need 2*N ) 00687 * CTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B') 00688 * (Integer Workspace: need N+2 ) 00689 * 00690 IF( ILV .OR. .NOT.WANTSN ) THEN 00691 IF( ILV ) THEN 00692 IF( ILVL ) THEN 00693 IF( ILVR ) THEN 00694 CHTEMP = 'B' 00695 ELSE 00696 CHTEMP = 'L' 00697 END IF 00698 ELSE 00699 CHTEMP = 'R' 00700 END IF 00701 * 00702 CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, 00703 $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK, 00704 $ IERR ) 00705 IF( IERR.NE.0 ) THEN 00706 INFO = N + 2 00707 GO TO 90 00708 END IF 00709 END IF 00710 * 00711 IF( .NOT.WANTSN ) THEN 00712 * 00713 * compute eigenvectors (STGEVC) and estimate condition 00714 * numbers (STGSNA). Note that the definition of the condition 00715 * number is not invariant under transformation (u,v) to 00716 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized 00717 * Schur form (S,T), Q and Z are orthogonal matrices. In order 00718 * to avoid using extra 2*N*N workspace, we have to 00719 * re-calculate eigenvectors and estimate the condition numbers 00720 * one at a time. 00721 * 00722 DO 20 I = 1, N 00723 * 00724 DO 10 J = 1, N 00725 BWORK( J ) = .FALSE. 00726 10 CONTINUE 00727 BWORK( I ) = .TRUE. 00728 * 00729 IWRK = N + 1 00730 IWRK1 = IWRK + N 00731 * 00732 IF( WANTSE .OR. WANTSB ) THEN 00733 CALL CTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB, 00734 $ WORK( 1 ), N, WORK( IWRK ), N, 1, M, 00735 $ WORK( IWRK1 ), RWORK, IERR ) 00736 IF( IERR.NE.0 ) THEN 00737 INFO = N + 2 00738 GO TO 90 00739 END IF 00740 END IF 00741 * 00742 CALL CTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB, 00743 $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ), 00744 $ RCONDV( I ), 1, M, WORK( IWRK1 ), 00745 $ LWORK-IWRK1+1, IWORK, IERR ) 00746 * 00747 20 CONTINUE 00748 END IF 00749 END IF 00750 * 00751 * Undo balancing on VL and VR and normalization 00752 * (Workspace: none needed) 00753 * 00754 IF( ILVL ) THEN 00755 CALL CGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL, 00756 $ LDVL, IERR ) 00757 * 00758 DO 50 JC = 1, N 00759 TEMP = ZERO 00760 DO 30 JR = 1, N 00761 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) ) 00762 30 CONTINUE 00763 IF( TEMP.LT.SMLNUM ) 00764 $ GO TO 50 00765 TEMP = ONE / TEMP 00766 DO 40 JR = 1, N 00767 VL( JR, JC ) = VL( JR, JC )*TEMP 00768 40 CONTINUE 00769 50 CONTINUE 00770 END IF 00771 * 00772 IF( ILVR ) THEN 00773 CALL CGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR, 00774 $ LDVR, IERR ) 00775 DO 80 JC = 1, N 00776 TEMP = ZERO 00777 DO 60 JR = 1, N 00778 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) ) 00779 60 CONTINUE 00780 IF( TEMP.LT.SMLNUM ) 00781 $ GO TO 80 00782 TEMP = ONE / TEMP 00783 DO 70 JR = 1, N 00784 VR( JR, JC ) = VR( JR, JC )*TEMP 00785 70 CONTINUE 00786 80 CONTINUE 00787 END IF 00788 * 00789 * Undo scaling if necessary 00790 * 00791 90 CONTINUE 00792 * 00793 IF( ILASCL ) 00794 $ CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 00795 * 00796 IF( ILBSCL ) 00797 $ CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 00798 * 00799 WORK( 1 ) = MAXWRK 00800 RETURN 00801 * 00802 * End of CGGEVX 00803 * 00804 END