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