1: *> \brief <b> DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors 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 DGEGS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegs.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegs.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegs.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
22: * ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
23: * LWORK, INFO )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER JOBVSL, JOBVSR
27: * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
28: * ..
29: * .. Array Arguments ..
30: * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
31: * $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
32: * $ VSR( LDVSR, * ), WORK( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> This routine is deprecated and has been replaced by routine DGGES.
42: *>
43: *> DGEGS computes the eigenvalues, real Schur form, and, optionally,
44: *> left and or/right Schur vectors of a real matrix pair (A,B).
45: *> Given two square matrices A and B, the generalized real Schur
46: *> factorization has the form
47: *>
48: *> A = Q*S*Z**T, B = Q*T*Z**T
49: *>
50: *> where Q and Z are orthogonal matrices, T is upper triangular, and S
51: *> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
52: *> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
53: *> of eigenvalues of (A,B). The columns of Q are the left Schur vectors
54: *> and the columns of Z are the right Schur vectors.
55: *>
56: *> If only the eigenvalues of (A,B) are needed, the driver routine
57: *> DGEGV should be used instead. See DGEGV for a description of the
58: *> eigenvalues of the generalized nonsymmetric eigenvalue problem
59: *> (GNEP).
60: *> \endverbatim
61: *
62: * Arguments:
63: * ==========
64: *
65: *> \param[in] JOBVSL
66: *> \verbatim
67: *> JOBVSL is CHARACTER*1
68: *> = 'N': do not compute the left Schur vectors;
69: *> = 'V': compute the left Schur vectors (returned in VSL).
70: *> \endverbatim
71: *>
72: *> \param[in] JOBVSR
73: *> \verbatim
74: *> JOBVSR is CHARACTER*1
75: *> = 'N': do not compute the right Schur vectors;
76: *> = 'V': compute the right Schur vectors (returned in VSR).
77: *> \endverbatim
78: *>
79: *> \param[in] N
80: *> \verbatim
81: *> N is INTEGER
82: *> The order of the matrices A, B, VSL, and VSR. N >= 0.
83: *> \endverbatim
84: *>
85: *> \param[in,out] A
86: *> \verbatim
87: *> A is DOUBLE PRECISION array, dimension (LDA, N)
88: *> On entry, the matrix A.
89: *> On exit, the upper quasi-triangular matrix S from the
90: *> generalized real Schur factorization.
91: *> \endverbatim
92: *>
93: *> \param[in] LDA
94: *> \verbatim
95: *> LDA is INTEGER
96: *> The leading dimension of A. LDA >= max(1,N).
97: *> \endverbatim
98: *>
99: *> \param[in,out] B
100: *> \verbatim
101: *> B is DOUBLE PRECISION array, dimension (LDB, N)
102: *> On entry, the matrix B.
103: *> On exit, the upper triangular matrix T from the generalized
104: *> real Schur factorization.
105: *> \endverbatim
106: *>
107: *> \param[in] LDB
108: *> \verbatim
109: *> LDB is INTEGER
110: *> The leading dimension of B. LDB >= max(1,N).
111: *> \endverbatim
112: *>
113: *> \param[out] ALPHAR
114: *> \verbatim
115: *> ALPHAR is DOUBLE PRECISION array, dimension (N)
116: *> The real parts of each scalar alpha defining an eigenvalue
117: *> of GNEP.
118: *> \endverbatim
119: *>
120: *> \param[out] ALPHAI
121: *> \verbatim
122: *> ALPHAI is DOUBLE PRECISION array, dimension (N)
123: *> The imaginary parts of each scalar alpha defining an
124: *> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
125: *> eigenvalue is real; if positive, then the j-th and (j+1)-st
126: *> eigenvalues are a complex conjugate pair, with
127: *> ALPHAI(j+1) = -ALPHAI(j).
128: *> \endverbatim
129: *>
130: *> \param[out] BETA
131: *> \verbatim
132: *> BETA is DOUBLE PRECISION array, dimension (N)
133: *> The scalars beta that define the eigenvalues of GNEP.
134: *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
135: *> beta = BETA(j) represent the j-th eigenvalue of the matrix
136: *> pair (A,B), in one of the forms lambda = alpha/beta or
137: *> mu = beta/alpha. Since either lambda or mu may overflow,
138: *> they should not, in general, be computed.
139: *> \endverbatim
140: *>
141: *> \param[out] VSL
142: *> \verbatim
143: *> VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
144: *> If JOBVSL = 'V', the matrix of left Schur vectors Q.
145: *> Not referenced if JOBVSL = 'N'.
146: *> \endverbatim
147: *>
148: *> \param[in] LDVSL
149: *> \verbatim
150: *> LDVSL is INTEGER
151: *> The leading dimension of the matrix VSL. LDVSL >=1, and
152: *> if JOBVSL = 'V', LDVSL >= N.
153: *> \endverbatim
154: *>
155: *> \param[out] VSR
156: *> \verbatim
157: *> VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
158: *> If JOBVSR = 'V', the matrix of right Schur vectors Z.
159: *> Not referenced if JOBVSR = 'N'.
160: *> \endverbatim
161: *>
162: *> \param[in] LDVSR
163: *> \verbatim
164: *> LDVSR is INTEGER
165: *> The leading dimension of the matrix VSR. LDVSR >= 1, and
166: *> if JOBVSR = 'V', LDVSR >= N.
167: *> \endverbatim
168: *>
169: *> \param[out] WORK
170: *> \verbatim
171: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
172: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
173: *> \endverbatim
174: *>
175: *> \param[in] LWORK
176: *> \verbatim
177: *> LWORK is INTEGER
178: *> The dimension of the array WORK. LWORK >= max(1,4*N).
179: *> For good performance, LWORK must generally be larger.
180: *> To compute the optimal value of LWORK, call ILAENV to get
181: *> blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
182: *> NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
183: *> The optimal LWORK is 2*N + N*(NB+1).
184: *>
185: *> If LWORK = -1, then a workspace query is assumed; the routine
186: *> only calculates the optimal size of the WORK array, returns
187: *> this value as the first entry of the WORK array, and no error
188: *> message related to LWORK is issued by XERBLA.
189: *> \endverbatim
190: *>
191: *> \param[out] INFO
192: *> \verbatim
193: *> INFO is INTEGER
194: *> = 0: successful exit
195: *> < 0: if INFO = -i, the i-th argument had an illegal value.
196: *> = 1,...,N:
197: *> The QZ iteration failed. (A,B) are not in Schur
198: *> form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
199: *> be correct for j=INFO+1,...,N.
200: *> > N: errors that usually indicate LAPACK problems:
201: *> =N+1: error return from DGGBAL
202: *> =N+2: error return from DGEQRF
203: *> =N+3: error return from DORMQR
204: *> =N+4: error return from DORGQR
205: *> =N+5: error return from DGGHRD
206: *> =N+6: error return from DHGEQZ (other than failed
207: *> iteration)
208: *> =N+7: error return from DGGBAK (computing VSL)
209: *> =N+8: error return from DGGBAK (computing VSR)
210: *> =N+9: error return from DLASCL (various places)
211: *> \endverbatim
212: *
213: * Authors:
214: * ========
215: *
216: *> \author Univ. of Tennessee
217: *> \author Univ. of California Berkeley
218: *> \author Univ. of Colorado Denver
219: *> \author NAG Ltd.
220: *
221: *> \ingroup doubleGEeigen
222: *
223: * =====================================================================
224: SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
225: $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
226: $ LWORK, INFO )
227: *
228: * -- LAPACK driver routine --
229: * -- LAPACK is a software package provided by Univ. of Tennessee, --
230: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231: *
232: * .. Scalar Arguments ..
233: CHARACTER JOBVSL, JOBVSR
234: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
235: * ..
236: * .. Array Arguments ..
237: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
238: $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
239: $ VSR( LDVSR, * ), WORK( * )
240: * ..
241: *
242: * =====================================================================
243: *
244: * .. Parameters ..
245: DOUBLE PRECISION ZERO, ONE
246: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
247: * ..
248: * .. Local Scalars ..
249: LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
250: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
251: $ IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
252: $ LWKOPT, NB, NB1, NB2, NB3
253: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254: $ SAFMIN, SMLNUM
255: * ..
256: * .. External Subroutines ..
257: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
258: $ DLASCL, DLASET, DORGQR, DORMQR, XERBLA
259: * ..
260: * .. External Functions ..
261: LOGICAL LSAME
262: INTEGER ILAENV
263: DOUBLE PRECISION DLAMCH, DLANGE
264: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
265: * ..
266: * .. Intrinsic Functions ..
267: INTRINSIC INT, MAX
268: * ..
269: * .. Executable Statements ..
270: *
271: * Decode the input arguments
272: *
273: IF( LSAME( JOBVSL, 'N' ) ) THEN
274: IJOBVL = 1
275: ILVSL = .FALSE.
276: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
277: IJOBVL = 2
278: ILVSL = .TRUE.
279: ELSE
280: IJOBVL = -1
281: ILVSL = .FALSE.
282: END IF
283: *
284: IF( LSAME( JOBVSR, 'N' ) ) THEN
285: IJOBVR = 1
286: ILVSR = .FALSE.
287: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
288: IJOBVR = 2
289: ILVSR = .TRUE.
290: ELSE
291: IJOBVR = -1
292: ILVSR = .FALSE.
293: END IF
294: *
295: * Test the input arguments
296: *
297: LWKMIN = MAX( 4*N, 1 )
298: LWKOPT = LWKMIN
299: WORK( 1 ) = LWKOPT
300: LQUERY = ( LWORK.EQ.-1 )
301: INFO = 0
302: IF( IJOBVL.LE.0 ) THEN
303: INFO = -1
304: ELSE IF( IJOBVR.LE.0 ) THEN
305: INFO = -2
306: ELSE IF( N.LT.0 ) THEN
307: INFO = -3
308: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
309: INFO = -5
310: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
311: INFO = -7
312: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
313: INFO = -12
314: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
315: INFO = -14
316: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
317: INFO = -16
318: END IF
319: *
320: IF( INFO.EQ.0 ) THEN
321: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
322: NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
323: NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
324: NB = MAX( NB1, NB2, NB3 )
325: LOPT = 2*N + N*( NB+1 )
326: WORK( 1 ) = LOPT
327: END IF
328: *
329: IF( INFO.NE.0 ) THEN
330: CALL XERBLA( 'DGEGS ', -INFO )
331: RETURN
332: ELSE IF( LQUERY ) THEN
333: RETURN
334: END IF
335: *
336: * Quick return if possible
337: *
338: IF( N.EQ.0 )
339: $ RETURN
340: *
341: * Get machine constants
342: *
343: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
344: SAFMIN = DLAMCH( 'S' )
345: SMLNUM = N*SAFMIN / 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, WORK )
351: ILASCL = .FALSE.
352: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
353: ANRMTO = SMLNUM
354: ILASCL = .TRUE.
355: ELSE IF( ANRM.GT.BIGNUM ) THEN
356: ANRMTO = BIGNUM
357: ILASCL = .TRUE.
358: END IF
359: *
360: IF( ILASCL ) THEN
361: CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
362: IF( IINFO.NE.0 ) THEN
363: INFO = N + 9
364: RETURN
365: END IF
366: END IF
367: *
368: * Scale B if max element outside range [SMLNUM,BIGNUM]
369: *
370: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
371: ILBSCL = .FALSE.
372: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
373: BNRMTO = SMLNUM
374: ILBSCL = .TRUE.
375: ELSE IF( BNRM.GT.BIGNUM ) THEN
376: BNRMTO = BIGNUM
377: ILBSCL = .TRUE.
378: END IF
379: *
380: IF( ILBSCL ) THEN
381: CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
382: IF( IINFO.NE.0 ) THEN
383: INFO = N + 9
384: RETURN
385: END IF
386: END IF
387: *
388: * Permute the matrix to make it more nearly triangular
389: * Workspace layout: (2*N words -- "work..." not actually used)
390: * left_permutation, right_permutation, work...
391: *
392: ILEFT = 1
393: IRIGHT = N + 1
394: IWORK = IRIGHT + N
395: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
396: $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
397: IF( IINFO.NE.0 ) THEN
398: INFO = N + 1
399: GO TO 10
400: END IF
401: *
402: * Reduce B to triangular form, and initialize VSL and/or VSR
403: * Workspace layout: ("work..." must have at least N words)
404: * left_permutation, right_permutation, tau, work...
405: *
406: IROWS = IHI + 1 - ILO
407: ICOLS = N + 1 - ILO
408: ITAU = IWORK
409: IWORK = ITAU + IROWS
410: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
411: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
412: IF( IINFO.GE.0 )
413: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
414: IF( IINFO.NE.0 ) THEN
415: INFO = N + 2
416: GO TO 10
417: END IF
418: *
419: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
420: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
421: $ LWORK+1-IWORK, IINFO )
422: IF( IINFO.GE.0 )
423: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
424: IF( IINFO.NE.0 ) THEN
425: INFO = N + 3
426: GO TO 10
427: END IF
428: *
429: IF( ILVSL ) THEN
430: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
431: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
432: $ VSL( ILO+1, ILO ), LDVSL )
433: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
434: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
435: $ IINFO )
436: IF( IINFO.GE.0 )
437: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
438: IF( IINFO.NE.0 ) THEN
439: INFO = N + 4
440: GO TO 10
441: END IF
442: END IF
443: *
444: IF( ILVSR )
445: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
446: *
447: * Reduce to generalized Hessenberg form
448: *
449: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
450: $ LDVSL, VSR, LDVSR, IINFO )
451: IF( IINFO.NE.0 ) THEN
452: INFO = N + 5
453: GO TO 10
454: END IF
455: *
456: * Perform QZ algorithm, computing Schur vectors if desired
457: * Workspace layout: ("work..." must have at least 1 word)
458: * left_permutation, right_permutation, work...
459: *
460: IWORK = ITAU
461: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
462: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
463: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
464: IF( IINFO.GE.0 )
465: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
466: IF( IINFO.NE.0 ) THEN
467: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
468: INFO = IINFO
469: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
470: INFO = IINFO - N
471: ELSE
472: INFO = N + 6
473: END IF
474: GO TO 10
475: END IF
476: *
477: * Apply permutation to VSL and VSR
478: *
479: IF( ILVSL ) THEN
480: CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
481: $ WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
482: IF( IINFO.NE.0 ) THEN
483: INFO = N + 7
484: GO TO 10
485: END IF
486: END IF
487: IF( ILVSR ) THEN
488: CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
489: $ WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
490: IF( IINFO.NE.0 ) THEN
491: INFO = N + 8
492: GO TO 10
493: END IF
494: END IF
495: *
496: * Undo scaling
497: *
498: IF( ILASCL ) THEN
499: CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
500: IF( IINFO.NE.0 ) THEN
501: INFO = N + 9
502: RETURN
503: END IF
504: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
505: $ IINFO )
506: IF( IINFO.NE.0 ) THEN
507: INFO = N + 9
508: RETURN
509: END IF
510: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
511: $ IINFO )
512: IF( IINFO.NE.0 ) THEN
513: INFO = N + 9
514: RETURN
515: END IF
516: END IF
517: *
518: IF( ILBSCL ) THEN
519: CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
520: IF( IINFO.NE.0 ) THEN
521: INFO = N + 9
522: RETURN
523: END IF
524: CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
525: IF( IINFO.NE.0 ) THEN
526: INFO = N + 9
527: RETURN
528: END IF
529: END IF
530: *
531: 10 CONTINUE
532: WORK( 1 ) = LWKOPT
533: *
534: RETURN
535: *
536: * End of DGEGS
537: *
538: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgegs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack