 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, 00002 $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBVL, JOBVR 00011 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00015 $ B( LDB, * ), BETA( * ), VL( LDVL, * ), 00016 $ VR( LDVR, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B) 00023 * the generalized eigenvalues, and optionally, the left and/or right 00024 * generalized eigenvectors. 00025 * 00026 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar 00027 * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is 00028 * singular. It is usually represented as the pair (alpha,beta), as 00029 * there is a reasonable interpretation for beta=0, and even for both 00030 * being zero. 00031 * 00032 * The right eigenvector v(j) corresponding to the eigenvalue lambda(j) 00033 * of (A,B) satisfies 00034 * 00035 * A * v(j) = lambda(j) * B * v(j). 00036 * 00037 * The left eigenvector u(j) corresponding to the eigenvalue lambda(j) 00038 * of (A,B) satisfies 00039 * 00040 * u(j)**H * A = lambda(j) * u(j)**H * B . 00041 * 00042 * where u(j)**H is the conjugate-transpose of u(j). 00043 * 00044 * 00045 * Arguments 00046 * ========= 00047 * 00048 * JOBVL (input) CHARACTER*1 00049 * = 'N': do not compute the left generalized eigenvectors; 00050 * = 'V': compute the left generalized eigenvectors. 00051 * 00052 * JOBVR (input) CHARACTER*1 00053 * = 'N': do not compute the right generalized eigenvectors; 00054 * = 'V': compute the right generalized eigenvectors. 00055 * 00056 * N (input) INTEGER 00057 * The order of the matrices A, B, VL, and VR. N >= 0. 00058 * 00059 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00060 * On entry, the matrix A in the pair (A,B). 00061 * On exit, A has been overwritten. 00062 * 00063 * LDA (input) INTEGER 00064 * The leading dimension of A. LDA >= max(1,N). 00065 * 00066 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00067 * On entry, the matrix B in the pair (A,B). 00068 * On exit, B has been overwritten. 00069 * 00070 * LDB (input) INTEGER 00071 * The leading dimension of B. LDB >= max(1,N). 00072 * 00073 * ALPHAR (output) DOUBLE PRECISION array, dimension (N) 00074 * ALPHAI (output) DOUBLE PRECISION array, dimension (N) 00075 * BETA (output) DOUBLE PRECISION array, dimension (N) 00076 * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will 00077 * be the generalized eigenvalues. If ALPHAI(j) is zero, then 00078 * the j-th eigenvalue is real; if positive, then the j-th and 00079 * (j+1)-st eigenvalues are a complex conjugate pair, with 00080 * ALPHAI(j+1) negative. 00081 * 00082 * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) 00083 * may easily over- or underflow, and BETA(j) may even be zero. 00084 * Thus, the user should avoid naively computing the ratio 00085 * alpha/beta. However, ALPHAR and ALPHAI will be always less 00086 * than and usually comparable with norm(A) in magnitude, and 00087 * BETA always less than and usually comparable with norm(B). 00088 * 00089 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) 00090 * If JOBVL = 'V', the left eigenvectors u(j) are stored one 00091 * after another in the columns of VL, in the same order as 00092 * their eigenvalues. If the j-th eigenvalue is real, then 00093 * u(j) = VL(:,j), the j-th column of VL. If the j-th and 00094 * (j+1)-th eigenvalues form a complex conjugate pair, then 00095 * u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1). 00096 * Each eigenvector is scaled so the largest component has 00097 * abs(real part)+abs(imag. part)=1. 00098 * Not referenced if JOBVL = 'N'. 00099 * 00100 * LDVL (input) INTEGER 00101 * The leading dimension of the matrix VL. LDVL >= 1, and 00102 * if JOBVL = 'V', LDVL >= N. 00103 * 00104 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) 00105 * If JOBVR = 'V', the right eigenvectors v(j) are stored one 00106 * after another in the columns of VR, in the same order as 00107 * their eigenvalues. If the j-th eigenvalue is real, then 00108 * v(j) = VR(:,j), the j-th column of VR. If the j-th and 00109 * (j+1)-th eigenvalues form a complex conjugate pair, then 00110 * v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1). 00111 * Each eigenvector is scaled so the largest component has 00112 * abs(real part)+abs(imag. part)=1. 00113 * Not referenced if JOBVR = 'N'. 00114 * 00115 * LDVR (input) INTEGER 00116 * The leading dimension of the matrix VR. LDVR >= 1, and 00117 * if JOBVR = 'V', LDVR >= N. 00118 * 00119 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00120 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00121 * 00122 * LWORK (input) INTEGER 00123 * The dimension of the array WORK. LWORK >= max(1,8*N). 00124 * For good performance, LWORK must generally be larger. 00125 * 00126 * If LWORK = -1, then a workspace query is assumed; the routine 00127 * only calculates the optimal size of the WORK array, returns 00128 * this value as the first entry of the WORK array, and no error 00129 * message related to LWORK is issued by XERBLA. 00130 * 00131 * INFO (output) INTEGER 00132 * = 0: successful exit 00133 * < 0: if INFO = -i, the i-th argument had an illegal value. 00134 * = 1,...,N: 00135 * The QZ iteration failed. No eigenvectors have been 00136 * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) 00137 * should be correct for j=INFO+1,...,N. 00138 * > N: =N+1: other than QZ iteration failed in DHGEQZ. 00139 * =N+2: error return from DTGEVC. 00140 * 00141 * ===================================================================== 00142 * 00143 * .. Parameters .. 00144 DOUBLE PRECISION ZERO, ONE 00145 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00146 * .. 00147 * .. Local Scalars .. 00148 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY 00149 CHARACTER CHTEMP 00150 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO, 00151 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK, 00152 $ MINWRK 00153 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, 00154 $ SMLNUM, TEMP 00155 * .. 00156 * .. Local Arrays .. 00157 LOGICAL LDUMMA( 1 ) 00158 * .. 00159 * .. External Subroutines .. 00160 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD, 00161 $ DLACPY,DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, 00162 $ XERBLA 00163 * .. 00164 * .. External Functions .. 00165 LOGICAL LSAME 00166 INTEGER ILAENV 00167 DOUBLE PRECISION DLAMCH, DLANGE 00168 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE 00169 * .. 00170 * .. Intrinsic Functions .. 00171 INTRINSIC ABS, MAX, SQRT 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Decode the input arguments 00176 * 00177 IF( LSAME( JOBVL, 'N' ) ) THEN 00178 IJOBVL = 1 00179 ILVL = .FALSE. 00180 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN 00181 IJOBVL = 2 00182 ILVL = .TRUE. 00183 ELSE 00184 IJOBVL = -1 00185 ILVL = .FALSE. 00186 END IF 00187 * 00188 IF( LSAME( JOBVR, 'N' ) ) THEN 00189 IJOBVR = 1 00190 ILVR = .FALSE. 00191 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN 00192 IJOBVR = 2 00193 ILVR = .TRUE. 00194 ELSE 00195 IJOBVR = -1 00196 ILVR = .FALSE. 00197 END IF 00198 ILV = ILVL .OR. ILVR 00199 * 00200 * Test the input arguments 00201 * 00202 INFO = 0 00203 LQUERY = ( LWORK.EQ.-1 ) 00204 IF( IJOBVL.LE.0 ) THEN 00205 INFO = -1 00206 ELSE IF( IJOBVR.LE.0 ) THEN 00207 INFO = -2 00208 ELSE IF( N.LT.0 ) THEN 00209 INFO = -3 00210 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00211 INFO = -5 00212 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00213 INFO = -7 00214 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN 00215 INFO = -12 00216 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN 00217 INFO = -14 00218 END IF 00219 * 00220 * Compute workspace 00221 * (Note: Comments in the code beginning "Workspace:" describe the 00222 * minimal amount of workspace needed at that point in the code, 00223 * as well as the preferred amount for good performance. 00224 * NB refers to the optimal block size for the immediately 00225 * following subroutine, as returned by ILAENV. The workspace is 00226 * computed assuming ILO = 1 and IHI = N, the worst case.) 00227 * 00228 IF( INFO.EQ.0 ) THEN 00229 MINWRK = MAX( 1, 8*N ) 00230 MAXWRK = MAX( 1, N*( 7 + 00231 $ ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) ) ) 00232 MAXWRK = MAX( MAXWRK, N*( 7 + 00233 $ ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) ) ) 00234 IF( ILVL ) THEN 00235 MAXWRK = MAX( MAXWRK, N*( 7 + 00236 $ ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) ) 00237 END IF 00238 WORK( 1 ) = MAXWRK 00239 * 00240 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 00241 $ INFO = -16 00242 END IF 00243 * 00244 IF( INFO.NE.0 ) THEN 00245 CALL XERBLA( 'DGGEV ', -INFO ) 00246 RETURN 00247 ELSE IF( LQUERY ) THEN 00248 RETURN 00249 END IF 00250 * 00251 * Quick return if possible 00252 * 00253 IF( N.EQ.0 ) 00254 $ RETURN 00255 * 00256 * Get machine constants 00257 * 00258 EPS = DLAMCH( 'P' ) 00259 SMLNUM = DLAMCH( 'S' ) 00260 BIGNUM = ONE / SMLNUM 00261 CALL DLABAD( SMLNUM, BIGNUM ) 00262 SMLNUM = SQRT( SMLNUM ) / EPS 00263 BIGNUM = ONE / SMLNUM 00264 * 00265 * Scale A if max element outside range [SMLNUM,BIGNUM] 00266 * 00267 ANRM = DLANGE( 'M', N, N, A, LDA, WORK ) 00268 ILASCL = .FALSE. 00269 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00270 ANRMTO = SMLNUM 00271 ILASCL = .TRUE. 00272 ELSE IF( ANRM.GT.BIGNUM ) THEN 00273 ANRMTO = BIGNUM 00274 ILASCL = .TRUE. 00275 END IF 00276 IF( ILASCL ) 00277 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 00278 * 00279 * Scale B if max element outside range [SMLNUM,BIGNUM] 00280 * 00281 BNRM = DLANGE( 'M', N, N, B, LDB, WORK ) 00282 ILBSCL = .FALSE. 00283 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00284 BNRMTO = SMLNUM 00285 ILBSCL = .TRUE. 00286 ELSE IF( BNRM.GT.BIGNUM ) THEN 00287 BNRMTO = BIGNUM 00288 ILBSCL = .TRUE. 00289 END IF 00290 IF( ILBSCL ) 00291 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 00292 * 00293 * Permute the matrices A, B to isolate eigenvalues if possible 00294 * (Workspace: need 6*N) 00295 * 00296 ILEFT = 1 00297 IRIGHT = N + 1 00298 IWRK = IRIGHT + N 00299 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ), 00300 $ WORK( IRIGHT ), WORK( IWRK ), IERR ) 00301 * 00302 * Reduce B to triangular form (QR decomposition of B) 00303 * (Workspace: need N, prefer N*NB) 00304 * 00305 IROWS = IHI + 1 - ILO 00306 IF( ILV ) THEN 00307 ICOLS = N + 1 - ILO 00308 ELSE 00309 ICOLS = IROWS 00310 END IF 00311 ITAU = IWRK 00312 IWRK = ITAU + IROWS 00313 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 00314 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 00315 * 00316 * Apply the orthogonal transformation to matrix A 00317 * (Workspace: need N, prefer N*NB) 00318 * 00319 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 00320 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 00321 $ LWORK+1-IWRK, IERR ) 00322 * 00323 * Initialize VL 00324 * (Workspace: need N, prefer N*NB) 00325 * 00326 IF( ILVL ) THEN 00327 CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL ) 00328 IF( IROWS.GT.1 ) THEN 00329 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 00330 $ VL( ILO+1, ILO ), LDVL ) 00331 END IF 00332 CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, 00333 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 00334 END IF 00335 * 00336 * Initialize VR 00337 * 00338 IF( ILVR ) 00339 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR ) 00340 * 00341 * Reduce to generalized Hessenberg form 00342 * (Workspace: none needed) 00343 * 00344 IF( ILV ) THEN 00345 * 00346 * Eigenvectors requested -- work on whole matrix. 00347 * 00348 CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, 00349 $ LDVL, VR, LDVR, IERR ) 00350 ELSE 00351 CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, 00352 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR ) 00353 END IF 00354 * 00355 * Perform QZ algorithm (Compute eigenvalues, and optionally, the 00356 * Schur forms and Schur vectors) 00357 * (Workspace: need N) 00358 * 00359 IWRK = ITAU 00360 IF( ILV ) THEN 00361 CHTEMP = 'S' 00362 ELSE 00363 CHTEMP = 'E' 00364 END IF 00365 CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, 00366 $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, 00367 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 00368 IF( IERR.NE.0 ) THEN 00369 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 00370 INFO = IERR 00371 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 00372 INFO = IERR - N 00373 ELSE 00374 INFO = N + 1 00375 END IF 00376 GO TO 110 00377 END IF 00378 * 00379 * Compute Eigenvectors 00380 * (Workspace: need 6*N) 00381 * 00382 IF( ILV ) THEN 00383 IF( ILVL ) THEN 00384 IF( ILVR ) THEN 00385 CHTEMP = 'B' 00386 ELSE 00387 CHTEMP = 'L' 00388 END IF 00389 ELSE 00390 CHTEMP = 'R' 00391 END IF 00392 CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, 00393 $ VR, LDVR, N, IN, WORK( IWRK ), IERR ) 00394 IF( IERR.NE.0 ) THEN 00395 INFO = N + 2 00396 GO TO 110 00397 END IF 00398 * 00399 * Undo balancing on VL and VR and normalization 00400 * (Workspace: none needed) 00401 * 00402 IF( ILVL ) THEN 00403 CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ), 00404 $ WORK( IRIGHT ), N, VL, LDVL, IERR ) 00405 DO 50 JC = 1, N 00406 IF( ALPHAI( JC ).LT.ZERO ) 00407 $ GO TO 50 00408 TEMP = ZERO 00409 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00410 DO 10 JR = 1, N 00411 TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) ) 00412 10 CONTINUE 00413 ELSE 00414 DO 20 JR = 1, N 00415 TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+ 00416 $ ABS( VL( JR, JC+1 ) ) ) 00417 20 CONTINUE 00418 END IF 00419 IF( TEMP.LT.SMLNUM ) 00420 $ GO TO 50 00421 TEMP = ONE / TEMP 00422 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00423 DO 30 JR = 1, N 00424 VL( JR, JC ) = VL( JR, JC )*TEMP 00425 30 CONTINUE 00426 ELSE 00427 DO 40 JR = 1, N 00428 VL( JR, JC ) = VL( JR, JC )*TEMP 00429 VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP 00430 40 CONTINUE 00431 END IF 00432 50 CONTINUE 00433 END IF 00434 IF( ILVR ) THEN 00435 CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ), 00436 $ WORK( IRIGHT ), N, VR, LDVR, IERR ) 00437 DO 100 JC = 1, N 00438 IF( ALPHAI( JC ).LT.ZERO ) 00439 $ GO TO 100 00440 TEMP = ZERO 00441 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00442 DO 60 JR = 1, N 00443 TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) ) 00444 60 CONTINUE 00445 ELSE 00446 DO 70 JR = 1, N 00447 TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+ 00448 $ ABS( VR( JR, JC+1 ) ) ) 00449 70 CONTINUE 00450 END IF 00451 IF( TEMP.LT.SMLNUM ) 00452 $ GO TO 100 00453 TEMP = ONE / TEMP 00454 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00455 DO 80 JR = 1, N 00456 VR( JR, JC ) = VR( JR, JC )*TEMP 00457 80 CONTINUE 00458 ELSE 00459 DO 90 JR = 1, N 00460 VR( JR, JC ) = VR( JR, JC )*TEMP 00461 VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP 00462 90 CONTINUE 00463 END IF 00464 100 CONTINUE 00465 END IF 00466 * 00467 * End of eigenvector calculation 00468 * 00469 END IF 00470 * 00471 * Undo scaling if necessary 00472 * 00473 IF( ILASCL ) THEN 00474 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR ) 00475 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR ) 00476 END IF 00477 * 00478 IF( ILBSCL ) THEN 00479 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 00480 END IF 00481 * 00482 110 CONTINUE 00483 * 00484 WORK( 1 ) = MAXWRK 00485 * 00486 RETURN 00487 * 00488 * End of DGGEV 00489 * 00490 END
1.7.3 
