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