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