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