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: *> \date June 2016
185: *
186: * @precisions fortran d -> s
187: *
188: *> \ingroup doubleGEeigen
189: *
190: * =====================================================================
191: SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
192: $ LDVR, WORK, LWORK, INFO )
193: implicit none
194: *
195: * -- LAPACK driver routine (version 3.6.1) --
196: * -- LAPACK is a software package provided by Univ. of Tennessee, --
197: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198: * June 2016
199: *
200: * .. Scalar Arguments ..
201: CHARACTER JOBVL, JOBVR
202: INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
203: * ..
204: * .. Array Arguments ..
205: DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
206: $ WI( * ), WORK( * ), WR( * )
207: * ..
208: *
209: * =====================================================================
210: *
211: * .. Parameters ..
212: DOUBLE PRECISION ZERO, ONE
213: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
214: * ..
215: * .. Local Scalars ..
216: LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
217: CHARACTER SIDE
218: INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
219: $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
220: DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
221: $ SN
222: * ..
223: * .. Local Arrays ..
224: LOGICAL SELECT( 1 )
225: DOUBLE PRECISION DUM( 1 )
226: * ..
227: * .. External Subroutines ..
228: EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
229: $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
230: $ XERBLA
231: * ..
232: * .. External Functions ..
233: LOGICAL LSAME
234: INTEGER IDAMAX, ILAENV
235: DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
236: EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
237: $ DNRM2
238: * ..
239: * .. Intrinsic Functions ..
240: INTRINSIC MAX, SQRT
241: * ..
242: * .. Executable Statements ..
243: *
244: * Test the input arguments
245: *
246: INFO = 0
247: LQUERY = ( LWORK.EQ.-1 )
248: WANTVL = LSAME( JOBVL, 'V' )
249: WANTVR = LSAME( JOBVR, 'V' )
250: IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
251: INFO = -1
252: ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
253: INFO = -2
254: ELSE IF( N.LT.0 ) THEN
255: INFO = -3
256: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
257: INFO = -5
258: ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
259: INFO = -9
260: ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
261: INFO = -11
262: END IF
263: *
264: * Compute workspace
265: * (Note: Comments in the code beginning "Workspace:" describe the
266: * minimal amount of workspace needed at that point in the code,
267: * as well as the preferred amount for good performance.
268: * NB refers to the optimal block size for the immediately
269: * following subroutine, as returned by ILAENV.
270: * HSWORK refers to the workspace preferred by DHSEQR, as
271: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
272: * the worst case.)
273: *
274: IF( INFO.EQ.0 ) THEN
275: IF( N.EQ.0 ) THEN
276: MINWRK = 1
277: MAXWRK = 1
278: ELSE
279: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
280: IF( WANTVL ) THEN
281: MINWRK = 4*N
282: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
283: $ 'DORGHR', ' ', N, 1, N, -1 ) )
284: CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
285: $ WORK, -1, INFO )
286: HSWORK = INT( WORK(1) )
287: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
288: CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
289: $ VL, LDVL, VR, LDVR, N, NOUT,
290: $ WORK, -1, IERR )
291: LWORK_TREVC = INT( WORK(1) )
292: MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
293: MAXWRK = MAX( MAXWRK, 4*N )
294: ELSE IF( WANTVR ) THEN
295: MINWRK = 4*N
296: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
297: $ 'DORGHR', ' ', N, 1, N, -1 ) )
298: CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
299: $ WORK, -1, INFO )
300: HSWORK = INT( WORK(1) )
301: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
302: CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
303: $ VL, LDVL, VR, LDVR, N, NOUT,
304: $ WORK, -1, IERR )
305: LWORK_TREVC = INT( WORK(1) )
306: MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
307: MAXWRK = MAX( MAXWRK, 4*N )
308: ELSE
309: MINWRK = 3*N
310: CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
311: $ WORK, -1, INFO )
312: HSWORK = INT( WORK(1) )
313: MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
314: END IF
315: MAXWRK = MAX( MAXWRK, MINWRK )
316: END IF
317: WORK( 1 ) = MAXWRK
318: *
319: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
320: INFO = -13
321: END IF
322: END IF
323: *
324: IF( INFO.NE.0 ) THEN
325: CALL XERBLA( 'DGEEV ', -INFO )
326: RETURN
327: ELSE IF( LQUERY ) THEN
328: RETURN
329: END IF
330: *
331: * Quick return if possible
332: *
333: IF( N.EQ.0 )
334: $ RETURN
335: *
336: * Get machine constants
337: *
338: EPS = DLAMCH( 'P' )
339: SMLNUM = DLAMCH( 'S' )
340: BIGNUM = ONE / SMLNUM
341: CALL DLABAD( SMLNUM, BIGNUM )
342: SMLNUM = SQRT( SMLNUM ) / EPS
343: BIGNUM = ONE / SMLNUM
344: *
345: * Scale A if max element outside range [SMLNUM,BIGNUM]
346: *
347: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
348: SCALEA = .FALSE.
349: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
350: SCALEA = .TRUE.
351: CSCALE = SMLNUM
352: ELSE IF( ANRM.GT.BIGNUM ) THEN
353: SCALEA = .TRUE.
354: CSCALE = BIGNUM
355: END IF
356: IF( SCALEA )
357: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
358: *
359: * Balance the matrix
360: * (Workspace: need N)
361: *
362: IBAL = 1
363: CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
364: *
365: * Reduce to upper Hessenberg form
366: * (Workspace: need 3*N, prefer 2*N+N*NB)
367: *
368: ITAU = IBAL + N
369: IWRK = ITAU + N
370: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
371: $ LWORK-IWRK+1, IERR )
372: *
373: IF( WANTVL ) THEN
374: *
375: * Want left eigenvectors
376: * Copy Householder vectors to VL
377: *
378: SIDE = 'L'
379: CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
380: *
381: * Generate orthogonal matrix in VL
382: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
383: *
384: CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
385: $ LWORK-IWRK+1, IERR )
386: *
387: * Perform QR iteration, accumulating Schur vectors in VL
388: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
389: *
390: IWRK = ITAU
391: CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
392: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
393: *
394: IF( WANTVR ) THEN
395: *
396: * Want left and right eigenvectors
397: * Copy Schur vectors to VR
398: *
399: SIDE = 'B'
400: CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
401: END IF
402: *
403: ELSE IF( WANTVR ) THEN
404: *
405: * Want right eigenvectors
406: * Copy Householder vectors to VR
407: *
408: SIDE = 'R'
409: CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
410: *
411: * Generate orthogonal matrix in VR
412: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
413: *
414: CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
415: $ LWORK-IWRK+1, IERR )
416: *
417: * Perform QR iteration, accumulating Schur vectors in VR
418: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
419: *
420: IWRK = ITAU
421: CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
422: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
423: *
424: ELSE
425: *
426: * Compute eigenvalues only
427: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
428: *
429: IWRK = ITAU
430: CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
431: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
432: END IF
433: *
434: * If INFO .NE. 0 from DHSEQR, then quit
435: *
436: IF( INFO.NE.0 )
437: $ GO TO 50
438: *
439: IF( WANTVL .OR. WANTVR ) THEN
440: *
441: * Compute left and/or right eigenvectors
442: * (Workspace: need 4*N, prefer N + N + 2*N*NB)
443: *
444: CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
445: $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
446: END IF
447: *
448: IF( WANTVL ) THEN
449: *
450: * Undo balancing of left eigenvectors
451: * (Workspace: need N)
452: *
453: CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
454: $ IERR )
455: *
456: * Normalize left eigenvectors and make largest component real
457: *
458: DO 20 I = 1, N
459: IF( WI( I ).EQ.ZERO ) THEN
460: SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
461: CALL DSCAL( N, SCL, VL( 1, I ), 1 )
462: ELSE IF( WI( I ).GT.ZERO ) THEN
463: SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
464: $ DNRM2( N, VL( 1, I+1 ), 1 ) )
465: CALL DSCAL( N, SCL, VL( 1, I ), 1 )
466: CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
467: DO 10 K = 1, N
468: WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
469: 10 CONTINUE
470: K = IDAMAX( N, WORK( IWRK ), 1 )
471: CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
472: CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
473: VL( K, I+1 ) = ZERO
474: END IF
475: 20 CONTINUE
476: END IF
477: *
478: IF( WANTVR ) THEN
479: *
480: * Undo balancing of right eigenvectors
481: * (Workspace: need N)
482: *
483: CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
484: $ IERR )
485: *
486: * Normalize right eigenvectors and make largest component real
487: *
488: DO 40 I = 1, N
489: IF( WI( I ).EQ.ZERO ) THEN
490: SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
491: CALL DSCAL( N, SCL, VR( 1, I ), 1 )
492: ELSE IF( WI( I ).GT.ZERO ) THEN
493: SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
494: $ DNRM2( N, VR( 1, I+1 ), 1 ) )
495: CALL DSCAL( N, SCL, VR( 1, I ), 1 )
496: CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
497: DO 30 K = 1, N
498: WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
499: 30 CONTINUE
500: K = IDAMAX( N, WORK( IWRK ), 1 )
501: CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
502: CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
503: VR( K, I+1 ) = ZERO
504: END IF
505: 40 CONTINUE
506: END IF
507: *
508: * Undo scaling if necessary
509: *
510: 50 CONTINUE
511: IF( SCALEA ) THEN
512: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
513: $ MAX( N-INFO, 1 ), IERR )
514: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
515: $ MAX( N-INFO, 1 ), IERR )
516: IF( INFO.GT.0 ) THEN
517: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
518: $ IERR )
519: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
520: $ IERR )
521: END IF
522: END IF
523: *
524: WORK( 1 ) = MAXWRK
525: RETURN
526: *
527: * End of DGEEV
528: *
529: 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