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