Annotation of rpl/lapack/lapack/dgeesx.f, revision 1.11
1.9 bertrand 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 procedure) 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: *> \date November 2011
276: *
277: *> \ingroup doubleGEeigen
278: *
279: * =====================================================================
1.1 bertrand 280: SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
281: $ WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
282: $ IWORK, LIWORK, BWORK, INFO )
283: *
1.9 bertrand 284: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 285: * -- LAPACK is a software package provided by Univ. of Tennessee, --
286: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 287: * November 2011
1.1 bertrand 288: *
289: * .. Scalar Arguments ..
290: CHARACTER JOBVS, SENSE, SORT
291: INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
292: DOUBLE PRECISION RCONDE, RCONDV
293: * ..
294: * .. Array Arguments ..
295: LOGICAL BWORK( * )
296: INTEGER IWORK( * )
297: DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
298: $ WR( * )
299: * ..
300: * .. Function Arguments ..
301: LOGICAL SELECT
302: EXTERNAL SELECT
303: * ..
304: *
305: * =====================================================================
306: *
307: * .. Parameters ..
308: DOUBLE PRECISION ZERO, ONE
309: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
310: * ..
311: * .. Local Scalars ..
312: LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
313: $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
314: INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
315: $ IHI, ILO, INXT, IP, ITAU, IWRK, LIWRK, LWRK,
316: $ MAXWRK, MINWRK
317: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
318: * ..
319: * .. Local Arrays ..
320: DOUBLE PRECISION DUM( 1 )
321: * ..
322: * .. External Subroutines ..
323: EXTERNAL DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
324: $ DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
325: * ..
326: * .. External Functions ..
327: LOGICAL LSAME
328: INTEGER ILAENV
329: DOUBLE PRECISION DLAMCH, DLANGE
330: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
331: * ..
332: * .. Intrinsic Functions ..
333: INTRINSIC MAX, SQRT
334: * ..
335: * .. Executable Statements ..
336: *
337: * Test the input arguments
338: *
339: INFO = 0
340: WANTVS = LSAME( JOBVS, 'V' )
341: WANTST = LSAME( SORT, 'S' )
342: WANTSN = LSAME( SENSE, 'N' )
343: WANTSE = LSAME( SENSE, 'E' )
344: WANTSV = LSAME( SENSE, 'V' )
345: WANTSB = LSAME( SENSE, 'B' )
346: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
1.5 bertrand 347: *
1.1 bertrand 348: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
349: INFO = -1
350: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
351: INFO = -2
352: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
353: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
354: INFO = -4
355: ELSE IF( N.LT.0 ) THEN
356: INFO = -5
357: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
358: INFO = -7
359: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
360: INFO = -12
361: END IF
362: *
363: * Compute workspace
364: * (Note: Comments in the code beginning "RWorkspace:" describe the
365: * minimal amount of real workspace needed at that point in the
366: * code, as well as the preferred amount for good performance.
367: * IWorkspace refers to integer workspace.
368: * NB refers to the optimal block size for the immediately
369: * following subroutine, as returned by ILAENV.
370: * HSWORK refers to the workspace preferred by DHSEQR, as
371: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
372: * the worst case.
373: * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
374: * depends on SDIM, which is computed by the routine DTRSEN later
375: * in the code.)
376: *
377: IF( INFO.EQ.0 ) THEN
378: LIWRK = 1
379: IF( N.EQ.0 ) THEN
380: MINWRK = 1
381: LWRK = 1
382: ELSE
383: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
384: MINWRK = 3*N
385: *
386: CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
387: $ WORK, -1, IEVAL )
388: HSWORK = WORK( 1 )
389: *
390: IF( .NOT.WANTVS ) THEN
391: MAXWRK = MAX( MAXWRK, N + HSWORK )
392: ELSE
393: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
394: $ 'DORGHR', ' ', N, 1, N, -1 ) )
395: MAXWRK = MAX( MAXWRK, N + HSWORK )
396: END IF
397: LWRK = MAXWRK
398: IF( .NOT.WANTSN )
399: $ LWRK = MAX( LWRK, N + ( N*N )/2 )
400: IF( WANTSV .OR. WANTSB )
401: $ LIWRK = ( N*N )/4
402: END IF
403: IWORK( 1 ) = LIWRK
404: WORK( 1 ) = LWRK
405: *
406: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
407: INFO = -16
408: ELSE IF( LIWORK.LT.1 .AND. .NOT.LQUERY ) THEN
409: INFO = -18
410: END IF
411: END IF
412: *
413: IF( INFO.NE.0 ) THEN
414: CALL XERBLA( 'DGEESX', -INFO )
415: RETURN
1.5 bertrand 416: ELSE IF( LQUERY ) THEN
417: RETURN
1.1 bertrand 418: END IF
419: *
420: * Quick return if possible
421: *
422: IF( N.EQ.0 ) THEN
423: SDIM = 0
424: RETURN
425: END IF
426: *
427: * Get machine constants
428: *
429: EPS = DLAMCH( 'P' )
430: SMLNUM = DLAMCH( 'S' )
431: BIGNUM = ONE / SMLNUM
432: CALL DLABAD( SMLNUM, BIGNUM )
433: SMLNUM = SQRT( SMLNUM ) / EPS
434: BIGNUM = ONE / SMLNUM
435: *
436: * Scale A if max element outside range [SMLNUM,BIGNUM]
437: *
438: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
439: SCALEA = .FALSE.
440: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
441: SCALEA = .TRUE.
442: CSCALE = SMLNUM
443: ELSE IF( ANRM.GT.BIGNUM ) THEN
444: SCALEA = .TRUE.
445: CSCALE = BIGNUM
446: END IF
447: IF( SCALEA )
448: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
449: *
450: * Permute the matrix to make it more nearly triangular
451: * (RWorkspace: need N)
452: *
453: IBAL = 1
454: CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
455: *
456: * Reduce to upper Hessenberg form
457: * (RWorkspace: need 3*N, prefer 2*N+N*NB)
458: *
459: ITAU = N + IBAL
460: IWRK = N + ITAU
461: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
462: $ LWORK-IWRK+1, IERR )
463: *
464: IF( WANTVS ) THEN
465: *
466: * Copy Householder vectors to VS
467: *
468: CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
469: *
470: * Generate orthogonal matrix in VS
471: * (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
472: *
473: CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
474: $ LWORK-IWRK+1, IERR )
475: END IF
476: *
477: SDIM = 0
478: *
479: * Perform QR iteration, accumulating Schur vectors in VS if desired
480: * (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
481: *
482: IWRK = ITAU
483: CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
484: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
485: IF( IEVAL.GT.0 )
486: $ INFO = IEVAL
487: *
488: * Sort eigenvalues if desired
489: *
490: IF( WANTST .AND. INFO.EQ.0 ) THEN
491: IF( SCALEA ) THEN
492: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
493: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
494: END IF
495: DO 10 I = 1, N
496: BWORK( I ) = SELECT( WR( I ), WI( I ) )
497: 10 CONTINUE
498: *
499: * Reorder eigenvalues, transform Schur vectors, and compute
500: * reciprocal condition numbers
501: * (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
502: * otherwise, need N )
503: * (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
504: * otherwise, need 0 )
505: *
506: CALL DTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
507: $ SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
508: $ IWORK, LIWORK, ICOND )
509: IF( .NOT.WANTSN )
510: $ MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
511: IF( ICOND.EQ.-15 ) THEN
512: *
513: * Not enough real workspace
514: *
515: INFO = -16
516: ELSE IF( ICOND.EQ.-17 ) THEN
517: *
518: * Not enough integer workspace
519: *
520: INFO = -18
521: ELSE IF( ICOND.GT.0 ) THEN
522: *
523: * DTRSEN failed to reorder or to restore standard Schur form
524: *
525: INFO = ICOND + N
526: END IF
527: END IF
528: *
529: IF( WANTVS ) THEN
530: *
531: * Undo balancing
532: * (RWorkspace: need N)
533: *
534: CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
535: $ IERR )
536: END IF
537: *
538: IF( SCALEA ) THEN
539: *
540: * Undo scaling for the Schur form of A
541: *
542: CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
543: CALL DCOPY( N, A, LDA+1, WR, 1 )
544: IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
545: DUM( 1 ) = RCONDV
546: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
547: RCONDV = DUM( 1 )
548: END IF
549: IF( CSCALE.EQ.SMLNUM ) THEN
550: *
551: * If scaling back towards underflow, adjust WI if an
552: * offdiagonal element of a 2-by-2 block in the Schur form
553: * underflows.
554: *
555: IF( IEVAL.GT.0 ) THEN
556: I1 = IEVAL + 1
557: I2 = IHI - 1
558: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
559: $ IERR )
560: ELSE IF( WANTST ) THEN
561: I1 = 1
562: I2 = N - 1
563: ELSE
564: I1 = ILO
565: I2 = IHI - 1
566: END IF
567: INXT = I1 - 1
568: DO 20 I = I1, I2
569: IF( I.LT.INXT )
570: $ GO TO 20
571: IF( WI( I ).EQ.ZERO ) THEN
572: INXT = I + 1
573: ELSE
574: IF( A( I+1, I ).EQ.ZERO ) THEN
575: WI( I ) = ZERO
576: WI( I+1 ) = ZERO
577: ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
578: $ ZERO ) THEN
579: WI( I ) = ZERO
580: WI( I+1 ) = ZERO
581: IF( I.GT.1 )
582: $ CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
583: IF( N.GT.I+1 )
584: $ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
585: $ A( I+1, I+2 ), LDA )
586: CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
587: A( I, I+1 ) = A( I+1, I )
588: A( I+1, I ) = ZERO
589: END IF
590: INXT = I + 2
591: END IF
592: 20 CONTINUE
593: END IF
594: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
595: $ WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
596: END IF
597: *
598: IF( WANTST .AND. INFO.EQ.0 ) THEN
599: *
600: * Check if reordering successful
601: *
602: LASTSL = .TRUE.
603: LST2SL = .TRUE.
604: SDIM = 0
605: IP = 0
606: DO 30 I = 1, N
607: CURSL = SELECT( WR( I ), WI( I ) )
608: IF( WI( I ).EQ.ZERO ) THEN
609: IF( CURSL )
610: $ SDIM = SDIM + 1
611: IP = 0
612: IF( CURSL .AND. .NOT.LASTSL )
613: $ INFO = N + 2
614: ELSE
615: IF( IP.EQ.1 ) THEN
616: *
617: * Last eigenvalue of conjugate pair
618: *
619: CURSL = CURSL .OR. LASTSL
620: LASTSL = CURSL
621: IF( CURSL )
622: $ SDIM = SDIM + 2
623: IP = -1
624: IF( CURSL .AND. .NOT.LST2SL )
625: $ INFO = N + 2
626: ELSE
627: *
628: * First eigenvalue of conjugate pair
629: *
630: IP = 1
631: END IF
632: END IF
633: LST2SL = LASTSL
634: LASTSL = CURSL
635: 30 CONTINUE
636: END IF
637: *
638: WORK( 1 ) = MAXWRK
639: IF( WANTSV .OR. WANTSB ) THEN
640: IWORK( 1 ) = MAX( 1, SDIM*( N-SDIM ) )
641: ELSE
642: IWORK( 1 ) = 1
643: END IF
644: *
645: RETURN
646: *
647: * End of DGEESX
648: *
649: 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