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: *> \date November 2011
212: *
213: *> \ingroup doubleGEeigen
214: *
215: * =====================================================================
216: SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
217: $ VS, LDVS, WORK, LWORK, BWORK, INFO )
218: *
219: * -- LAPACK driver routine (version 3.4.0) --
220: * -- LAPACK is a software package provided by Univ. of Tennessee, --
221: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222: * November 2011
223: *
224: * .. Scalar Arguments ..
225: CHARACTER JOBVS, SORT
226: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
227: * ..
228: * .. Array Arguments ..
229: LOGICAL BWORK( * )
230: DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
231: $ WR( * )
232: * ..
233: * .. Function Arguments ..
234: LOGICAL SELECT
235: EXTERNAL SELECT
236: * ..
237: *
238: * =====================================================================
239: *
240: * .. Parameters ..
241: DOUBLE PRECISION ZERO, ONE
242: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
243: * ..
244: * .. Local Scalars ..
245: LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
246: $ WANTVS
247: INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
248: $ IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
249: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
250: * ..
251: * .. Local Arrays ..
252: INTEGER IDUM( 1 )
253: DOUBLE PRECISION DUM( 1 )
254: * ..
255: * .. External Subroutines ..
256: EXTERNAL DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
257: $ DLABAD, DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
258: * ..
259: * .. External Functions ..
260: LOGICAL LSAME
261: INTEGER ILAENV
262: DOUBLE PRECISION DLAMCH, DLANGE
263: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
264: * ..
265: * .. Intrinsic Functions ..
266: INTRINSIC MAX, SQRT
267: * ..
268: * .. Executable Statements ..
269: *
270: * Test the input arguments
271: *
272: INFO = 0
273: LQUERY = ( LWORK.EQ.-1 )
274: WANTVS = LSAME( JOBVS, 'V' )
275: WANTST = LSAME( SORT, 'S' )
276: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
277: INFO = -1
278: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
279: INFO = -2
280: ELSE IF( N.LT.0 ) THEN
281: INFO = -4
282: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
283: INFO = -6
284: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
285: INFO = -11
286: END IF
287: *
288: * Compute workspace
289: * (Note: Comments in the code beginning "Workspace:" describe the
290: * minimal amount of workspace needed at that point in the code,
291: * as well as the preferred amount for good performance.
292: * NB refers to the optimal block size for the immediately
293: * following subroutine, as returned by ILAENV.
294: * HSWORK refers to the workspace preferred by DHSEQR, as
295: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
296: * the worst case.)
297: *
298: IF( INFO.EQ.0 ) THEN
299: IF( N.EQ.0 ) THEN
300: MINWRK = 1
301: MAXWRK = 1
302: ELSE
303: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
304: MINWRK = 3*N
305: *
306: CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
307: $ WORK, -1, IEVAL )
308: HSWORK = WORK( 1 )
309: *
310: IF( .NOT.WANTVS ) THEN
311: MAXWRK = MAX( MAXWRK, N + HSWORK )
312: ELSE
313: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
314: $ 'DORGHR', ' ', N, 1, N, -1 ) )
315: MAXWRK = MAX( MAXWRK, N + HSWORK )
316: END IF
317: END IF
318: WORK( 1 ) = MAXWRK
319: *
320: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
321: INFO = -13
322: END IF
323: END IF
324: *
325: IF( INFO.NE.0 ) THEN
326: CALL XERBLA( 'DGEES ', -INFO )
327: RETURN
328: ELSE IF( LQUERY ) THEN
329: RETURN
330: END IF
331: *
332: * Quick return if possible
333: *
334: IF( N.EQ.0 ) THEN
335: SDIM = 0
336: RETURN
337: END IF
338: *
339: * Get machine constants
340: *
341: EPS = DLAMCH( 'P' )
342: SMLNUM = DLAMCH( 'S' )
343: BIGNUM = ONE / SMLNUM
344: CALL DLABAD( SMLNUM, BIGNUM )
345: SMLNUM = SQRT( SMLNUM ) / EPS
346: BIGNUM = ONE / SMLNUM
347: *
348: * Scale A if max element outside range [SMLNUM,BIGNUM]
349: *
350: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
351: SCALEA = .FALSE.
352: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
353: SCALEA = .TRUE.
354: CSCALE = SMLNUM
355: ELSE IF( ANRM.GT.BIGNUM ) THEN
356: SCALEA = .TRUE.
357: CSCALE = BIGNUM
358: END IF
359: IF( SCALEA )
360: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
361: *
362: * Permute the matrix to make it more nearly triangular
363: * (Workspace: need N)
364: *
365: IBAL = 1
366: CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
367: *
368: * Reduce to upper Hessenberg form
369: * (Workspace: need 3*N, prefer 2*N+N*NB)
370: *
371: ITAU = N + IBAL
372: IWRK = N + ITAU
373: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
374: $ LWORK-IWRK+1, IERR )
375: *
376: IF( WANTVS ) THEN
377: *
378: * Copy Householder vectors to VS
379: *
380: CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
381: *
382: * Generate orthogonal matrix in VS
383: * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
384: *
385: CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
386: $ LWORK-IWRK+1, IERR )
387: END IF
388: *
389: SDIM = 0
390: *
391: * Perform QR iteration, accumulating Schur vectors in VS if desired
392: * (Workspace: need N+1, prefer N+HSWORK (see comments) )
393: *
394: IWRK = ITAU
395: CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
396: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
397: IF( IEVAL.GT.0 )
398: $ INFO = IEVAL
399: *
400: * Sort eigenvalues if desired
401: *
402: IF( WANTST .AND. INFO.EQ.0 ) THEN
403: IF( SCALEA ) THEN
404: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
405: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
406: END IF
407: DO 10 I = 1, N
408: BWORK( I ) = SELECT( WR( I ), WI( I ) )
409: 10 CONTINUE
410: *
411: * Reorder eigenvalues and transform Schur vectors
412: * (Workspace: none needed)
413: *
414: CALL DTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
415: $ SDIM, S, SEP, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
416: $ ICOND )
417: IF( ICOND.GT.0 )
418: $ INFO = N + ICOND
419: END IF
420: *
421: IF( WANTVS ) THEN
422: *
423: * Undo balancing
424: * (Workspace: need N)
425: *
426: CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
427: $ IERR )
428: END IF
429: *
430: IF( SCALEA ) THEN
431: *
432: * Undo scaling for the Schur form of A
433: *
434: CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
435: CALL DCOPY( N, A, LDA+1, WR, 1 )
436: IF( CSCALE.EQ.SMLNUM ) THEN
437: *
438: * If scaling back towards underflow, adjust WI if an
439: * offdiagonal element of a 2-by-2 block in the Schur form
440: * underflows.
441: *
442: IF( IEVAL.GT.0 ) THEN
443: I1 = IEVAL + 1
444: I2 = IHI - 1
445: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI,
446: $ MAX( ILO-1, 1 ), IERR )
447: ELSE IF( WANTST ) THEN
448: I1 = 1
449: I2 = N - 1
450: ELSE
451: I1 = ILO
452: I2 = IHI - 1
453: END IF
454: INXT = I1 - 1
455: DO 20 I = I1, I2
456: IF( I.LT.INXT )
457: $ GO TO 20
458: IF( WI( I ).EQ.ZERO ) THEN
459: INXT = I + 1
460: ELSE
461: IF( A( I+1, I ).EQ.ZERO ) THEN
462: WI( I ) = ZERO
463: WI( I+1 ) = ZERO
464: ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
465: $ ZERO ) THEN
466: WI( I ) = ZERO
467: WI( I+1 ) = ZERO
468: IF( I.GT.1 )
469: $ CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
470: IF( N.GT.I+1 )
471: $ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
472: $ A( I+1, I+2 ), LDA )
473: IF( WANTVS ) THEN
474: CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
475: END IF
476: A( I, I+1 ) = A( I+1, I )
477: A( I+1, I ) = ZERO
478: END IF
479: INXT = I + 2
480: END IF
481: 20 CONTINUE
482: END IF
483: *
484: * Undo scaling for the imaginary part of the eigenvalues
485: *
486: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
487: $ WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
488: END IF
489: *
490: IF( WANTST .AND. INFO.EQ.0 ) THEN
491: *
492: * Check if reordering successful
493: *
494: LASTSL = .TRUE.
495: LST2SL = .TRUE.
496: SDIM = 0
497: IP = 0
498: DO 30 I = 1, N
499: CURSL = SELECT( WR( I ), WI( I ) )
500: IF( WI( I ).EQ.ZERO ) THEN
501: IF( CURSL )
502: $ SDIM = SDIM + 1
503: IP = 0
504: IF( CURSL .AND. .NOT.LASTSL )
505: $ INFO = N + 2
506: ELSE
507: IF( IP.EQ.1 ) THEN
508: *
509: * Last eigenvalue of conjugate pair
510: *
511: CURSL = CURSL .OR. LASTSL
512: LASTSL = CURSL
513: IF( CURSL )
514: $ SDIM = SDIM + 2
515: IP = -1
516: IF( CURSL .AND. .NOT.LST2SL )
517: $ INFO = N + 2
518: ELSE
519: *
520: * First eigenvalue of conjugate pair
521: *
522: IP = 1
523: END IF
524: END IF
525: LST2SL = LASTSL
526: LASTSL = CURSL
527: 30 CONTINUE
528: END IF
529: *
530: WORK( 1 ) = MAXWRK
531: RETURN
532: *
533: * End of DGEES
534: *
535: 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