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