Annotation of rpl/lapack/lapack/zgees.f, revision 1.7
1.1 bertrand 1: SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
2: $ LDVS, WORK, LWORK, RWORK, BWORK, 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 JOBVS, SORT
11: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
12: * ..
13: * .. Array Arguments ..
14: LOGICAL BWORK( * )
15: DOUBLE PRECISION RWORK( * )
16: COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
17: * ..
18: * .. Function Arguments ..
19: LOGICAL SELECT
20: EXTERNAL SELECT
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
27: * eigenvalues, the Schur form T, and, optionally, the matrix of Schur
28: * vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
29: *
30: * Optionally, it also orders the eigenvalues on the diagonal of the
31: * Schur form so that selected eigenvalues are at the top left.
32: * The leading columns of Z then form an orthonormal basis for the
33: * invariant subspace corresponding to the selected eigenvalues.
34: *
35: * A complex matrix is in Schur form if it is upper triangular.
36: *
37: * Arguments
38: * =========
39: *
40: * JOBVS (input) CHARACTER*1
41: * = 'N': Schur vectors are not computed;
42: * = 'V': Schur vectors are computed.
43: *
44: * SORT (input) CHARACTER*1
45: * Specifies whether or not to order the eigenvalues on the
46: * diagonal of the Schur form.
47: * = 'N': Eigenvalues are not ordered:
48: * = 'S': Eigenvalues are ordered (see SELECT).
49: *
50: * SELECT (external procedure) LOGICAL FUNCTION of one COMPLEX*16 argument
51: * SELECT must be declared EXTERNAL in the calling subroutine.
52: * If SORT = 'S', SELECT is used to select eigenvalues to order
53: * to the top left of the Schur form.
54: * IF SORT = 'N', SELECT is not referenced.
55: * The eigenvalue W(j) is selected if SELECT(W(j)) is true.
56: *
57: * N (input) INTEGER
58: * The order of the matrix A. N >= 0.
59: *
60: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
61: * On entry, the N-by-N matrix A.
62: * On exit, A has been overwritten by its Schur form T.
63: *
64: * LDA (input) INTEGER
65: * The leading dimension of the array A. LDA >= max(1,N).
66: *
67: * SDIM (output) INTEGER
68: * If SORT = 'N', SDIM = 0.
69: * If SORT = 'S', SDIM = number of eigenvalues for which
70: * SELECT is true.
71: *
72: * W (output) COMPLEX*16 array, dimension (N)
73: * W contains the computed eigenvalues, in the same order that
74: * they appear on the diagonal of the output Schur form T.
75: *
76: * VS (output) COMPLEX*16 array, dimension (LDVS,N)
77: * If JOBVS = 'V', VS contains the unitary matrix Z of Schur
78: * vectors.
79: * If JOBVS = 'N', VS is not referenced.
80: *
81: * LDVS (input) INTEGER
82: * The leading dimension of the array VS. LDVS >= 1; if
83: * JOBVS = 'V', LDVS >= N.
84: *
85: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
86: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
87: *
88: * LWORK (input) INTEGER
89: * The dimension of the array WORK. LWORK >= max(1,2*N).
90: * For good performance, LWORK must generally be larger.
91: *
92: * If LWORK = -1, then a workspace query is assumed; the routine
93: * only calculates the optimal size of the WORK array, returns
94: * this value as the first entry of the WORK array, and no error
95: * message related to LWORK is issued by XERBLA.
96: *
97: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
98: *
99: * BWORK (workspace) LOGICAL array, dimension (N)
100: * Not referenced if SORT = 'N'.
101: *
102: * INFO (output) INTEGER
103: * = 0: successful exit
104: * < 0: if INFO = -i, the i-th argument had an illegal value.
105: * > 0: if INFO = i, and i is
106: * <= N: the QR algorithm failed to compute all the
107: * eigenvalues; elements 1:ILO-1 and i+1:N of W
108: * contain those eigenvalues which have converged;
109: * if JOBVS = 'V', VS contains the matrix which
110: * reduces A to its partially converged Schur form.
111: * = N+1: the eigenvalues could not be reordered because
112: * some eigenvalues were too close to separate (the
113: * problem is very ill-conditioned);
114: * = N+2: after reordering, roundoff changed values of
115: * some complex eigenvalues so that leading
116: * eigenvalues in the Schur form no longer satisfy
117: * SELECT = .TRUE.. This could also be caused by
118: * underflow due to scaling.
119: *
120: * =====================================================================
121: *
122: * .. Parameters ..
123: DOUBLE PRECISION ZERO, ONE
124: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
125: * ..
126: * .. Local Scalars ..
127: LOGICAL LQUERY, SCALEA, WANTST, WANTVS
128: INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
129: $ ITAU, IWRK, MAXWRK, MINWRK
130: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
131: * ..
132: * .. Local Arrays ..
133: DOUBLE PRECISION DUM( 1 )
134: * ..
135: * .. External Subroutines ..
136: EXTERNAL DLABAD, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, ZGEHRD,
137: $ ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
138: * ..
139: * .. External Functions ..
140: LOGICAL LSAME
141: INTEGER ILAENV
142: DOUBLE PRECISION DLAMCH, ZLANGE
143: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
144: * ..
145: * .. Intrinsic Functions ..
146: INTRINSIC MAX, SQRT
147: * ..
148: * .. Executable Statements ..
149: *
150: * Test the input arguments
151: *
152: INFO = 0
153: LQUERY = ( LWORK.EQ.-1 )
154: WANTVS = LSAME( JOBVS, 'V' )
155: WANTST = LSAME( SORT, 'S' )
156: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
157: INFO = -1
158: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
159: INFO = -2
160: ELSE IF( N.LT.0 ) THEN
161: INFO = -4
162: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
163: INFO = -6
164: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
165: INFO = -10
166: END IF
167: *
168: * Compute workspace
169: * (Note: Comments in the code beginning "Workspace:" describe the
170: * minimal amount of workspace needed at that point in the code,
171: * as well as the preferred amount for good performance.
172: * CWorkspace refers to complex workspace, and RWorkspace to real
173: * workspace. NB refers to the optimal block size for the
174: * immediately following subroutine, as returned by ILAENV.
175: * HSWORK refers to the workspace preferred by ZHSEQR, 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 = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
185: MINWRK = 2*N
186: *
187: CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
188: $ WORK, -1, IEVAL )
189: HSWORK = WORK( 1 )
190: *
191: IF( .NOT.WANTVS ) THEN
192: MAXWRK = MAX( MAXWRK, HSWORK )
193: ELSE
194: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
195: $ ' ', N, 1, N, -1 ) )
196: MAXWRK = MAX( MAXWRK, HSWORK )
197: END IF
198: END IF
199: WORK( 1 ) = MAXWRK
200: *
201: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
202: INFO = -12
203: END IF
204: END IF
205: *
206: IF( INFO.NE.0 ) THEN
207: CALL XERBLA( 'ZGEES ', -INFO )
208: RETURN
209: ELSE IF( LQUERY ) THEN
210: RETURN
211: END IF
212: *
213: * Quick return if possible
214: *
215: IF( N.EQ.0 ) THEN
216: SDIM = 0
217: RETURN
218: END IF
219: *
220: * Get machine constants
221: *
222: EPS = DLAMCH( 'P' )
223: SMLNUM = DLAMCH( 'S' )
224: BIGNUM = ONE / SMLNUM
225: CALL DLABAD( SMLNUM, BIGNUM )
226: SMLNUM = SQRT( SMLNUM ) / EPS
227: BIGNUM = ONE / SMLNUM
228: *
229: * Scale A if max element outside range [SMLNUM,BIGNUM]
230: *
231: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
232: SCALEA = .FALSE.
233: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
234: SCALEA = .TRUE.
235: CSCALE = SMLNUM
236: ELSE IF( ANRM.GT.BIGNUM ) THEN
237: SCALEA = .TRUE.
238: CSCALE = BIGNUM
239: END IF
240: IF( SCALEA )
241: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
242: *
243: * Permute the matrix to make it more nearly triangular
244: * (CWorkspace: none)
245: * (RWorkspace: need N)
246: *
247: IBAL = 1
248: CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
249: *
250: * Reduce to upper Hessenberg form
251: * (CWorkspace: need 2*N, prefer N+N*NB)
252: * (RWorkspace: none)
253: *
254: ITAU = 1
255: IWRK = N + ITAU
256: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
257: $ LWORK-IWRK+1, IERR )
258: *
259: IF( WANTVS ) THEN
260: *
261: * Copy Householder vectors to VS
262: *
263: CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
264: *
265: * Generate unitary matrix in VS
266: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
267: * (RWorkspace: none)
268: *
269: CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
270: $ LWORK-IWRK+1, IERR )
271: END IF
272: *
273: SDIM = 0
274: *
275: * Perform QR iteration, accumulating Schur vectors in VS if desired
276: * (CWorkspace: need 1, prefer HSWORK (see comments) )
277: * (RWorkspace: none)
278: *
279: IWRK = ITAU
280: CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
281: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
282: IF( IEVAL.GT.0 )
283: $ INFO = IEVAL
284: *
285: * Sort eigenvalues if desired
286: *
287: IF( WANTST .AND. INFO.EQ.0 ) THEN
288: IF( SCALEA )
289: $ CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
290: DO 10 I = 1, N
291: BWORK( I ) = SELECT( W( I ) )
292: 10 CONTINUE
293: *
294: * Reorder eigenvalues and transform Schur vectors
295: * (CWorkspace: none)
296: * (RWorkspace: none)
297: *
298: CALL ZTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
299: $ S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
300: END IF
301: *
302: IF( WANTVS ) THEN
303: *
304: * Undo balancing
305: * (CWorkspace: none)
306: * (RWorkspace: need N)
307: *
308: CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
309: $ IERR )
310: END IF
311: *
312: IF( SCALEA ) THEN
313: *
314: * Undo scaling for the Schur form of A
315: *
316: CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
317: CALL ZCOPY( N, A, LDA+1, W, 1 )
318: END IF
319: *
320: WORK( 1 ) = MAXWRK
321: RETURN
322: *
323: * End of ZGEES
324: *
325: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgees.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack