1: *> \brief <b> ZGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors 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 ZGEES + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgees.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgees.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgees.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
22: * LDVS, WORK, LWORK, RWORK, BWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBVS, SORT
26: * INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
27: * ..
28: * .. Array Arguments ..
29: * LOGICAL BWORK( * )
30: * DOUBLE PRECISION RWORK( * )
31: * COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
32: * ..
33: * .. Function Arguments ..
34: * LOGICAL SELECT
35: * EXTERNAL SELECT
36: * ..
37: *
38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> ZGEES computes for an N-by-N complex nonsymmetric matrix A, the
45: *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
46: *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
47: *>
48: *> Optionally, it also orders the eigenvalues on the diagonal of the
49: *> Schur form so that selected eigenvalues are at the top left.
50: *> The leading columns of Z then form an orthonormal basis for the
51: *> invariant subspace corresponding to the selected eigenvalues.
52: *>
53: *> A complex matrix is in Schur form if it is upper triangular.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] JOBVS
60: *> \verbatim
61: *> JOBVS is CHARACTER*1
62: *> = 'N': Schur vectors are not computed;
63: *> = 'V': Schur vectors are computed.
64: *> \endverbatim
65: *>
66: *> \param[in] SORT
67: *> \verbatim
68: *> SORT is CHARACTER*1
69: *> Specifies whether or not to order the eigenvalues on the
70: *> diagonal of the Schur form.
71: *> = 'N': Eigenvalues are not ordered:
72: *> = 'S': Eigenvalues are ordered (see SELECT).
73: *> \endverbatim
74: *>
75: *> \param[in] SELECT
76: *> \verbatim
77: *> SELECT is a LOGICAL FUNCTION of one COMPLEX*16 argument
78: *> SELECT must be declared EXTERNAL in the calling subroutine.
79: *> If SORT = 'S', SELECT is used to select eigenvalues to order
80: *> to the top left of the Schur form.
81: *> IF SORT = 'N', SELECT is not referenced.
82: *> The eigenvalue W(j) is selected if SELECT(W(j)) is true.
83: *> \endverbatim
84: *>
85: *> \param[in] N
86: *> \verbatim
87: *> N is INTEGER
88: *> The order of the matrix A. N >= 0.
89: *> \endverbatim
90: *>
91: *> \param[in,out] A
92: *> \verbatim
93: *> A is COMPLEX*16 array, dimension (LDA,N)
94: *> On entry, the N-by-N matrix A.
95: *> On exit, A has been overwritten by its Schur form T.
96: *> \endverbatim
97: *>
98: *> \param[in] LDA
99: *> \verbatim
100: *> LDA is INTEGER
101: *> The leading dimension of the array A. LDA >= max(1,N).
102: *> \endverbatim
103: *>
104: *> \param[out] SDIM
105: *> \verbatim
106: *> SDIM is INTEGER
107: *> If SORT = 'N', SDIM = 0.
108: *> If SORT = 'S', SDIM = number of eigenvalues for which
109: *> SELECT is true.
110: *> \endverbatim
111: *>
112: *> \param[out] W
113: *> \verbatim
114: *> W is COMPLEX*16 array, dimension (N)
115: *> W contains the computed eigenvalues, in the same order that
116: *> they appear on the diagonal of the output Schur form T.
117: *> \endverbatim
118: *>
119: *> \param[out] VS
120: *> \verbatim
121: *> VS is COMPLEX*16 array, dimension (LDVS,N)
122: *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
123: *> vectors.
124: *> If JOBVS = 'N', VS is not referenced.
125: *> \endverbatim
126: *>
127: *> \param[in] LDVS
128: *> \verbatim
129: *> LDVS is INTEGER
130: *> The leading dimension of the array VS. LDVS >= 1; if
131: *> JOBVS = 'V', LDVS >= N.
132: *> \endverbatim
133: *>
134: *> \param[out] WORK
135: *> \verbatim
136: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
137: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
138: *> \endverbatim
139: *>
140: *> \param[in] LWORK
141: *> \verbatim
142: *> LWORK is INTEGER
143: *> The dimension of the array WORK. LWORK >= max(1,2*N).
144: *> For good performance, LWORK must generally be larger.
145: *>
146: *> If LWORK = -1, then a workspace query is assumed; the routine
147: *> only calculates the optimal size of the WORK array, returns
148: *> this value as the first entry of the WORK array, and no error
149: *> message related to LWORK is issued by XERBLA.
150: *> \endverbatim
151: *>
152: *> \param[out] RWORK
153: *> \verbatim
154: *> RWORK is DOUBLE PRECISION array, dimension (N)
155: *> \endverbatim
156: *>
157: *> \param[out] BWORK
158: *> \verbatim
159: *> BWORK is LOGICAL array, dimension (N)
160: *> Not referenced if SORT = 'N'.
161: *> \endverbatim
162: *>
163: *> \param[out] INFO
164: *> \verbatim
165: *> INFO is INTEGER
166: *> = 0: successful exit
167: *> < 0: if INFO = -i, the i-th argument had an illegal value.
168: *> > 0: if INFO = i, and i is
169: *> <= N: the QR algorithm failed to compute all the
170: *> eigenvalues; elements 1:ILO-1 and i+1:N of W
171: *> contain those eigenvalues which have converged;
172: *> if JOBVS = 'V', VS contains the matrix which
173: *> reduces A to its partially converged Schur form.
174: *> = N+1: the eigenvalues could not be reordered because
175: *> some eigenvalues were too close to separate (the
176: *> problem is very ill-conditioned);
177: *> = N+2: after reordering, roundoff changed values of
178: *> some complex eigenvalues so that leading
179: *> eigenvalues in the Schur form no longer satisfy
180: *> SELECT = .TRUE.. This could also be caused by
181: *> underflow due to scaling.
182: *> \endverbatim
183: *
184: * Authors:
185: * ========
186: *
187: *> \author Univ. of Tennessee
188: *> \author Univ. of California Berkeley
189: *> \author Univ. of Colorado Denver
190: *> \author NAG Ltd.
191: *
192: *> \ingroup complex16GEeigen
193: *
194: * =====================================================================
195: SUBROUTINE ZGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
196: $ LDVS, WORK, LWORK, RWORK, BWORK, INFO )
197: *
198: * -- LAPACK driver routine --
199: * -- LAPACK is a software package provided by Univ. of Tennessee, --
200: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
201: *
202: * .. Scalar Arguments ..
203: CHARACTER JOBVS, SORT
204: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
205: * ..
206: * .. Array Arguments ..
207: LOGICAL BWORK( * )
208: DOUBLE PRECISION RWORK( * )
209: COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
210: * ..
211: * .. Function Arguments ..
212: LOGICAL SELECT
213: EXTERNAL SELECT
214: * ..
215: *
216: * =====================================================================
217: *
218: * .. Parameters ..
219: DOUBLE PRECISION ZERO, ONE
220: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
221: * ..
222: * .. Local Scalars ..
223: LOGICAL LQUERY, SCALEA, WANTST, WANTVS
224: INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
225: $ ITAU, IWRK, MAXWRK, MINWRK
226: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
227: * ..
228: * .. Local Arrays ..
229: DOUBLE PRECISION DUM( 1 )
230: * ..
231: * .. External Subroutines ..
232: EXTERNAL DLABAD, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, ZGEHRD,
233: $ ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
234: * ..
235: * .. External Functions ..
236: LOGICAL LSAME
237: INTEGER ILAENV
238: DOUBLE PRECISION DLAMCH, ZLANGE
239: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
240: * ..
241: * .. Intrinsic Functions ..
242: INTRINSIC MAX, SQRT
243: * ..
244: * .. Executable Statements ..
245: *
246: * Test the input arguments
247: *
248: INFO = 0
249: LQUERY = ( LWORK.EQ.-1 )
250: WANTVS = LSAME( JOBVS, 'V' )
251: WANTST = LSAME( SORT, 'S' )
252: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
253: INFO = -1
254: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
255: INFO = -2
256: ELSE IF( N.LT.0 ) THEN
257: INFO = -4
258: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
259: INFO = -6
260: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
261: INFO = -10
262: END IF
263: *
264: * Compute workspace
265: * (Note: Comments in the code beginning "Workspace:" describe the
266: * minimal amount of workspace needed at that point in the code,
267: * as well as the preferred amount for good performance.
268: * CWorkspace refers to complex workspace, and RWorkspace to real
269: * workspace. NB refers to the optimal block size for the
270: * immediately following subroutine, as returned by ILAENV.
271: * HSWORK refers to the workspace preferred by ZHSEQR, as
272: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
273: * the worst case.)
274: *
275: IF( INFO.EQ.0 ) THEN
276: IF( N.EQ.0 ) THEN
277: MINWRK = 1
278: MAXWRK = 1
279: ELSE
280: MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
281: MINWRK = 2*N
282: *
283: CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
284: $ WORK, -1, IEVAL )
285: HSWORK = INT( WORK( 1 ) )
286: *
287: IF( .NOT.WANTVS ) THEN
288: MAXWRK = MAX( MAXWRK, HSWORK )
289: ELSE
290: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
291: $ ' ', N, 1, N, -1 ) )
292: MAXWRK = MAX( MAXWRK, HSWORK )
293: END IF
294: END IF
295: WORK( 1 ) = MAXWRK
296: *
297: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
298: INFO = -12
299: END IF
300: END IF
301: *
302: IF( INFO.NE.0 ) THEN
303: CALL XERBLA( 'ZGEES ', -INFO )
304: RETURN
305: ELSE IF( LQUERY ) THEN
306: RETURN
307: END IF
308: *
309: * Quick return if possible
310: *
311: IF( N.EQ.0 ) THEN
312: SDIM = 0
313: RETURN
314: END IF
315: *
316: * Get machine constants
317: *
318: EPS = DLAMCH( 'P' )
319: SMLNUM = DLAMCH( 'S' )
320: BIGNUM = ONE / SMLNUM
321: CALL DLABAD( SMLNUM, BIGNUM )
322: SMLNUM = SQRT( SMLNUM ) / EPS
323: BIGNUM = ONE / SMLNUM
324: *
325: * Scale A if max element outside range [SMLNUM,BIGNUM]
326: *
327: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
328: SCALEA = .FALSE.
329: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
330: SCALEA = .TRUE.
331: CSCALE = SMLNUM
332: ELSE IF( ANRM.GT.BIGNUM ) THEN
333: SCALEA = .TRUE.
334: CSCALE = BIGNUM
335: END IF
336: IF( SCALEA )
337: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
338: *
339: * Permute the matrix to make it more nearly triangular
340: * (CWorkspace: none)
341: * (RWorkspace: need N)
342: *
343: IBAL = 1
344: CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
345: *
346: * Reduce to upper Hessenberg form
347: * (CWorkspace: need 2*N, prefer N+N*NB)
348: * (RWorkspace: none)
349: *
350: ITAU = 1
351: IWRK = N + ITAU
352: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
353: $ LWORK-IWRK+1, IERR )
354: *
355: IF( WANTVS ) THEN
356: *
357: * Copy Householder vectors to VS
358: *
359: CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
360: *
361: * Generate unitary matrix in VS
362: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
363: * (RWorkspace: none)
364: *
365: CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
366: $ LWORK-IWRK+1, IERR )
367: END IF
368: *
369: SDIM = 0
370: *
371: * Perform QR iteration, accumulating Schur vectors in VS if desired
372: * (CWorkspace: need 1, prefer HSWORK (see comments) )
373: * (RWorkspace: none)
374: *
375: IWRK = ITAU
376: CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
377: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
378: IF( IEVAL.GT.0 )
379: $ INFO = IEVAL
380: *
381: * Sort eigenvalues if desired
382: *
383: IF( WANTST .AND. INFO.EQ.0 ) THEN
384: IF( SCALEA )
385: $ CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
386: DO 10 I = 1, N
387: BWORK( I ) = SELECT( W( I ) )
388: 10 CONTINUE
389: *
390: * Reorder eigenvalues and transform Schur vectors
391: * (CWorkspace: none)
392: * (RWorkspace: none)
393: *
394: CALL ZTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
395: $ S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
396: END IF
397: *
398: IF( WANTVS ) THEN
399: *
400: * Undo balancing
401: * (CWorkspace: none)
402: * (RWorkspace: need N)
403: *
404: CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
405: $ IERR )
406: END IF
407: *
408: IF( SCALEA ) THEN
409: *
410: * Undo scaling for the Schur form of A
411: *
412: CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
413: CALL ZCOPY( N, A, LDA+1, W, 1 )
414: END IF
415: *
416: WORK( 1 ) = MAXWRK
417: RETURN
418: *
419: * End of ZGEES
420: *
421: 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