1: *> \brief <b> DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGEEV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeev.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeev.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeev.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
22: * LDVR, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBVL, JOBVR
26: * INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
30: * $ WI( * ), WORK( * ), WR( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DGEEV computes for an N-by-N real nonsymmetric matrix A, the
40: *> eigenvalues and, optionally, the left and/or right eigenvectors.
41: *>
42: *> The right eigenvector v(j) of A satisfies
43: *> A * v(j) = lambda(j) * v(j)
44: *> where lambda(j) is its eigenvalue.
45: *> The left eigenvector u(j) of A satisfies
46: *> u(j)**H * A = lambda(j) * u(j)**H
47: *> where u(j)**H denotes the conjugate-transpose of u(j).
48: *>
49: *> The computed eigenvectors are normalized to have Euclidean norm
50: *> equal to 1 and largest component real.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] JOBVL
57: *> \verbatim
58: *> JOBVL is CHARACTER*1
59: *> = 'N': left eigenvectors of A are not computed;
60: *> = 'V': left eigenvectors of A are computed.
61: *> \endverbatim
62: *>
63: *> \param[in] JOBVR
64: *> \verbatim
65: *> JOBVR is CHARACTER*1
66: *> = 'N': right eigenvectors of A are not computed;
67: *> = 'V': right eigenvectors of A are computed.
68: *> \endverbatim
69: *>
70: *> \param[in] N
71: *> \verbatim
72: *> N is INTEGER
73: *> The order of the matrix A. N >= 0.
74: *> \endverbatim
75: *>
76: *> \param[in,out] A
77: *> \verbatim
78: *> A is DOUBLE PRECISION array, dimension (LDA,N)
79: *> On entry, the N-by-N matrix A.
80: *> On exit, A has been overwritten.
81: *> \endverbatim
82: *>
83: *> \param[in] LDA
84: *> \verbatim
85: *> LDA is INTEGER
86: *> The leading dimension of the array A. LDA >= max(1,N).
87: *> \endverbatim
88: *>
89: *> \param[out] WR
90: *> \verbatim
91: *> WR is DOUBLE PRECISION array, dimension (N)
92: *> \endverbatim
93: *>
94: *> \param[out] WI
95: *> \verbatim
96: *> WI is DOUBLE PRECISION array, dimension (N)
97: *> WR and WI contain the real and imaginary parts,
98: *> respectively, of the computed eigenvalues. Complex
99: *> conjugate pairs of eigenvalues appear consecutively
100: *> with the eigenvalue having the positive imaginary part
101: *> first.
102: *> \endverbatim
103: *>
104: *> \param[out] VL
105: *> \verbatim
106: *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
107: *> If JOBVL = 'V', the left eigenvectors u(j) are stored one
108: *> after another in the columns of VL, in the same order
109: *> as their eigenvalues.
110: *> If JOBVL = 'N', VL is not referenced.
111: *> If the j-th eigenvalue is real, then u(j) = VL(:,j),
112: *> the j-th column of VL.
113: *> If the j-th and (j+1)-st eigenvalues form a complex
114: *> conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
115: *> u(j+1) = VL(:,j) - i*VL(:,j+1).
116: *> \endverbatim
117: *>
118: *> \param[in] LDVL
119: *> \verbatim
120: *> LDVL is INTEGER
121: *> The leading dimension of the array VL. LDVL >= 1; if
122: *> JOBVL = 'V', LDVL >= N.
123: *> \endverbatim
124: *>
125: *> \param[out] VR
126: *> \verbatim
127: *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
128: *> If JOBVR = 'V', the right eigenvectors v(j) are stored one
129: *> after another in the columns of VR, in the same order
130: *> as their eigenvalues.
131: *> If JOBVR = 'N', VR is not referenced.
132: *> If the j-th eigenvalue is real, then v(j) = VR(:,j),
133: *> the j-th column of VR.
134: *> If the j-th and (j+1)-st eigenvalues form a complex
135: *> conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
136: *> v(j+1) = VR(:,j) - i*VR(:,j+1).
137: *> \endverbatim
138: *>
139: *> \param[in] LDVR
140: *> \verbatim
141: *> LDVR is INTEGER
142: *> The leading dimension of the array VR. LDVR >= 1; if
143: *> JOBVR = 'V', LDVR >= N.
144: *> \endverbatim
145: *>
146: *> \param[out] WORK
147: *> \verbatim
148: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
149: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
150: *> \endverbatim
151: *>
152: *> \param[in] LWORK
153: *> \verbatim
154: *> LWORK is INTEGER
155: *> The dimension of the array WORK. LWORK >= max(1,3*N), and
156: *> if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
157: *> performance, LWORK must generally be larger.
158: *>
159: *> If LWORK = -1, then a workspace query is assumed; the routine
160: *> only calculates the optimal size of the WORK array, returns
161: *> this value as the first entry of the WORK array, and no error
162: *> message related to LWORK is issued by XERBLA.
163: *> \endverbatim
164: *>
165: *> \param[out] INFO
166: *> \verbatim
167: *> INFO is INTEGER
168: *> = 0: successful exit
169: *> < 0: if INFO = -i, the i-th argument had an illegal value.
170: *> > 0: if INFO = i, the QR algorithm failed to compute all the
171: *> eigenvalues, and no eigenvectors have been computed;
172: *> elements i+1:N of WR and WI contain eigenvalues which
173: *> have converged.
174: *> \endverbatim
175: *
176: * Authors:
177: * ========
178: *
179: *> \author Univ. of Tennessee
180: *> \author Univ. of California Berkeley
181: *> \author Univ. of Colorado Denver
182: *> \author NAG Ltd.
183: *
184: *
185: * @precisions fortran d -> s
186: *
187: *> \ingroup doubleGEeigen
188: *
189: * =====================================================================
190: SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
191: $ LDVR, WORK, LWORK, INFO )
192: implicit none
193: *
194: * -- LAPACK driver routine --
195: * -- LAPACK is a software package provided by Univ. of Tennessee, --
196: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197: *
198: * .. Scalar Arguments ..
199: CHARACTER JOBVL, JOBVR
200: INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
201: * ..
202: * .. Array Arguments ..
203: DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
204: $ WI( * ), WORK( * ), WR( * )
205: * ..
206: *
207: * =====================================================================
208: *
209: * .. Parameters ..
210: DOUBLE PRECISION ZERO, ONE
211: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
212: * ..
213: * .. Local Scalars ..
214: LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
215: CHARACTER SIDE
216: INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
217: $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
218: DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
219: $ SN
220: * ..
221: * .. Local Arrays ..
222: LOGICAL SELECT( 1 )
223: DOUBLE PRECISION DUM( 1 )
224: * ..
225: * .. External Subroutines ..
226: EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
227: $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
228: $ XERBLA
229: * ..
230: * .. External Functions ..
231: LOGICAL LSAME
232: INTEGER IDAMAX, ILAENV
233: DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
234: EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
235: $ DNRM2
236: * ..
237: * .. Intrinsic Functions ..
238: INTRINSIC MAX, SQRT
239: * ..
240: * .. Executable Statements ..
241: *
242: * Test the input arguments
243: *
244: INFO = 0
245: LQUERY = ( LWORK.EQ.-1 )
246: WANTVL = LSAME( JOBVL, 'V' )
247: WANTVR = LSAME( JOBVR, 'V' )
248: IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
249: INFO = -1
250: ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
251: INFO = -2
252: ELSE IF( N.LT.0 ) THEN
253: INFO = -3
254: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
255: INFO = -5
256: ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
257: INFO = -9
258: ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
259: INFO = -11
260: END IF
261: *
262: * Compute workspace
263: * (Note: Comments in the code beginning "Workspace:" describe the
264: * minimal amount of workspace needed at that point in the code,
265: * as well as the preferred amount for good performance.
266: * NB refers to the optimal block size for the immediately
267: * following subroutine, as returned by ILAENV.
268: * HSWORK refers to the workspace preferred by DHSEQR, as
269: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
270: * the worst case.)
271: *
272: IF( INFO.EQ.0 ) THEN
273: IF( N.EQ.0 ) THEN
274: MINWRK = 1
275: MAXWRK = 1
276: ELSE
277: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
278: IF( WANTVL ) THEN
279: MINWRK = 4*N
280: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
281: $ 'DORGHR', ' ', N, 1, N, -1 ) )
282: CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
283: $ WORK, -1, INFO )
284: HSWORK = INT( WORK(1) )
285: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
286: CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
287: $ VL, LDVL, VR, LDVR, N, NOUT,
288: $ WORK, -1, IERR )
289: LWORK_TREVC = INT( WORK(1) )
290: MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
291: MAXWRK = MAX( MAXWRK, 4*N )
292: ELSE IF( WANTVR ) THEN
293: MINWRK = 4*N
294: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
295: $ 'DORGHR', ' ', N, 1, N, -1 ) )
296: CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
297: $ WORK, -1, INFO )
298: HSWORK = INT( WORK(1) )
299: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
300: CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
301: $ VL, LDVL, VR, LDVR, N, NOUT,
302: $ WORK, -1, IERR )
303: LWORK_TREVC = INT( WORK(1) )
304: MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
305: MAXWRK = MAX( MAXWRK, 4*N )
306: ELSE
307: MINWRK = 3*N
308: CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
309: $ WORK, -1, INFO )
310: HSWORK = INT( WORK(1) )
311: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
312: END IF
313: MAXWRK = MAX( MAXWRK, MINWRK )
314: END IF
315: WORK( 1 ) = MAXWRK
316: *
317: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
318: INFO = -13
319: END IF
320: END IF
321: *
322: IF( INFO.NE.0 ) THEN
323: CALL XERBLA( 'DGEEV ', -INFO )
324: RETURN
325: ELSE IF( LQUERY ) THEN
326: RETURN
327: END IF
328: *
329: * Quick return if possible
330: *
331: IF( N.EQ.0 )
332: $ RETURN
333: *
334: * Get machine constants
335: *
336: EPS = DLAMCH( 'P' )
337: SMLNUM = DLAMCH( 'S' )
338: BIGNUM = ONE / SMLNUM
339: CALL DLABAD( SMLNUM, BIGNUM )
340: SMLNUM = SQRT( SMLNUM ) / EPS
341: BIGNUM = ONE / SMLNUM
342: *
343: * Scale A if max element outside range [SMLNUM,BIGNUM]
344: *
345: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
346: SCALEA = .FALSE.
347: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
348: SCALEA = .TRUE.
349: CSCALE = SMLNUM
350: ELSE IF( ANRM.GT.BIGNUM ) THEN
351: SCALEA = .TRUE.
352: CSCALE = BIGNUM
353: END IF
354: IF( SCALEA )
355: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
356: *
357: * Balance the matrix
358: * (Workspace: need N)
359: *
360: IBAL = 1
361: CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
362: *
363: * Reduce to upper Hessenberg form
364: * (Workspace: need 3*N, prefer 2*N+N*NB)
365: *
366: ITAU = IBAL + N
367: IWRK = ITAU + N
368: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
369: $ LWORK-IWRK+1, IERR )
370: *
371: IF( WANTVL ) THEN
372: *
373: * Want left eigenvectors
374: * Copy Householder vectors to VL
375: *
376: SIDE = 'L'
377: CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
378: *
379: * Generate orthogonal matrix in VL
380: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381: *
382: CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
383: $ LWORK-IWRK+1, IERR )
384: *
385: * Perform QR iteration, accumulating Schur vectors in VL
386: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
387: *
388: IWRK = ITAU
389: CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
390: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
391: *
392: IF( WANTVR ) THEN
393: *
394: * Want left and right eigenvectors
395: * Copy Schur vectors to VR
396: *
397: SIDE = 'B'
398: CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
399: END IF
400: *
401: ELSE IF( WANTVR ) THEN
402: *
403: * Want right eigenvectors
404: * Copy Householder vectors to VR
405: *
406: SIDE = 'R'
407: CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
408: *
409: * Generate orthogonal matrix in VR
410: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
411: *
412: CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
413: $ LWORK-IWRK+1, IERR )
414: *
415: * Perform QR iteration, accumulating Schur vectors in VR
416: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
417: *
418: IWRK = ITAU
419: CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
420: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
421: *
422: ELSE
423: *
424: * Compute eigenvalues only
425: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
426: *
427: IWRK = ITAU
428: CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
429: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
430: END IF
431: *
432: * If INFO .NE. 0 from DHSEQR, then quit
433: *
434: IF( INFO.NE.0 )
435: $ GO TO 50
436: *
437: IF( WANTVL .OR. WANTVR ) THEN
438: *
439: * Compute left and/or right eigenvectors
440: * (Workspace: need 4*N, prefer N + N + 2*N*NB)
441: *
442: CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
443: $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
444: END IF
445: *
446: IF( WANTVL ) THEN
447: *
448: * Undo balancing of left eigenvectors
449: * (Workspace: need N)
450: *
451: CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
452: $ IERR )
453: *
454: * Normalize left eigenvectors and make largest component real
455: *
456: DO 20 I = 1, N
457: IF( WI( I ).EQ.ZERO ) THEN
458: SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
459: CALL DSCAL( N, SCL, VL( 1, I ), 1 )
460: ELSE IF( WI( I ).GT.ZERO ) THEN
461: SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
462: $ DNRM2( N, VL( 1, I+1 ), 1 ) )
463: CALL DSCAL( N, SCL, VL( 1, I ), 1 )
464: CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
465: DO 10 K = 1, N
466: WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
467: 10 CONTINUE
468: K = IDAMAX( N, WORK( IWRK ), 1 )
469: CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
470: CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
471: VL( K, I+1 ) = ZERO
472: END IF
473: 20 CONTINUE
474: END IF
475: *
476: IF( WANTVR ) THEN
477: *
478: * Undo balancing of right eigenvectors
479: * (Workspace: need N)
480: *
481: CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
482: $ IERR )
483: *
484: * Normalize right eigenvectors and make largest component real
485: *
486: DO 40 I = 1, N
487: IF( WI( I ).EQ.ZERO ) THEN
488: SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
489: CALL DSCAL( N, SCL, VR( 1, I ), 1 )
490: ELSE IF( WI( I ).GT.ZERO ) THEN
491: SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
492: $ DNRM2( N, VR( 1, I+1 ), 1 ) )
493: CALL DSCAL( N, SCL, VR( 1, I ), 1 )
494: CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
495: DO 30 K = 1, N
496: WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
497: 30 CONTINUE
498: K = IDAMAX( N, WORK( IWRK ), 1 )
499: CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
500: CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
501: VR( K, I+1 ) = ZERO
502: END IF
503: 40 CONTINUE
504: END IF
505: *
506: * Undo scaling if necessary
507: *
508: 50 CONTINUE
509: IF( SCALEA ) THEN
510: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
511: $ MAX( N-INFO, 1 ), IERR )
512: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
513: $ MAX( N-INFO, 1 ), IERR )
514: IF( INFO.GT.0 ) THEN
515: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
516: $ IERR )
517: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
518: $ IERR )
519: END IF
520: END IF
521: *
522: WORK( 1 ) = MAXWRK
523: RETURN
524: *
525: * End of DGEEV
526: *
527: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgeev.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack