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