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