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