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