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