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