1 SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
2 $ LDVR, WORK, LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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)**T * A = lambda(j) * u(j)**T
29 * where u(j)**T denotes the 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 (MAX(1,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 *
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 *
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 ..
121 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
122 CHARACTER SIDE
123 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
124 $ 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 ..
133 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
134 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
135 $ XERBLA
136 * ..
137 * .. External Functions ..
138 LOGICAL LSAME
139 INTEGER IDAMAX, ILAENV
140 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
141 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
142 $ DNRM2
143 * ..
144 * .. Intrinsic Functions ..
145 INTRINSIC MAX, SQRT
146 * ..
147 * .. Executable Statements ..
148 *
149 * Test the input arguments
150 *
151 INFO = 0
152 LQUERY = ( LWORK.EQ.-1 )
153 WANTVL = LSAME( JOBVL, 'V' )
154 WANTVR = LSAME( JOBVR, 'V' )
155 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
156 INFO = -1
157 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
158 INFO = -2
159 ELSE IF( N.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162 INFO = -5
163 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
164 INFO = -9
165 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
166 INFO = -11
167 END IF
168 *
169 * Compute workspace
170 * (Note: Comments in the code beginning "Workspace:" describe the
171 * minimal amount of workspace needed at that point in the code,
172 * as well as the preferred amount for good performance.
173 * NB refers to the optimal block size for the immediately
174 * following subroutine, as returned by ILAENV.
175 * HSWORK refers to the workspace preferred by DHSEQR, as
176 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
177 * the worst case.)
178 *
179 IF( INFO.EQ.0 ) THEN
180 IF( N.EQ.0 ) THEN
181 MINWRK = 1
182 MAXWRK = 1
183 ELSE
184 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
185 IF( WANTVL ) THEN
186 MINWRK = 4*N
187 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
188 $ 'DORGHR', ' ', N, 1, N, -1 ) )
189 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
190 $ WORK, -1, INFO )
191 HSWORK = WORK( 1 )
192 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
193 MAXWRK = MAX( MAXWRK, 4*N )
194 ELSE IF( WANTVR ) THEN
195 MINWRK = 4*N
196 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
197 $ 'DORGHR', ' ', N, 1, N, -1 ) )
198 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
199 $ WORK, -1, INFO )
200 HSWORK = WORK( 1 )
201 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
202 MAXWRK = MAX( MAXWRK, 4*N )
203 ELSE
204 MINWRK = 3*N
205 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
206 $ WORK, -1, INFO )
207 HSWORK = WORK( 1 )
208 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
209 END IF
210 MAXWRK = MAX( MAXWRK, MINWRK )
211 END IF
212 WORK( 1 ) = MAXWRK
213 *
214 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
215 INFO = -13
216 END IF
217 END IF
218 *
219 IF( INFO.NE.0 ) THEN
220 CALL XERBLA( 'DGEEV ', -INFO )
221 RETURN
222 ELSE IF( LQUERY ) THEN
223 RETURN
224 END IF
225 *
226 * Quick return if possible
227 *
228 IF( N.EQ.0 )
229 $ RETURN
230 *
231 * Get machine constants
232 *
233 EPS = DLAMCH( 'P' )
234 SMLNUM = DLAMCH( 'S' )
235 BIGNUM = ONE / SMLNUM
236 CALL DLABAD( SMLNUM, BIGNUM )
237 SMLNUM = SQRT( SMLNUM ) / EPS
238 BIGNUM = ONE / SMLNUM
239 *
240 * Scale A if max element outside range [SMLNUM,BIGNUM]
241 *
242 ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
243 SCALEA = .FALSE.
244 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
245 SCALEA = .TRUE.
246 CSCALE = SMLNUM
247 ELSE IF( ANRM.GT.BIGNUM ) THEN
248 SCALEA = .TRUE.
249 CSCALE = BIGNUM
250 END IF
251 IF( SCALEA )
252 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
253 *
254 * Balance the matrix
255 * (Workspace: need N)
256 *
257 IBAL = 1
258 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
259 *
260 * Reduce to upper Hessenberg form
261 * (Workspace: need 3*N, prefer 2*N+N*NB)
262 *
263 ITAU = IBAL + N
264 IWRK = ITAU + N
265 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
266 $ LWORK-IWRK+1, IERR )
267 *
268 IF( WANTVL ) THEN
269 *
270 * Want left eigenvectors
271 * Copy Householder vectors to VL
272 *
273 SIDE = 'L'
274 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
275 *
276 * Generate orthogonal matrix in VL
277 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
278 *
279 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
280 $ LWORK-IWRK+1, IERR )
281 *
282 * Perform QR iteration, accumulating Schur vectors in VL
283 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
284 *
285 IWRK = ITAU
286 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
287 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
288 *
289 IF( WANTVR ) THEN
290 *
291 * Want left and right eigenvectors
292 * Copy Schur vectors to VR
293 *
294 SIDE = 'B'
295 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
296 END IF
297 *
298 ELSE IF( WANTVR ) THEN
299 *
300 * Want right eigenvectors
301 * Copy Householder vectors to VR
302 *
303 SIDE = 'R'
304 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
305 *
306 * Generate orthogonal matrix in VR
307 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
308 *
309 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
310 $ LWORK-IWRK+1, IERR )
311 *
312 * Perform QR iteration, accumulating Schur vectors in VR
313 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
314 *
315 IWRK = ITAU
316 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
317 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
318 *
319 ELSE
320 *
321 * Compute eigenvalues only
322 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
323 *
324 IWRK = ITAU
325 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
326 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
327 END IF
328 *
329 * If INFO > 0 from DHSEQR, then quit
330 *
331 IF( INFO.GT.0 )
332 $ GO TO 50
333 *
334 IF( WANTVL .OR. WANTVR ) THEN
335 *
336 * Compute left and/or right eigenvectors
337 * (Workspace: need 4*N)
338 *
339 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
340 $ N, NOUT, WORK( IWRK ), IERR )
341 END IF
342 *
343 IF( WANTVL ) THEN
344 *
345 * Undo balancing of left eigenvectors
346 * (Workspace: need N)
347 *
348 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
349 $ IERR )
350 *
351 * Normalize left eigenvectors and make largest component real
352 *
353 DO 20 I = 1, N
354 IF( WI( I ).EQ.ZERO ) THEN
355 SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
356 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
357 ELSE IF( WI( I ).GT.ZERO ) THEN
358 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
359 $ DNRM2( N, VL( 1, I+1 ), 1 ) )
360 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
361 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
362 DO 10 K = 1, N
363 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
364 10 CONTINUE
365 K = IDAMAX( N, WORK( IWRK ), 1 )
366 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
367 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
368 VL( K, I+1 ) = ZERO
369 END IF
370 20 CONTINUE
371 END IF
372 *
373 IF( WANTVR ) THEN
374 *
375 * Undo balancing of right eigenvectors
376 * (Workspace: need N)
377 *
378 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
379 $ IERR )
380 *
381 * Normalize right eigenvectors and make largest component real
382 *
383 DO 40 I = 1, N
384 IF( WI( I ).EQ.ZERO ) THEN
385 SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
386 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
387 ELSE IF( WI( I ).GT.ZERO ) THEN
388 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
389 $ DNRM2( N, VR( 1, I+1 ), 1 ) )
390 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
391 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
392 DO 30 K = 1, N
393 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
394 30 CONTINUE
395 K = IDAMAX( N, WORK( IWRK ), 1 )
396 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
397 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
398 VR( K, I+1 ) = ZERO
399 END IF
400 40 CONTINUE
401 END IF
402 *
403 * Undo scaling if necessary
404 *
405 50 CONTINUE
406 IF( SCALEA ) THEN
407 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
408 $ MAX( N-INFO, 1 ), IERR )
409 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
410 $ MAX( N-INFO, 1 ), IERR )
411 IF( INFO.GT.0 ) THEN
412 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
413 $ IERR )
414 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
415 $ IERR )
416 END IF
417 END IF
418 *
419 WORK( 1 ) = MAXWRK
420 RETURN
421 *
422 * End of DGEEV
423 *
424 END
2 $ LDVR, WORK, LWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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)**T * A = lambda(j) * u(j)**T
29 * where u(j)**T denotes the 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 (MAX(1,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 *
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 *
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 ..
121 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
122 CHARACTER SIDE
123 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
124 $ 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 ..
133 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
134 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
135 $ XERBLA
136 * ..
137 * .. External Functions ..
138 LOGICAL LSAME
139 INTEGER IDAMAX, ILAENV
140 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
141 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
142 $ DNRM2
143 * ..
144 * .. Intrinsic Functions ..
145 INTRINSIC MAX, SQRT
146 * ..
147 * .. Executable Statements ..
148 *
149 * Test the input arguments
150 *
151 INFO = 0
152 LQUERY = ( LWORK.EQ.-1 )
153 WANTVL = LSAME( JOBVL, 'V' )
154 WANTVR = LSAME( JOBVR, 'V' )
155 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
156 INFO = -1
157 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
158 INFO = -2
159 ELSE IF( N.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
162 INFO = -5
163 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
164 INFO = -9
165 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
166 INFO = -11
167 END IF
168 *
169 * Compute workspace
170 * (Note: Comments in the code beginning "Workspace:" describe the
171 * minimal amount of workspace needed at that point in the code,
172 * as well as the preferred amount for good performance.
173 * NB refers to the optimal block size for the immediately
174 * following subroutine, as returned by ILAENV.
175 * HSWORK refers to the workspace preferred by DHSEQR, as
176 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
177 * the worst case.)
178 *
179 IF( INFO.EQ.0 ) THEN
180 IF( N.EQ.0 ) THEN
181 MINWRK = 1
182 MAXWRK = 1
183 ELSE
184 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
185 IF( WANTVL ) THEN
186 MINWRK = 4*N
187 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
188 $ 'DORGHR', ' ', N, 1, N, -1 ) )
189 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
190 $ WORK, -1, INFO )
191 HSWORK = WORK( 1 )
192 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
193 MAXWRK = MAX( MAXWRK, 4*N )
194 ELSE IF( WANTVR ) THEN
195 MINWRK = 4*N
196 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
197 $ 'DORGHR', ' ', N, 1, N, -1 ) )
198 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
199 $ WORK, -1, INFO )
200 HSWORK = WORK( 1 )
201 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
202 MAXWRK = MAX( MAXWRK, 4*N )
203 ELSE
204 MINWRK = 3*N
205 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
206 $ WORK, -1, INFO )
207 HSWORK = WORK( 1 )
208 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
209 END IF
210 MAXWRK = MAX( MAXWRK, MINWRK )
211 END IF
212 WORK( 1 ) = MAXWRK
213 *
214 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
215 INFO = -13
216 END IF
217 END IF
218 *
219 IF( INFO.NE.0 ) THEN
220 CALL XERBLA( 'DGEEV ', -INFO )
221 RETURN
222 ELSE IF( LQUERY ) THEN
223 RETURN
224 END IF
225 *
226 * Quick return if possible
227 *
228 IF( N.EQ.0 )
229 $ RETURN
230 *
231 * Get machine constants
232 *
233 EPS = DLAMCH( 'P' )
234 SMLNUM = DLAMCH( 'S' )
235 BIGNUM = ONE / SMLNUM
236 CALL DLABAD( SMLNUM, BIGNUM )
237 SMLNUM = SQRT( SMLNUM ) / EPS
238 BIGNUM = ONE / SMLNUM
239 *
240 * Scale A if max element outside range [SMLNUM,BIGNUM]
241 *
242 ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
243 SCALEA = .FALSE.
244 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
245 SCALEA = .TRUE.
246 CSCALE = SMLNUM
247 ELSE IF( ANRM.GT.BIGNUM ) THEN
248 SCALEA = .TRUE.
249 CSCALE = BIGNUM
250 END IF
251 IF( SCALEA )
252 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
253 *
254 * Balance the matrix
255 * (Workspace: need N)
256 *
257 IBAL = 1
258 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
259 *
260 * Reduce to upper Hessenberg form
261 * (Workspace: need 3*N, prefer 2*N+N*NB)
262 *
263 ITAU = IBAL + N
264 IWRK = ITAU + N
265 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
266 $ LWORK-IWRK+1, IERR )
267 *
268 IF( WANTVL ) THEN
269 *
270 * Want left eigenvectors
271 * Copy Householder vectors to VL
272 *
273 SIDE = 'L'
274 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
275 *
276 * Generate orthogonal matrix in VL
277 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
278 *
279 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
280 $ LWORK-IWRK+1, IERR )
281 *
282 * Perform QR iteration, accumulating Schur vectors in VL
283 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
284 *
285 IWRK = ITAU
286 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
287 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
288 *
289 IF( WANTVR ) THEN
290 *
291 * Want left and right eigenvectors
292 * Copy Schur vectors to VR
293 *
294 SIDE = 'B'
295 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
296 END IF
297 *
298 ELSE IF( WANTVR ) THEN
299 *
300 * Want right eigenvectors
301 * Copy Householder vectors to VR
302 *
303 SIDE = 'R'
304 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
305 *
306 * Generate orthogonal matrix in VR
307 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
308 *
309 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
310 $ LWORK-IWRK+1, IERR )
311 *
312 * Perform QR iteration, accumulating Schur vectors in VR
313 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
314 *
315 IWRK = ITAU
316 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
317 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
318 *
319 ELSE
320 *
321 * Compute eigenvalues only
322 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
323 *
324 IWRK = ITAU
325 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
326 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
327 END IF
328 *
329 * If INFO > 0 from DHSEQR, then quit
330 *
331 IF( INFO.GT.0 )
332 $ GO TO 50
333 *
334 IF( WANTVL .OR. WANTVR ) THEN
335 *
336 * Compute left and/or right eigenvectors
337 * (Workspace: need 4*N)
338 *
339 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
340 $ N, NOUT, WORK( IWRK ), IERR )
341 END IF
342 *
343 IF( WANTVL ) THEN
344 *
345 * Undo balancing of left eigenvectors
346 * (Workspace: need N)
347 *
348 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
349 $ IERR )
350 *
351 * Normalize left eigenvectors and make largest component real
352 *
353 DO 20 I = 1, N
354 IF( WI( I ).EQ.ZERO ) THEN
355 SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
356 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
357 ELSE IF( WI( I ).GT.ZERO ) THEN
358 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
359 $ DNRM2( N, VL( 1, I+1 ), 1 ) )
360 CALL DSCAL( N, SCL, VL( 1, I ), 1 )
361 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
362 DO 10 K = 1, N
363 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
364 10 CONTINUE
365 K = IDAMAX( N, WORK( IWRK ), 1 )
366 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
367 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
368 VL( K, I+1 ) = ZERO
369 END IF
370 20 CONTINUE
371 END IF
372 *
373 IF( WANTVR ) THEN
374 *
375 * Undo balancing of right eigenvectors
376 * (Workspace: need N)
377 *
378 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
379 $ IERR )
380 *
381 * Normalize right eigenvectors and make largest component real
382 *
383 DO 40 I = 1, N
384 IF( WI( I ).EQ.ZERO ) THEN
385 SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
386 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
387 ELSE IF( WI( I ).GT.ZERO ) THEN
388 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
389 $ DNRM2( N, VR( 1, I+1 ), 1 ) )
390 CALL DSCAL( N, SCL, VR( 1, I ), 1 )
391 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
392 DO 30 K = 1, N
393 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
394 30 CONTINUE
395 K = IDAMAX( N, WORK( IWRK ), 1 )
396 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
397 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
398 VR( K, I+1 ) = ZERO
399 END IF
400 40 CONTINUE
401 END IF
402 *
403 * Undo scaling if necessary
404 *
405 50 CONTINUE
406 IF( SCALEA ) THEN
407 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
408 $ MAX( N-INFO, 1 ), IERR )
409 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
410 $ MAX( N-INFO, 1 ), IERR )
411 IF( INFO.GT.0 ) THEN
412 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
413 $ IERR )
414 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
415 $ IERR )
416 END IF
417 END IF
418 *
419 WORK( 1 ) = MAXWRK
420 RETURN
421 *
422 * End of DGEEV
423 *
424 END