Annotation of rpl/lapack/lapack/dgeev.f, revision 1.7
1.1 bertrand 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.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
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)**H * A = lambda(j) * u(j)**H
29: * where u(j)**H denotes the conjugate 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
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