Annotation of rpl/lapack/lapack/zgeev.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
2: $ WORK, LWORK, RWORK, 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 RWORK( * )
15: COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
16: $ W( * ), WORK( * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
23: * eigenvalues and, optionally, the left and/or right eigenvectors.
24: *
25: * The right eigenvector v(j) of A satisfies
26: * A * v(j) = lambda(j) * v(j)
27: * where lambda(j) is its eigenvalue.
28: * The left eigenvector u(j) of A satisfies
29: * u(j)**H * A = lambda(j) * u(j)**H
30: * where u(j)**H denotes the conjugate transpose of u(j).
31: *
32: * The computed eigenvectors are normalized to have Euclidean norm
33: * equal to 1 and largest component real.
34: *
35: * Arguments
36: * =========
37: *
38: * JOBVL (input) CHARACTER*1
39: * = 'N': left eigenvectors of A are not computed;
40: * = 'V': left eigenvectors of are computed.
41: *
42: * JOBVR (input) CHARACTER*1
43: * = 'N': right eigenvectors of A are not computed;
44: * = 'V': right eigenvectors of A are computed.
45: *
46: * N (input) INTEGER
47: * The order of the matrix A. N >= 0.
48: *
49: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
50: * On entry, the N-by-N matrix A.
51: * On exit, A has been overwritten.
52: *
53: * LDA (input) INTEGER
54: * The leading dimension of the array A. LDA >= max(1,N).
55: *
56: * W (output) COMPLEX*16 array, dimension (N)
57: * W contains the computed eigenvalues.
58: *
59: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
60: * If JOBVL = 'V', the left eigenvectors u(j) are stored one
61: * after another in the columns of VL, in the same order
62: * as their eigenvalues.
63: * If JOBVL = 'N', VL is not referenced.
64: * u(j) = VL(:,j), the j-th column of VL.
65: *
66: * LDVL (input) INTEGER
67: * The leading dimension of the array VL. LDVL >= 1; if
68: * JOBVL = 'V', LDVL >= N.
69: *
70: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
71: * If JOBVR = 'V', the right eigenvectors v(j) are stored one
72: * after another in the columns of VR, in the same order
73: * as their eigenvalues.
74: * If JOBVR = 'N', VR is not referenced.
75: * v(j) = VR(:,j), the j-th column of VR.
76: *
77: * LDVR (input) INTEGER
78: * The leading dimension of the array VR. LDVR >= 1; if
79: * JOBVR = 'V', LDVR >= N.
80: *
81: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
82: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
83: *
84: * LWORK (input) INTEGER
85: * The dimension of the array WORK. LWORK >= max(1,2*N).
86: * For good performance, LWORK must generally be larger.
87: *
88: * If LWORK = -1, then a workspace query is assumed; the routine
89: * only calculates the optimal size of the WORK array, returns
90: * this value as the first entry of the WORK array, and no error
91: * message related to LWORK is issued by XERBLA.
92: *
93: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
94: *
95: * INFO (output) INTEGER
96: * = 0: successful exit
97: * < 0: if INFO = -i, the i-th argument had an illegal value.
98: * > 0: if INFO = i, the QR algorithm failed to compute all the
99: * eigenvalues, and no eigenvectors have been computed;
100: * elements and i+1:N of W contain eigenvalues which have
101: * converged.
102: *
103: * =====================================================================
104: *
105: * .. Parameters ..
106: DOUBLE PRECISION ZERO, ONE
107: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
108: * ..
109: * .. Local Scalars ..
110: LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
111: CHARACTER SIDE
112: INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
113: $ IWRK, K, MAXWRK, MINWRK, NOUT
114: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
115: COMPLEX*16 TMP
116: * ..
117: * .. Local Arrays ..
118: LOGICAL SELECT( 1 )
119: DOUBLE PRECISION DUM( 1 )
120: * ..
121: * .. External Subroutines ..
122: EXTERNAL DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
123: $ ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
124: * ..
125: * .. External Functions ..
126: LOGICAL LSAME
127: INTEGER IDAMAX, ILAENV
128: DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE
129: EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
130: * ..
131: * .. Intrinsic Functions ..
132: INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
133: * ..
134: * .. Executable Statements ..
135: *
136: * Test the input arguments
137: *
138: INFO = 0
139: LQUERY = ( LWORK.EQ.-1 )
140: WANTVL = LSAME( JOBVL, 'V' )
141: WANTVR = LSAME( JOBVR, 'V' )
142: IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
143: INFO = -1
144: ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
145: INFO = -2
146: ELSE IF( N.LT.0 ) THEN
147: INFO = -3
148: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
149: INFO = -5
150: ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
151: INFO = -8
152: ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
153: INFO = -10
154: END IF
155: *
156: * Compute workspace
157: * (Note: Comments in the code beginning "Workspace:" describe the
158: * minimal amount of workspace needed at that point in the code,
159: * as well as the preferred amount for good performance.
160: * CWorkspace refers to complex workspace, and RWorkspace to real
161: * workspace. NB refers to the optimal block size for the
162: * immediately following subroutine, as returned by ILAENV.
163: * HSWORK refers to the workspace preferred by ZHSEQR, as
164: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
165: * the worst case.)
166: *
167: IF( INFO.EQ.0 ) THEN
168: IF( N.EQ.0 ) THEN
169: MINWRK = 1
170: MAXWRK = 1
171: ELSE
172: MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
173: MINWRK = 2*N
174: IF( WANTVL ) THEN
175: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
176: $ ' ', N, 1, N, -1 ) )
177: CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
178: $ WORK, -1, INFO )
179: ELSE IF( WANTVR ) THEN
180: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
181: $ ' ', N, 1, N, -1 ) )
182: CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
183: $ WORK, -1, INFO )
184: ELSE
185: CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
186: $ WORK, -1, INFO )
187: END IF
188: HSWORK = WORK( 1 )
189: MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
190: END IF
191: WORK( 1 ) = MAXWRK
192: *
193: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
194: INFO = -12
195: END IF
196: END IF
197: *
198: IF( INFO.NE.0 ) THEN
199: CALL XERBLA( 'ZGEEV ', -INFO )
200: RETURN
201: ELSE IF( LQUERY ) THEN
202: RETURN
203: END IF
204: *
205: * Quick return if possible
206: *
207: IF( N.EQ.0 )
208: $ RETURN
209: *
210: * Get machine constants
211: *
212: EPS = DLAMCH( 'P' )
213: SMLNUM = DLAMCH( 'S' )
214: BIGNUM = ONE / SMLNUM
215: CALL DLABAD( SMLNUM, BIGNUM )
216: SMLNUM = SQRT( SMLNUM ) / EPS
217: BIGNUM = ONE / SMLNUM
218: *
219: * Scale A if max element outside range [SMLNUM,BIGNUM]
220: *
221: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
222: SCALEA = .FALSE.
223: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
224: SCALEA = .TRUE.
225: CSCALE = SMLNUM
226: ELSE IF( ANRM.GT.BIGNUM ) THEN
227: SCALEA = .TRUE.
228: CSCALE = BIGNUM
229: END IF
230: IF( SCALEA )
231: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
232: *
233: * Balance the matrix
234: * (CWorkspace: none)
235: * (RWorkspace: need N)
236: *
237: IBAL = 1
238: CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
239: *
240: * Reduce to upper Hessenberg form
241: * (CWorkspace: need 2*N, prefer N+N*NB)
242: * (RWorkspace: none)
243: *
244: ITAU = 1
245: IWRK = ITAU + N
246: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
247: $ LWORK-IWRK+1, IERR )
248: *
249: IF( WANTVL ) THEN
250: *
251: * Want left eigenvectors
252: * Copy Householder vectors to VL
253: *
254: SIDE = 'L'
255: CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
256: *
257: * Generate unitary matrix in VL
258: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
259: * (RWorkspace: none)
260: *
261: CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
262: $ LWORK-IWRK+1, IERR )
263: *
264: * Perform QR iteration, accumulating Schur vectors in VL
265: * (CWorkspace: need 1, prefer HSWORK (see comments) )
266: * (RWorkspace: none)
267: *
268: IWRK = ITAU
269: CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
270: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
271: *
272: IF( WANTVR ) THEN
273: *
274: * Want left and right eigenvectors
275: * Copy Schur vectors to VR
276: *
277: SIDE = 'B'
278: CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
279: END IF
280: *
281: ELSE IF( WANTVR ) THEN
282: *
283: * Want right eigenvectors
284: * Copy Householder vectors to VR
285: *
286: SIDE = 'R'
287: CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
288: *
289: * Generate unitary matrix in VR
290: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
291: * (RWorkspace: none)
292: *
293: CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
294: $ LWORK-IWRK+1, IERR )
295: *
296: * Perform QR iteration, accumulating Schur vectors in VR
297: * (CWorkspace: need 1, prefer HSWORK (see comments) )
298: * (RWorkspace: none)
299: *
300: IWRK = ITAU
301: CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
302: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
303: *
304: ELSE
305: *
306: * Compute eigenvalues only
307: * (CWorkspace: need 1, prefer HSWORK (see comments) )
308: * (RWorkspace: none)
309: *
310: IWRK = ITAU
311: CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
312: $ WORK( IWRK ), LWORK-IWRK+1, INFO )
313: END IF
314: *
315: * If INFO > 0 from ZHSEQR, then quit
316: *
317: IF( INFO.GT.0 )
318: $ GO TO 50
319: *
320: IF( WANTVL .OR. WANTVR ) THEN
321: *
322: * Compute left and/or right eigenvectors
323: * (CWorkspace: need 2*N)
324: * (RWorkspace: need 2*N)
325: *
326: IRWORK = IBAL + N
327: CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
328: $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
329: END IF
330: *
331: IF( WANTVL ) THEN
332: *
333: * Undo balancing of left eigenvectors
334: * (CWorkspace: none)
335: * (RWorkspace: need N)
336: *
337: CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
338: $ IERR )
339: *
340: * Normalize left eigenvectors and make largest component real
341: *
342: DO 20 I = 1, N
343: SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
344: CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
345: DO 10 K = 1, N
346: RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
347: $ DIMAG( VL( K, I ) )**2
348: 10 CONTINUE
349: K = IDAMAX( N, RWORK( IRWORK ), 1 )
350: TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
351: CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
352: VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
353: 20 CONTINUE
354: END IF
355: *
356: IF( WANTVR ) THEN
357: *
358: * Undo balancing of right eigenvectors
359: * (CWorkspace: none)
360: * (RWorkspace: need N)
361: *
362: CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
363: $ IERR )
364: *
365: * Normalize right eigenvectors and make largest component real
366: *
367: DO 40 I = 1, N
368: SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
369: CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
370: DO 30 K = 1, N
371: RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
372: $ DIMAG( VR( K, I ) )**2
373: 30 CONTINUE
374: K = IDAMAX( N, RWORK( IRWORK ), 1 )
375: TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
376: CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
377: VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
378: 40 CONTINUE
379: END IF
380: *
381: * Undo scaling if necessary
382: *
383: 50 CONTINUE
384: IF( SCALEA ) THEN
385: CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
386: $ MAX( N-INFO, 1 ), IERR )
387: IF( INFO.GT.0 ) THEN
388: CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
389: END IF
390: END IF
391: *
392: WORK( 1 ) = MAXWRK
393: RETURN
394: *
395: * End of ZGEEV
396: *
397: 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