1: *> \brief <b> DGEESX 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 DGEESX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeesx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeesx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeesx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
22: * WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
23: * IWORK, LIWORK, BWORK, INFO )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER JOBVS, SENSE, SORT
27: * INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
28: * DOUBLE PRECISION RCONDE, RCONDV
29: * ..
30: * .. Array Arguments ..
31: * LOGICAL BWORK( * )
32: * INTEGER IWORK( * )
33: * DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
34: * $ WR( * )
35: * ..
36: * .. Function Arguments ..
37: * LOGICAL SELECT
38: * EXTERNAL SELECT
39: * ..
40: *
41: *
42: *> \par Purpose:
43: * =============
44: *>
45: *> \verbatim
46: *>
47: *> DGEESX computes for an N-by-N real nonsymmetric matrix A, the
48: *> eigenvalues, the real Schur form T, and, optionally, the matrix of
49: *> Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
50: *>
51: *> Optionally, it also orders the eigenvalues on the diagonal of the
52: *> real Schur form so that selected eigenvalues are at the top left;
53: *> computes a reciprocal condition number for the average of the
54: *> selected eigenvalues (RCONDE); and computes a reciprocal condition
55: *> number for the right invariant subspace corresponding to the
56: *> selected eigenvalues (RCONDV). The leading columns of Z form an
57: *> orthonormal basis for this invariant subspace.
58: *>
59: *> For further explanation of the reciprocal condition numbers RCONDE
60: *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
61: *> these quantities are called s and sep respectively).
62: *>
63: *> A real matrix is in real Schur form if it is upper quasi-triangular
64: *> with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
65: *> the form
66: *> [ a b ]
67: *> [ c a ]
68: *>
69: *> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
70: *> \endverbatim
71: *
72: * Arguments:
73: * ==========
74: *
75: *> \param[in] JOBVS
76: *> \verbatim
77: *> JOBVS is CHARACTER*1
78: *> = 'N': Schur vectors are not computed;
79: *> = 'V': Schur vectors are computed.
80: *> \endverbatim
81: *>
82: *> \param[in] SORT
83: *> \verbatim
84: *> SORT is CHARACTER*1
85: *> Specifies whether or not to order the eigenvalues on the
86: *> diagonal of the Schur form.
87: *> = 'N': Eigenvalues are not ordered;
88: *> = 'S': Eigenvalues are ordered (see SELECT).
89: *> \endverbatim
90: *>
91: *> \param[in] SELECT
92: *> \verbatim
93: *> SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments
94: *> SELECT must be declared EXTERNAL in the calling subroutine.
95: *> If SORT = 'S', SELECT is used to select eigenvalues to sort
96: *> to the top left of the Schur form.
97: *> If SORT = 'N', SELECT is not referenced.
98: *> An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
99: *> SELECT(WR(j),WI(j)) is true; i.e., if either one of a
100: *> complex conjugate pair of eigenvalues is selected, then both
101: *> are. Note that a selected complex eigenvalue may no longer
102: *> satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
103: *> ordering may change the value of complex eigenvalues
104: *> (especially if the eigenvalue is ill-conditioned); in this
105: *> case INFO may be set to N+3 (see INFO below).
106: *> \endverbatim
107: *>
108: *> \param[in] SENSE
109: *> \verbatim
110: *> SENSE is CHARACTER*1
111: *> Determines which reciprocal condition numbers are computed.
112: *> = 'N': None are computed;
113: *> = 'E': Computed for average of selected eigenvalues only;
114: *> = 'V': Computed for selected right invariant subspace only;
115: *> = 'B': Computed for both.
116: *> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
117: *> \endverbatim
118: *>
119: *> \param[in] N
120: *> \verbatim
121: *> N is INTEGER
122: *> The order of the matrix A. N >= 0.
123: *> \endverbatim
124: *>
125: *> \param[in,out] A
126: *> \verbatim
127: *> A is DOUBLE PRECISION array, dimension (LDA, N)
128: *> On entry, the N-by-N matrix A.
129: *> On exit, A is overwritten by its real Schur form T.
130: *> \endverbatim
131: *>
132: *> \param[in] LDA
133: *> \verbatim
134: *> LDA is INTEGER
135: *> The leading dimension of the array A. LDA >= max(1,N).
136: *> \endverbatim
137: *>
138: *> \param[out] SDIM
139: *> \verbatim
140: *> SDIM is INTEGER
141: *> If SORT = 'N', SDIM = 0.
142: *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
143: *> for which SELECT is true. (Complex conjugate
144: *> pairs for which SELECT is true for either
145: *> eigenvalue count as 2.)
146: *> \endverbatim
147: *>
148: *> \param[out] WR
149: *> \verbatim
150: *> WR is DOUBLE PRECISION array, dimension (N)
151: *> \endverbatim
152: *>
153: *> \param[out] WI
154: *> \verbatim
155: *> WI is DOUBLE PRECISION array, dimension (N)
156: *> WR and WI contain the real and imaginary parts, respectively,
157: *> of the computed eigenvalues, in the same order that they
158: *> appear on the diagonal of the output Schur form T. Complex
159: *> conjugate pairs of eigenvalues appear consecutively with the
160: *> eigenvalue having the positive imaginary part first.
161: *> \endverbatim
162: *>
163: *> \param[out] VS
164: *> \verbatim
165: *> VS is DOUBLE PRECISION array, dimension (LDVS,N)
166: *> If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
167: *> vectors.
168: *> If JOBVS = 'N', VS is not referenced.
169: *> \endverbatim
170: *>
171: *> \param[in] LDVS
172: *> \verbatim
173: *> LDVS is INTEGER
174: *> The leading dimension of the array VS. LDVS >= 1, and if
175: *> JOBVS = 'V', LDVS >= N.
176: *> \endverbatim
177: *>
178: *> \param[out] RCONDE
179: *> \verbatim
180: *> RCONDE is DOUBLE PRECISION
181: *> If SENSE = 'E' or 'B', RCONDE contains the reciprocal
182: *> condition number for the average of the selected eigenvalues.
183: *> Not referenced if SENSE = 'N' or 'V'.
184: *> \endverbatim
185: *>
186: *> \param[out] RCONDV
187: *> \verbatim
188: *> RCONDV is DOUBLE PRECISION
189: *> If SENSE = 'V' or 'B', RCONDV contains the reciprocal
190: *> condition number for the selected right invariant subspace.
191: *> Not referenced if SENSE = 'N' or 'E'.
192: *> \endverbatim
193: *>
194: *> \param[out] WORK
195: *> \verbatim
196: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
197: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
198: *> \endverbatim
199: *>
200: *> \param[in] LWORK
201: *> \verbatim
202: *> LWORK is INTEGER
203: *> The dimension of the array WORK. LWORK >= max(1,3*N).
204: *> Also, if SENSE = 'E' or 'V' or 'B',
205: *> LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
206: *> selected eigenvalues computed by this routine. Note that
207: *> N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
208: *> returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
209: *> 'B' this may not be large enough.
210: *> For good performance, LWORK must generally be larger.
211: *>
212: *> If LWORK = -1, then a workspace query is assumed; the routine
213: *> only calculates upper bounds on the optimal sizes of the
214: *> arrays WORK and IWORK, returns these values as the first
215: *> entries of the WORK and IWORK arrays, and no error messages
216: *> related to LWORK or LIWORK are issued by XERBLA.
217: *> \endverbatim
218: *>
219: *> \param[out] IWORK
220: *> \verbatim
221: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
222: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
223: *> \endverbatim
224: *>
225: *> \param[in] LIWORK
226: *> \verbatim
227: *> LIWORK is INTEGER
228: *> The dimension of the array IWORK.
229: *> LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
230: *> Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
231: *> only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
232: *> may not be large enough.
233: *>
234: *> If LIWORK = -1, then a workspace query is assumed; the
235: *> routine only calculates upper bounds on the optimal sizes of
236: *> the arrays WORK and IWORK, returns these values as the first
237: *> entries of the WORK and IWORK arrays, and no error messages
238: *> related to LWORK or LIWORK are issued by XERBLA.
239: *> \endverbatim
240: *>
241: *> \param[out] BWORK
242: *> \verbatim
243: *> BWORK is LOGICAL array, dimension (N)
244: *> Not referenced if SORT = 'N'.
245: *> \endverbatim
246: *>
247: *> \param[out] INFO
248: *> \verbatim
249: *> INFO is INTEGER
250: *> = 0: successful exit
251: *> < 0: if INFO = -i, the i-th argument had an illegal value.
252: *> > 0: if INFO = i, and i is
253: *> <= N: the QR algorithm failed to compute all the
254: *> eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
255: *> contain those eigenvalues which have converged; if
256: *> JOBVS = 'V', VS contains the transformation which
257: *> reduces A to its partially converged Schur form.
258: *> = N+1: the eigenvalues could not be reordered because some
259: *> eigenvalues were too close to separate (the problem
260: *> is very ill-conditioned);
261: *> = N+2: after reordering, roundoff changed values of some
262: *> complex eigenvalues so that leading eigenvalues in
263: *> the Schur form no longer satisfy SELECT=.TRUE. This
264: *> could also be caused by underflow due to scaling.
265: *> \endverbatim
266: *
267: * Authors:
268: * ========
269: *
270: *> \author Univ. of Tennessee
271: *> \author Univ. of California Berkeley
272: *> \author Univ. of Colorado Denver
273: *> \author NAG Ltd.
274: *
275: *> \ingroup doubleGEeigen
276: *
277: * =====================================================================
278: SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
279: $ WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
280: $ IWORK, LIWORK, BWORK, INFO )
281: *
282: * -- LAPACK driver routine --
283: * -- LAPACK is a software package provided by Univ. of Tennessee, --
284: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
285: *
286: * .. Scalar Arguments ..
287: CHARACTER JOBVS, SENSE, SORT
288: INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
289: DOUBLE PRECISION RCONDE, RCONDV
290: * ..
291: * .. Array Arguments ..
292: LOGICAL BWORK( * )
293: INTEGER IWORK( * )
294: DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
295: $ WR( * )
296: * ..
297: * .. Function Arguments ..
298: LOGICAL SELECT
299: EXTERNAL SELECT
300: * ..
301: *
302: * =====================================================================
303: *
304: * .. Parameters ..
305: DOUBLE PRECISION ZERO, ONE
306: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
307: * ..
308: * .. Local Scalars ..
309: LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
310: $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
311: INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
312: $ IHI, ILO, INXT, IP, ITAU, IWRK, LIWRK, LWRK,
313: $ MAXWRK, MINWRK
314: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315: * ..
316: * .. Local Arrays ..
317: DOUBLE PRECISION DUM( 1 )
318: * ..
319: * .. External Subroutines ..
320: EXTERNAL DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
321: $ DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
322: * ..
323: * .. External Functions ..
324: LOGICAL LSAME
325: INTEGER ILAENV
326: DOUBLE PRECISION DLAMCH, DLANGE
327: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
328: * ..
329: * .. Intrinsic Functions ..
330: INTRINSIC MAX, SQRT
331: * ..
332: * .. Executable Statements ..
333: *
334: * Test the input arguments
335: *
336: INFO = 0
337: WANTVS = LSAME( JOBVS, 'V' )
338: WANTST = LSAME( SORT, 'S' )
339: WANTSN = LSAME( SENSE, 'N' )
340: WANTSE = LSAME( SENSE, 'E' )
341: WANTSV = LSAME( SENSE, 'V' )
342: WANTSB = LSAME( SENSE, 'B' )
343: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
344: *
345: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
346: INFO = -1
347: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
348: INFO = -2
349: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
350: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
351: INFO = -4
352: ELSE IF( N.LT.0 ) THEN
353: INFO = -5
354: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
355: INFO = -7
356: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
357: INFO = -12
358: END IF
359: *
360: * Compute workspace
361: * (Note: Comments in the code beginning "RWorkspace:" describe the
362: * minimal amount of real workspace needed at that point in the
363: * code, as well as the preferred amount for good performance.
364: * IWorkspace refers to integer workspace.
365: * NB refers to the optimal block size for the immediately
366: * following subroutine, as returned by ILAENV.
367: * HSWORK refers to the workspace preferred by DHSEQR, as
368: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
369: * the worst case.
370: * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
371: * depends on SDIM, which is computed by the routine DTRSEN later
372: * in the code.)
373: *
374: IF( INFO.EQ.0 ) THEN
375: LIWRK = 1
376: IF( N.EQ.0 ) THEN
377: MINWRK = 1
378: LWRK = 1
379: ELSE
380: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
381: MINWRK = 3*N
382: *
383: CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
384: $ WORK, -1, IEVAL )
385: HSWORK = INT( WORK( 1 ) )
386: *
387: IF( .NOT.WANTVS ) THEN
388: MAXWRK = MAX( MAXWRK, N + HSWORK )
389: ELSE
390: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
391: $ 'DORGHR', ' ', N, 1, N, -1 ) )
392: MAXWRK = MAX( MAXWRK, N + HSWORK )
393: END IF
394: LWRK = MAXWRK
395: IF( .NOT.WANTSN )
396: $ LWRK = MAX( LWRK, N + ( N*N )/2 )
397: IF( WANTSV .OR. WANTSB )
398: $ LIWRK = ( N*N )/4
399: END IF
400: IWORK( 1 ) = LIWRK
401: WORK( 1 ) = LWRK
402: *
403: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
404: INFO = -16
405: ELSE IF( LIWORK.LT.1 .AND. .NOT.LQUERY ) THEN
406: INFO = -18
407: END IF
408: END IF
409: *
410: IF( INFO.NE.0 ) THEN
411: CALL XERBLA( 'DGEESX', -INFO )
412: RETURN
413: ELSE IF( LQUERY ) THEN
414: RETURN
415: END IF
416: *
417: * Quick return if possible
418: *
419: IF( N.EQ.0 ) THEN
420: SDIM = 0
421: RETURN
422: END IF
423: *
424: * Get machine constants
425: *
426: EPS = DLAMCH( 'P' )
427: SMLNUM = DLAMCH( 'S' )
428: BIGNUM = ONE / SMLNUM
429: CALL DLABAD( SMLNUM, BIGNUM )
430: SMLNUM = SQRT( SMLNUM ) / EPS
431: BIGNUM = ONE / SMLNUM
432: *
433: * Scale A if max element outside range [SMLNUM,BIGNUM]
434: *
435: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
436: SCALEA = .FALSE.
437: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
438: SCALEA = .TRUE.
439: CSCALE = SMLNUM
440: ELSE IF( ANRM.GT.BIGNUM ) THEN
441: SCALEA = .TRUE.
442: CSCALE = BIGNUM
443: END IF
444: IF( SCALEA )
445: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
446: *
447: * Permute the matrix to make it more nearly triangular
448: * (RWorkspace: need N)
449: *
450: IBAL = 1
451: CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
452: *
453: * Reduce to upper Hessenberg form
454: * (RWorkspace: need 3*N, prefer 2*N+N*NB)
455: *
456: ITAU = N + IBAL
457: IWRK = N + ITAU
458: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
459: $ LWORK-IWRK+1, IERR )
460: *
461: IF( WANTVS ) THEN
462: *
463: * Copy Householder vectors to VS
464: *
465: CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
466: *
467: * Generate orthogonal matrix in VS
468: * (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
469: *
470: CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
471: $ LWORK-IWRK+1, IERR )
472: END IF
473: *
474: SDIM = 0
475: *
476: * Perform QR iteration, accumulating Schur vectors in VS if desired
477: * (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
478: *
479: IWRK = ITAU
480: CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
481: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
482: IF( IEVAL.GT.0 )
483: $ INFO = IEVAL
484: *
485: * Sort eigenvalues if desired
486: *
487: IF( WANTST .AND. INFO.EQ.0 ) THEN
488: IF( SCALEA ) THEN
489: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
490: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
491: END IF
492: DO 10 I = 1, N
493: BWORK( I ) = SELECT( WR( I ), WI( I ) )
494: 10 CONTINUE
495: *
496: * Reorder eigenvalues, transform Schur vectors, and compute
497: * reciprocal condition numbers
498: * (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
499: * otherwise, need N )
500: * (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
501: * otherwise, need 0 )
502: *
503: CALL DTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
504: $ SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
505: $ IWORK, LIWORK, ICOND )
506: IF( .NOT.WANTSN )
507: $ MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
508: IF( ICOND.EQ.-15 ) THEN
509: *
510: * Not enough real workspace
511: *
512: INFO = -16
513: ELSE IF( ICOND.EQ.-17 ) THEN
514: *
515: * Not enough integer workspace
516: *
517: INFO = -18
518: ELSE IF( ICOND.GT.0 ) THEN
519: *
520: * DTRSEN failed to reorder or to restore standard Schur form
521: *
522: INFO = ICOND + N
523: END IF
524: END IF
525: *
526: IF( WANTVS ) THEN
527: *
528: * Undo balancing
529: * (RWorkspace: need N)
530: *
531: CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
532: $ IERR )
533: END IF
534: *
535: IF( SCALEA ) THEN
536: *
537: * Undo scaling for the Schur form of A
538: *
539: CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
540: CALL DCOPY( N, A, LDA+1, WR, 1 )
541: IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
542: DUM( 1 ) = RCONDV
543: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
544: RCONDV = DUM( 1 )
545: END IF
546: IF( CSCALE.EQ.SMLNUM ) THEN
547: *
548: * If scaling back towards underflow, adjust WI if an
549: * offdiagonal element of a 2-by-2 block in the Schur form
550: * underflows.
551: *
552: IF( IEVAL.GT.0 ) THEN
553: I1 = IEVAL + 1
554: I2 = IHI - 1
555: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
556: $ IERR )
557: ELSE IF( WANTST ) THEN
558: I1 = 1
559: I2 = N - 1
560: ELSE
561: I1 = ILO
562: I2 = IHI - 1
563: END IF
564: INXT = I1 - 1
565: DO 20 I = I1, I2
566: IF( I.LT.INXT )
567: $ GO TO 20
568: IF( WI( I ).EQ.ZERO ) THEN
569: INXT = I + 1
570: ELSE
571: IF( A( I+1, I ).EQ.ZERO ) THEN
572: WI( I ) = ZERO
573: WI( I+1 ) = ZERO
574: ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
575: $ ZERO ) THEN
576: WI( I ) = ZERO
577: WI( I+1 ) = ZERO
578: IF( I.GT.1 )
579: $ CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
580: IF( N.GT.I+1 )
581: $ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
582: $ A( I+1, I+2 ), LDA )
583: IF( WANTVS ) THEN
584: CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
585: END IF
586: A( I, I+1 ) = A( I+1, I )
587: A( I+1, I ) = ZERO
588: END IF
589: INXT = I + 2
590: END IF
591: 20 CONTINUE
592: END IF
593: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
594: $ WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
595: END IF
596: *
597: IF( WANTST .AND. INFO.EQ.0 ) THEN
598: *
599: * Check if reordering successful
600: *
601: LASTSL = .TRUE.
602: LST2SL = .TRUE.
603: SDIM = 0
604: IP = 0
605: DO 30 I = 1, N
606: CURSL = SELECT( WR( I ), WI( I ) )
607: IF( WI( I ).EQ.ZERO ) THEN
608: IF( CURSL )
609: $ SDIM = SDIM + 1
610: IP = 0
611: IF( CURSL .AND. .NOT.LASTSL )
612: $ INFO = N + 2
613: ELSE
614: IF( IP.EQ.1 ) THEN
615: *
616: * Last eigenvalue of conjugate pair
617: *
618: CURSL = CURSL .OR. LASTSL
619: LASTSL = CURSL
620: IF( CURSL )
621: $ SDIM = SDIM + 2
622: IP = -1
623: IF( CURSL .AND. .NOT.LST2SL )
624: $ INFO = N + 2
625: ELSE
626: *
627: * First eigenvalue of conjugate pair
628: *
629: IP = 1
630: END IF
631: END IF
632: LST2SL = LASTSL
633: LASTSL = CURSL
634: 30 CONTINUE
635: END IF
636: *
637: WORK( 1 ) = MAXWRK
638: IF( WANTSV .OR. WANTSB ) THEN
639: IWORK( 1 ) = MAX( 1, SDIM*( N-SDIM ) )
640: ELSE
641: IWORK( 1 ) = 1
642: END IF
643: *
644: RETURN
645: *
646: * End of DGEESX
647: *
648: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgeesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack