|
2329
|
1 SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, |
|
|
2 $ LDVR, WORK, LWORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK driver routine (version 3.0) -- |
|
2329
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
7 * June 30, 1999 |
|
2329
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER JOBVL, JOBVR |
|
|
11 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), |
|
|
15 $ WI( * ), WORK( * ), WR( * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * DGEEV computes for an N-by-N real nonsymmetric matrix A, the |
|
|
22 * eigenvalues and, optionally, the left and/or right eigenvectors. |
|
|
23 * |
|
|
24 * The right eigenvector v(j) of A satisfies |
|
|
25 * A * v(j) = lambda(j) * v(j) |
|
|
26 * where lambda(j) is its eigenvalue. |
|
|
27 * The left eigenvector u(j) of A satisfies |
|
|
28 * u(j)**H * A = lambda(j) * u(j)**H |
|
|
29 * where u(j)**H denotes the conjugate transpose of u(j). |
|
|
30 * |
|
|
31 * The computed eigenvectors are normalized to have Euclidean norm |
|
|
32 * equal to 1 and largest component real. |
|
|
33 * |
|
|
34 * Arguments |
|
|
35 * ========= |
|
|
36 * |
|
|
37 * JOBVL (input) CHARACTER*1 |
|
|
38 * = 'N': left eigenvectors of A are not computed; |
|
|
39 * = 'V': left eigenvectors of A are computed. |
|
|
40 * |
|
|
41 * JOBVR (input) CHARACTER*1 |
|
|
42 * = 'N': right eigenvectors of A are not computed; |
|
|
43 * = 'V': right eigenvectors of A are computed. |
|
|
44 * |
|
|
45 * N (input) INTEGER |
|
|
46 * The order of the matrix A. N >= 0. |
|
|
47 * |
|
|
48 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) |
|
|
49 * On entry, the N-by-N matrix A. |
|
|
50 * On exit, A has been overwritten. |
|
|
51 * |
|
|
52 * LDA (input) INTEGER |
|
|
53 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
54 * |
|
|
55 * WR (output) DOUBLE PRECISION array, dimension (N) |
|
|
56 * WI (output) DOUBLE PRECISION array, dimension (N) |
|
|
57 * WR and WI contain the real and imaginary parts, |
|
|
58 * respectively, of the computed eigenvalues. Complex |
|
|
59 * conjugate pairs of eigenvalues appear consecutively |
|
|
60 * with the eigenvalue having the positive imaginary part |
|
|
61 * first. |
|
|
62 * |
|
|
63 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) |
|
|
64 * If JOBVL = 'V', the left eigenvectors u(j) are stored one |
|
|
65 * after another in the columns of VL, in the same order |
|
|
66 * as their eigenvalues. |
|
|
67 * If JOBVL = 'N', VL is not referenced. |
|
|
68 * If the j-th eigenvalue is real, then u(j) = VL(:,j), |
|
|
69 * the j-th column of VL. |
|
|
70 * If the j-th and (j+1)-st eigenvalues form a complex |
|
|
71 * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and |
|
|
72 * u(j+1) = VL(:,j) - i*VL(:,j+1). |
|
|
73 * |
|
|
74 * LDVL (input) INTEGER |
|
|
75 * The leading dimension of the array VL. LDVL >= 1; if |
|
|
76 * JOBVL = 'V', LDVL >= N. |
|
|
77 * |
|
|
78 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) |
|
|
79 * If JOBVR = 'V', the right eigenvectors v(j) are stored one |
|
|
80 * after another in the columns of VR, in the same order |
|
|
81 * as their eigenvalues. |
|
|
82 * If JOBVR = 'N', VR is not referenced. |
|
|
83 * If the j-th eigenvalue is real, then v(j) = VR(:,j), |
|
|
84 * the j-th column of VR. |
|
|
85 * If the j-th and (j+1)-st eigenvalues form a complex |
|
|
86 * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and |
|
|
87 * v(j+1) = VR(:,j) - i*VR(:,j+1). |
|
|
88 * |
|
|
89 * LDVR (input) INTEGER |
|
|
90 * The leading dimension of the array VR. LDVR >= 1; if |
|
|
91 * JOBVR = 'V', LDVR >= N. |
|
|
92 * |
|
|
93 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) |
|
|
94 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
|
95 * |
|
|
96 * LWORK (input) INTEGER |
|
|
97 * The dimension of the array WORK. LWORK >= max(1,3*N), and |
|
|
98 * if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good |
|
|
99 * performance, LWORK must generally be larger. |
|
|
100 * |
|
3333
|
101 * If LWORK = -1, then a workspace query is assumed; the routine |
|
|
102 * only calculates the optimal size of the WORK array, returns |
|
|
103 * this value as the first entry of the WORK array, and no error |
|
|
104 * message related to LWORK is issued by XERBLA. |
|
|
105 * |
|
2329
|
106 * INFO (output) INTEGER |
|
|
107 * = 0: successful exit |
|
|
108 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
109 * > 0: if INFO = i, the QR algorithm failed to compute all the |
|
|
110 * eigenvalues, and no eigenvectors have been computed; |
|
|
111 * elements i+1:N of WR and WI contain eigenvalues which |
|
|
112 * have converged. |
|
|
113 * |
|
|
114 * ===================================================================== |
|
|
115 * |
|
|
116 * .. Parameters .. |
|
|
117 DOUBLE PRECISION ZERO, ONE |
|
|
118 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) |
|
|
119 * .. |
|
|
120 * .. Local Scalars .. |
|
3333
|
121 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR |
|
2329
|
122 CHARACTER SIDE |
|
|
123 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K, |
|
|
124 $ MAXB, MAXWRK, MINWRK, NOUT |
|
|
125 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, |
|
|
126 $ SN |
|
|
127 * .. |
|
|
128 * .. Local Arrays .. |
|
|
129 LOGICAL SELECT( 1 ) |
|
|
130 DOUBLE PRECISION DUM( 1 ) |
|
|
131 * .. |
|
|
132 * .. External Subroutines .. |
|
3333
|
133 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY, DLARTG, |
|
|
134 $ DLASCL, DORGHR, DROT, DSCAL, DTREVC, XERBLA |
|
2329
|
135 * .. |
|
|
136 * .. External Functions .. |
|
|
137 LOGICAL LSAME |
|
|
138 INTEGER IDAMAX, ILAENV |
|
|
139 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2 |
|
|
140 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2, |
|
|
141 $ DNRM2 |
|
|
142 * .. |
|
|
143 * .. Intrinsic Functions .. |
|
|
144 INTRINSIC MAX, MIN, SQRT |
|
|
145 * .. |
|
|
146 * .. Executable Statements .. |
|
|
147 * |
|
|
148 * Test the input arguments |
|
|
149 * |
|
|
150 INFO = 0 |
|
3333
|
151 LQUERY = ( LWORK.EQ.-1 ) |
|
2329
|
152 WANTVL = LSAME( JOBVL, 'V' ) |
|
|
153 WANTVR = LSAME( JOBVR, 'V' ) |
|
|
154 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN |
|
|
155 INFO = -1 |
|
|
156 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN |
|
|
157 INFO = -2 |
|
|
158 ELSE IF( N.LT.0 ) THEN |
|
|
159 INFO = -3 |
|
|
160 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
161 INFO = -5 |
|
|
162 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN |
|
|
163 INFO = -9 |
|
|
164 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN |
|
|
165 INFO = -11 |
|
|
166 END IF |
|
|
167 * |
|
|
168 * Compute workspace |
|
|
169 * (Note: Comments in the code beginning "Workspace:" describe the |
|
|
170 * minimal amount of workspace needed at that point in the code, |
|
|
171 * as well as the preferred amount for good performance. |
|
|
172 * NB refers to the optimal block size for the immediately |
|
|
173 * following subroutine, as returned by ILAENV. |
|
|
174 * HSWORK refers to the workspace preferred by DHSEQR, as |
|
|
175 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, |
|
|
176 * the worst case.) |
|
|
177 * |
|
|
178 MINWRK = 1 |
|
|
179 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN |
|
|
180 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 ) |
|
|
181 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN |
|
|
182 MINWRK = MAX( 1, 3*N ) |
|
|
183 MAXB = MAX( ILAENV( 8, 'DHSEQR', 'EN', N, 1, N, -1 ), 2 ) |
|
|
184 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'EN', N, 1, |
|
|
185 $ N, -1 ) ) ) |
|
|
186 HSWORK = MAX( K*( K+2 ), 2*N ) |
|
|
187 MAXWRK = MAX( MAXWRK, N+1, N+HSWORK ) |
|
|
188 ELSE |
|
|
189 MINWRK = MAX( 1, 4*N ) |
|
|
190 MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* |
|
|
191 $ ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) ) |
|
|
192 MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 ) |
|
|
193 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1, |
|
|
194 $ N, -1 ) ) ) |
|
|
195 HSWORK = MAX( K*( K+2 ), 2*N ) |
|
|
196 MAXWRK = MAX( MAXWRK, N+1, N+HSWORK ) |
|
|
197 MAXWRK = MAX( MAXWRK, 4*N ) |
|
|
198 END IF |
|
|
199 WORK( 1 ) = MAXWRK |
|
|
200 END IF |
|
3333
|
201 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN |
|
2329
|
202 INFO = -13 |
|
|
203 END IF |
|
|
204 IF( INFO.NE.0 ) THEN |
|
|
205 CALL XERBLA( 'DGEEV ', -INFO ) |
|
|
206 RETURN |
|
3333
|
207 ELSE IF( LQUERY ) THEN |
|
|
208 RETURN |
|
2329
|
209 END IF |
|
|
210 * |
|
|
211 * Quick return if possible |
|
|
212 * |
|
|
213 IF( N.EQ.0 ) |
|
|
214 $ RETURN |
|
|
215 * |
|
|
216 * Get machine constants |
|
|
217 * |
|
|
218 EPS = DLAMCH( 'P' ) |
|
|
219 SMLNUM = DLAMCH( 'S' ) |
|
|
220 BIGNUM = ONE / SMLNUM |
|
|
221 CALL DLABAD( SMLNUM, BIGNUM ) |
|
|
222 SMLNUM = SQRT( SMLNUM ) / EPS |
|
|
223 BIGNUM = ONE / SMLNUM |
|
|
224 * |
|
|
225 * Scale A if max element outside range [SMLNUM,BIGNUM] |
|
|
226 * |
|
|
227 ANRM = DLANGE( 'M', N, N, A, LDA, DUM ) |
|
|
228 SCALEA = .FALSE. |
|
|
229 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN |
|
|
230 SCALEA = .TRUE. |
|
|
231 CSCALE = SMLNUM |
|
|
232 ELSE IF( ANRM.GT.BIGNUM ) THEN |
|
|
233 SCALEA = .TRUE. |
|
|
234 CSCALE = BIGNUM |
|
|
235 END IF |
|
|
236 IF( SCALEA ) |
|
|
237 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) |
|
|
238 * |
|
|
239 * Balance the matrix |
|
|
240 * (Workspace: need N) |
|
|
241 * |
|
|
242 IBAL = 1 |
|
|
243 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR ) |
|
|
244 * |
|
|
245 * Reduce to upper Hessenberg form |
|
|
246 * (Workspace: need 3*N, prefer 2*N+N*NB) |
|
|
247 * |
|
|
248 ITAU = IBAL + N |
|
|
249 IWRK = ITAU + N |
|
|
250 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), |
|
|
251 $ LWORK-IWRK+1, IERR ) |
|
|
252 * |
|
|
253 IF( WANTVL ) THEN |
|
|
254 * |
|
|
255 * Want left eigenvectors |
|
|
256 * Copy Householder vectors to VL |
|
|
257 * |
|
|
258 SIDE = 'L' |
|
|
259 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL ) |
|
|
260 * |
|
|
261 * Generate orthogonal matrix in VL |
|
|
262 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) |
|
|
263 * |
|
|
264 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), |
|
|
265 $ LWORK-IWRK+1, IERR ) |
|
|
266 * |
|
|
267 * Perform QR iteration, accumulating Schur vectors in VL |
|
|
268 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) |
|
|
269 * |
|
|
270 IWRK = ITAU |
|
|
271 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, |
|
|
272 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) |
|
|
273 * |
|
|
274 IF( WANTVR ) THEN |
|
|
275 * |
|
|
276 * Want left and right eigenvectors |
|
|
277 * Copy Schur vectors to VR |
|
|
278 * |
|
|
279 SIDE = 'B' |
|
|
280 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) |
|
|
281 END IF |
|
|
282 * |
|
|
283 ELSE IF( WANTVR ) THEN |
|
|
284 * |
|
|
285 * Want right eigenvectors |
|
|
286 * Copy Householder vectors to VR |
|
|
287 * |
|
|
288 SIDE = 'R' |
|
|
289 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR ) |
|
|
290 * |
|
|
291 * Generate orthogonal matrix in VR |
|
|
292 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) |
|
|
293 * |
|
|
294 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), |
|
|
295 $ LWORK-IWRK+1, IERR ) |
|
|
296 * |
|
|
297 * Perform QR iteration, accumulating Schur vectors in VR |
|
|
298 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) |
|
|
299 * |
|
|
300 IWRK = ITAU |
|
|
301 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, |
|
|
302 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) |
|
|
303 * |
|
|
304 ELSE |
|
|
305 * |
|
|
306 * Compute eigenvalues only |
|
|
307 * (Workspace: need N+1, prefer N+HSWORK (see comments) ) |
|
|
308 * |
|
|
309 IWRK = ITAU |
|
|
310 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, |
|
|
311 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) |
|
|
312 END IF |
|
|
313 * |
|
|
314 * If INFO > 0 from DHSEQR, then quit |
|
|
315 * |
|
|
316 IF( INFO.GT.0 ) |
|
|
317 $ GO TO 50 |
|
|
318 * |
|
|
319 IF( WANTVL .OR. WANTVR ) THEN |
|
|
320 * |
|
|
321 * Compute left and/or right eigenvectors |
|
|
322 * (Workspace: need 4*N) |
|
|
323 * |
|
|
324 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, |
|
|
325 $ N, NOUT, WORK( IWRK ), IERR ) |
|
|
326 END IF |
|
|
327 * |
|
|
328 IF( WANTVL ) THEN |
|
|
329 * |
|
|
330 * Undo balancing of left eigenvectors |
|
|
331 * (Workspace: need N) |
|
|
332 * |
|
|
333 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL, |
|
|
334 $ IERR ) |
|
|
335 * |
|
|
336 * Normalize left eigenvectors and make largest component real |
|
|
337 * |
|
|
338 DO 20 I = 1, N |
|
|
339 IF( WI( I ).EQ.ZERO ) THEN |
|
|
340 SCL = ONE / DNRM2( N, VL( 1, I ), 1 ) |
|
|
341 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) |
|
|
342 ELSE IF( WI( I ).GT.ZERO ) THEN |
|
|
343 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ), |
|
|
344 $ DNRM2( N, VL( 1, I+1 ), 1 ) ) |
|
|
345 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) |
|
|
346 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 ) |
|
|
347 DO 10 K = 1, N |
|
|
348 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2 |
|
|
349 10 CONTINUE |
|
|
350 K = IDAMAX( N, WORK( IWRK ), 1 ) |
|
|
351 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) |
|
|
352 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) |
|
|
353 VL( K, I+1 ) = ZERO |
|
|
354 END IF |
|
|
355 20 CONTINUE |
|
|
356 END IF |
|
|
357 * |
|
|
358 IF( WANTVR ) THEN |
|
|
359 * |
|
|
360 * Undo balancing of right eigenvectors |
|
|
361 * (Workspace: need N) |
|
|
362 * |
|
|
363 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR, |
|
|
364 $ IERR ) |
|
|
365 * |
|
|
366 * Normalize right eigenvectors and make largest component real |
|
|
367 * |
|
|
368 DO 40 I = 1, N |
|
|
369 IF( WI( I ).EQ.ZERO ) THEN |
|
|
370 SCL = ONE / DNRM2( N, VR( 1, I ), 1 ) |
|
|
371 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) |
|
|
372 ELSE IF( WI( I ).GT.ZERO ) THEN |
|
|
373 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ), |
|
|
374 $ DNRM2( N, VR( 1, I+1 ), 1 ) ) |
|
|
375 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) |
|
|
376 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 ) |
|
|
377 DO 30 K = 1, N |
|
|
378 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2 |
|
|
379 30 CONTINUE |
|
|
380 K = IDAMAX( N, WORK( IWRK ), 1 ) |
|
|
381 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) |
|
|
382 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) |
|
|
383 VR( K, I+1 ) = ZERO |
|
|
384 END IF |
|
|
385 40 CONTINUE |
|
|
386 END IF |
|
|
387 * |
|
|
388 * Undo scaling if necessary |
|
|
389 * |
|
|
390 50 CONTINUE |
|
|
391 IF( SCALEA ) THEN |
|
|
392 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), |
|
|
393 $ MAX( N-INFO, 1 ), IERR ) |
|
|
394 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), |
|
|
395 $ MAX( N-INFO, 1 ), IERR ) |
|
|
396 IF( INFO.GT.0 ) THEN |
|
|
397 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, |
|
|
398 $ IERR ) |
|
|
399 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, |
|
|
400 $ IERR ) |
|
|
401 END IF |
|
|
402 END IF |
|
|
403 * |
|
|
404 WORK( 1 ) = MAXWRK |
|
|
405 RETURN |
|
|
406 * |
|
|
407 * End of DGEEV |
|
|
408 * |
|
|
409 END |
