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