 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, 00002 $ LDVR, WORK, LWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBVL, JOBVR 00011 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 00015 $ WI( * ), WORK( * ), WR( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGEEV computes for an N-by-N real nonsymmetric matrix A, the 00022 * eigenvalues and, optionally, the left and/or right eigenvectors. 00023 * 00024 * The right eigenvector v(j) of A satisfies 00025 * A * v(j) = lambda(j) * v(j) 00026 * where lambda(j) is its eigenvalue. 00027 * The left eigenvector u(j) of A satisfies 00028 * u(j)**T * A = lambda(j) * u(j)**T 00029 * where u(j)**T denotes the transpose of u(j). 00030 * 00031 * The computed eigenvectors are normalized to have Euclidean norm 00032 * equal to 1 and largest component real. 00033 * 00034 * Arguments 00035 * ========= 00036 * 00037 * JOBVL (input) CHARACTER*1 00038 * = 'N': left eigenvectors of A are not computed; 00039 * = 'V': left eigenvectors of A are computed. 00040 * 00041 * JOBVR (input) CHARACTER*1 00042 * = 'N': right eigenvectors of A are not computed; 00043 * = 'V': right eigenvectors of A are computed. 00044 * 00045 * N (input) INTEGER 00046 * The order of the matrix A. N >= 0. 00047 * 00048 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00049 * On entry, the N-by-N matrix A. 00050 * On exit, A has been overwritten. 00051 * 00052 * LDA (input) INTEGER 00053 * The leading dimension of the array A. LDA >= max(1,N). 00054 * 00055 * WR (output) DOUBLE PRECISION array, dimension (N) 00056 * WI (output) DOUBLE PRECISION array, dimension (N) 00057 * WR and WI contain the real and imaginary parts, 00058 * respectively, of the computed eigenvalues. Complex 00059 * conjugate pairs of eigenvalues appear consecutively 00060 * with the eigenvalue having the positive imaginary part 00061 * first. 00062 * 00063 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) 00064 * If JOBVL = 'V', the left eigenvectors u(j) are stored one 00065 * after another in the columns of VL, in the same order 00066 * as their eigenvalues. 00067 * If JOBVL = 'N', VL is not referenced. 00068 * If the j-th eigenvalue is real, then u(j) = VL(:,j), 00069 * the j-th column of VL. 00070 * If the j-th and (j+1)-st eigenvalues form a complex 00071 * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and 00072 * u(j+1) = VL(:,j) - i*VL(:,j+1). 00073 * 00074 * LDVL (input) INTEGER 00075 * The leading dimension of the array VL. LDVL >= 1; if 00076 * JOBVL = 'V', LDVL >= N. 00077 * 00078 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) 00079 * If JOBVR = 'V', the right eigenvectors v(j) are stored one 00080 * after another in the columns of VR, in the same order 00081 * as their eigenvalues. 00082 * If JOBVR = 'N', VR is not referenced. 00083 * If the j-th eigenvalue is real, then v(j) = VR(:,j), 00084 * the j-th column of VR. 00085 * If the j-th and (j+1)-st eigenvalues form a complex 00086 * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and 00087 * v(j+1) = VR(:,j) - i*VR(:,j+1). 00088 * 00089 * LDVR (input) INTEGER 00090 * The leading dimension of the array VR. LDVR >= 1; if 00091 * JOBVR = 'V', LDVR >= N. 00092 * 00093 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00094 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00095 * 00096 * LWORK (input) INTEGER 00097 * The dimension of the array WORK. LWORK >= max(1,3*N), and 00098 * if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good 00099 * performance, LWORK must generally be larger. 00100 * 00101 * If LWORK = -1, then a workspace query is assumed; the routine 00102 * only calculates the optimal size of the WORK array, returns 00103 * this value as the first entry of the WORK array, and no error 00104 * message related to LWORK is issued by XERBLA. 00105 * 00106 * INFO (output) INTEGER 00107 * = 0: successful exit 00108 * < 0: if INFO = -i, the i-th argument had an illegal value. 00109 * > 0: if INFO = i, the QR algorithm failed to compute all the 00110 * eigenvalues, and no eigenvectors have been computed; 00111 * elements i+1:N of WR and WI contain eigenvalues which 00112 * have converged. 00113 * 00114 * ===================================================================== 00115 * 00116 * .. Parameters .. 00117 DOUBLE PRECISION ZERO, ONE 00118 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00119 * .. 00120 * .. Local Scalars .. 00121 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR 00122 CHARACTER SIDE 00123 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K, 00124 $ MAXWRK, MINWRK, NOUT 00125 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 00126 $ SN 00127 * .. 00128 * .. Local Arrays .. 00129 LOGICAL SELECT( 1 ) 00130 DOUBLE PRECISION DUM( 1 ) 00131 * .. 00132 * .. External Subroutines .. 00133 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY, 00134 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC, 00135 $ XERBLA 00136 * .. 00137 * .. External Functions .. 00138 LOGICAL LSAME 00139 INTEGER IDAMAX, ILAENV 00140 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2 00141 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2, 00142 $ DNRM2 00143 * .. 00144 * .. Intrinsic Functions .. 00145 INTRINSIC MAX, SQRT 00146 * .. 00147 * .. Executable Statements .. 00148 * 00149 * Test the input arguments 00150 * 00151 INFO = 0 00152 LQUERY = ( LWORK.EQ.-1 ) 00153 WANTVL = LSAME( JOBVL, 'V' ) 00154 WANTVR = LSAME( JOBVR, 'V' ) 00155 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 00156 INFO = -1 00157 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 00158 INFO = -2 00159 ELSE IF( N.LT.0 ) THEN 00160 INFO = -3 00161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00162 INFO = -5 00163 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 00164 INFO = -9 00165 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 00166 INFO = -11 00167 END IF 00168 * 00169 * Compute workspace 00170 * (Note: Comments in the code beginning "Workspace:" describe the 00171 * minimal amount of workspace needed at that point in the code, 00172 * as well as the preferred amount for good performance. 00173 * NB refers to the optimal block size for the immediately 00174 * following subroutine, as returned by ILAENV. 00175 * HSWORK refers to the workspace preferred by DHSEQR, as 00176 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 00177 * the worst case.) 00178 * 00179 IF( INFO.EQ.0 ) THEN 00180 IF( N.EQ.0 ) THEN 00181 MINWRK = 1 00182 MAXWRK = 1 00183 ELSE 00184 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 ) 00185 IF( WANTVL ) THEN 00186 MINWRK = 4*N 00187 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00188 $ 'DORGHR', ' ', N, 1, N, -1 ) ) 00189 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, 00190 $ WORK, -1, INFO ) 00191 HSWORK = WORK( 1 ) 00192 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 00193 MAXWRK = MAX( MAXWRK, 4*N ) 00194 ELSE IF( WANTVR ) THEN 00195 MINWRK = 4*N 00196 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00197 $ 'DORGHR', ' ', N, 1, N, -1 ) ) 00198 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, 00199 $ WORK, -1, INFO ) 00200 HSWORK = WORK( 1 ) 00201 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 00202 MAXWRK = MAX( MAXWRK, 4*N ) 00203 ELSE 00204 MINWRK = 3*N 00205 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR, 00206 $ WORK, -1, INFO ) 00207 HSWORK = WORK( 1 ) 00208 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 00209 END IF 00210 MAXWRK = MAX( MAXWRK, MINWRK ) 00211 END IF 00212 WORK( 1 ) = MAXWRK 00213 * 00214 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00215 INFO = -13 00216 END IF 00217 END IF 00218 * 00219 IF( INFO.NE.0 ) THEN 00220 CALL XERBLA( 'DGEEV ', -INFO ) 00221 RETURN 00222 ELSE IF( LQUERY ) THEN 00223 RETURN 00224 END IF 00225 * 00226 * Quick return if possible 00227 * 00228 IF( N.EQ.0 ) 00229 $ RETURN 00230 * 00231 * Get machine constants 00232 * 00233 EPS = DLAMCH( 'P' ) 00234 SMLNUM = DLAMCH( 'S' ) 00235 BIGNUM = ONE / SMLNUM 00236 CALL DLABAD( SMLNUM, BIGNUM ) 00237 SMLNUM = SQRT( SMLNUM ) / EPS 00238 BIGNUM = ONE / SMLNUM 00239 * 00240 * Scale A if max element outside range [SMLNUM,BIGNUM] 00241 * 00242 ANRM = DLANGE( 'M', N, N, A, LDA, DUM ) 00243 SCALEA = .FALSE. 00244 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00245 SCALEA = .TRUE. 00246 CSCALE = SMLNUM 00247 ELSE IF( ANRM.GT.BIGNUM ) THEN 00248 SCALEA = .TRUE. 00249 CSCALE = BIGNUM 00250 END IF 00251 IF( SCALEA ) 00252 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 00253 * 00254 * Balance the matrix 00255 * (Workspace: need N) 00256 * 00257 IBAL = 1 00258 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR ) 00259 * 00260 * Reduce to upper Hessenberg form 00261 * (Workspace: need 3*N, prefer 2*N+N*NB) 00262 * 00263 ITAU = IBAL + N 00264 IWRK = ITAU + N 00265 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 00266 $ LWORK-IWRK+1, IERR ) 00267 * 00268 IF( WANTVL ) THEN 00269 * 00270 * Want left eigenvectors 00271 * Copy Householder vectors to VL 00272 * 00273 SIDE = 'L' 00274 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL ) 00275 * 00276 * Generate orthogonal matrix in VL 00277 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 00278 * 00279 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 00280 $ LWORK-IWRK+1, IERR ) 00281 * 00282 * Perform QR iteration, accumulating Schur vectors in VL 00283 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) 00284 * 00285 IWRK = ITAU 00286 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 00287 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00288 * 00289 IF( WANTVR ) THEN 00290 * 00291 * Want left and right eigenvectors 00292 * Copy Schur vectors to VR 00293 * 00294 SIDE = 'B' 00295 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 00296 END IF 00297 * 00298 ELSE IF( WANTVR ) THEN 00299 * 00300 * Want right eigenvectors 00301 * Copy Householder vectors to VR 00302 * 00303 SIDE = 'R' 00304 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR ) 00305 * 00306 * Generate orthogonal matrix in VR 00307 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 00308 * 00309 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 00310 $ LWORK-IWRK+1, IERR ) 00311 * 00312 * Perform QR iteration, accumulating Schur vectors in VR 00313 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) 00314 * 00315 IWRK = ITAU 00316 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00317 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00318 * 00319 ELSE 00320 * 00321 * Compute eigenvalues only 00322 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) 00323 * 00324 IWRK = ITAU 00325 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00326 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00327 END IF 00328 * 00329 * If INFO > 0 from DHSEQR, then quit 00330 * 00331 IF( INFO.GT.0 ) 00332 $ GO TO 50 00333 * 00334 IF( WANTVL .OR. WANTVR ) THEN 00335 * 00336 * Compute left and/or right eigenvectors 00337 * (Workspace: need 4*N) 00338 * 00339 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00340 $ N, NOUT, WORK( IWRK ), IERR ) 00341 END IF 00342 * 00343 IF( WANTVL ) THEN 00344 * 00345 * Undo balancing of left eigenvectors 00346 * (Workspace: need N) 00347 * 00348 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL, 00349 $ IERR ) 00350 * 00351 * Normalize left eigenvectors and make largest component real 00352 * 00353 DO 20 I = 1, N 00354 IF( WI( I ).EQ.ZERO ) THEN 00355 SCL = ONE / DNRM2( N, VL( 1, I ), 1 ) 00356 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 00357 ELSE IF( WI( I ).GT.ZERO ) THEN 00358 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ), 00359 $ DNRM2( N, VL( 1, I+1 ), 1 ) ) 00360 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 00361 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 ) 00362 DO 10 K = 1, N 00363 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2 00364 10 CONTINUE 00365 K = IDAMAX( N, WORK( IWRK ), 1 ) 00366 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 00367 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 00368 VL( K, I+1 ) = ZERO 00369 END IF 00370 20 CONTINUE 00371 END IF 00372 * 00373 IF( WANTVR ) THEN 00374 * 00375 * Undo balancing of right eigenvectors 00376 * (Workspace: need N) 00377 * 00378 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR, 00379 $ IERR ) 00380 * 00381 * Normalize right eigenvectors and make largest component real 00382 * 00383 DO 40 I = 1, N 00384 IF( WI( I ).EQ.ZERO ) THEN 00385 SCL = ONE / DNRM2( N, VR( 1, I ), 1 ) 00386 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 00387 ELSE IF( WI( I ).GT.ZERO ) THEN 00388 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ), 00389 $ DNRM2( N, VR( 1, I+1 ), 1 ) ) 00390 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 00391 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 ) 00392 DO 30 K = 1, N 00393 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2 00394 30 CONTINUE 00395 K = IDAMAX( N, WORK( IWRK ), 1 ) 00396 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 00397 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 00398 VR( K, I+1 ) = ZERO 00399 END IF 00400 40 CONTINUE 00401 END IF 00402 * 00403 * Undo scaling if necessary 00404 * 00405 50 CONTINUE 00406 IF( SCALEA ) THEN 00407 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 00408 $ MAX( N-INFO, 1 ), IERR ) 00409 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 00410 $ MAX( N-INFO, 1 ), IERR ) 00411 IF( INFO.GT.0 ) THEN 00412 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 00413 $ IERR ) 00414 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 00415 $ IERR ) 00416 END IF 00417 END IF 00418 * 00419 WORK( 1 ) = MAXWRK 00420 RETURN 00421 * 00422 * End of DGEEV 00423 * 00424 END
1.7.3 
