Annotation of rpl/lapack/lapack/dgegs.f, revision 1.10
1.8 bertrand 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: *> \date November 2011
222: *
223: *> \ingroup doubleGEeigen
224: *
225: * =====================================================================
1.1 bertrand 226: SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
227: $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
228: $ LWORK, INFO )
229: *
1.8 bertrand 230: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 231: * -- LAPACK is a software package provided by Univ. of Tennessee, --
232: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 bertrand 233: * November 2011
1.1 bertrand 234: *
235: * .. Scalar Arguments ..
236: CHARACTER JOBVSL, JOBVSR
237: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
238: * ..
239: * .. Array Arguments ..
240: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
241: $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
242: $ VSR( LDVSR, * ), WORK( * )
243: * ..
244: *
245: * =====================================================================
246: *
247: * .. Parameters ..
248: DOUBLE PRECISION ZERO, ONE
249: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
250: * ..
251: * .. Local Scalars ..
252: LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
253: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
254: $ IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
255: $ LWKOPT, NB, NB1, NB2, NB3
256: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
257: $ SAFMIN, SMLNUM
258: * ..
259: * .. External Subroutines ..
260: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
261: $ DLASCL, DLASET, DORGQR, DORMQR, XERBLA
262: * ..
263: * .. External Functions ..
264: LOGICAL LSAME
265: INTEGER ILAENV
266: DOUBLE PRECISION DLAMCH, DLANGE
267: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
268: * ..
269: * .. Intrinsic Functions ..
270: INTRINSIC INT, MAX
271: * ..
272: * .. Executable Statements ..
273: *
274: * Decode the input arguments
275: *
276: IF( LSAME( JOBVSL, 'N' ) ) THEN
277: IJOBVL = 1
278: ILVSL = .FALSE.
279: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
280: IJOBVL = 2
281: ILVSL = .TRUE.
282: ELSE
283: IJOBVL = -1
284: ILVSL = .FALSE.
285: END IF
286: *
287: IF( LSAME( JOBVSR, 'N' ) ) THEN
288: IJOBVR = 1
289: ILVSR = .FALSE.
290: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
291: IJOBVR = 2
292: ILVSR = .TRUE.
293: ELSE
294: IJOBVR = -1
295: ILVSR = .FALSE.
296: END IF
297: *
298: * Test the input arguments
299: *
300: LWKMIN = MAX( 4*N, 1 )
301: LWKOPT = LWKMIN
302: WORK( 1 ) = LWKOPT
303: LQUERY = ( LWORK.EQ.-1 )
304: INFO = 0
305: IF( IJOBVL.LE.0 ) THEN
306: INFO = -1
307: ELSE IF( IJOBVR.LE.0 ) THEN
308: INFO = -2
309: ELSE IF( N.LT.0 ) THEN
310: INFO = -3
311: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
312: INFO = -5
313: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
314: INFO = -7
315: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
316: INFO = -12
317: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
318: INFO = -14
319: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
320: INFO = -16
321: END IF
322: *
323: IF( INFO.EQ.0 ) THEN
324: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
325: NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
326: NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
327: NB = MAX( NB1, NB2, NB3 )
328: LOPT = 2*N + N*( NB+1 )
329: WORK( 1 ) = LOPT
330: END IF
331: *
332: IF( INFO.NE.0 ) THEN
333: CALL XERBLA( 'DGEGS ', -INFO )
334: RETURN
335: ELSE IF( LQUERY ) THEN
336: RETURN
337: END IF
338: *
339: * Quick return if possible
340: *
341: IF( N.EQ.0 )
342: $ RETURN
343: *
344: * Get machine constants
345: *
346: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
347: SAFMIN = DLAMCH( 'S' )
348: SMLNUM = N*SAFMIN / EPS
349: BIGNUM = ONE / SMLNUM
350: *
351: * Scale A if max element outside range [SMLNUM,BIGNUM]
352: *
353: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
354: ILASCL = .FALSE.
355: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
356: ANRMTO = SMLNUM
357: ILASCL = .TRUE.
358: ELSE IF( ANRM.GT.BIGNUM ) THEN
359: ANRMTO = BIGNUM
360: ILASCL = .TRUE.
361: END IF
362: *
363: IF( ILASCL ) THEN
364: CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
365: IF( IINFO.NE.0 ) THEN
366: INFO = N + 9
367: RETURN
368: END IF
369: END IF
370: *
371: * Scale B if max element outside range [SMLNUM,BIGNUM]
372: *
373: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
374: ILBSCL = .FALSE.
375: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
376: BNRMTO = SMLNUM
377: ILBSCL = .TRUE.
378: ELSE IF( BNRM.GT.BIGNUM ) THEN
379: BNRMTO = BIGNUM
380: ILBSCL = .TRUE.
381: END IF
382: *
383: IF( ILBSCL ) THEN
384: CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
385: IF( IINFO.NE.0 ) THEN
386: INFO = N + 9
387: RETURN
388: END IF
389: END IF
390: *
391: * Permute the matrix to make it more nearly triangular
392: * Workspace layout: (2*N words -- "work..." not actually used)
393: * left_permutation, right_permutation, work...
394: *
395: ILEFT = 1
396: IRIGHT = N + 1
397: IWORK = IRIGHT + N
398: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
399: $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
400: IF( IINFO.NE.0 ) THEN
401: INFO = N + 1
402: GO TO 10
403: END IF
404: *
405: * Reduce B to triangular form, and initialize VSL and/or VSR
406: * Workspace layout: ("work..." must have at least N words)
407: * left_permutation, right_permutation, tau, work...
408: *
409: IROWS = IHI + 1 - ILO
410: ICOLS = N + 1 - ILO
411: ITAU = IWORK
412: IWORK = ITAU + IROWS
413: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
414: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
415: IF( IINFO.GE.0 )
416: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
417: IF( IINFO.NE.0 ) THEN
418: INFO = N + 2
419: GO TO 10
420: END IF
421: *
422: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
423: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
424: $ LWORK+1-IWORK, IINFO )
425: IF( IINFO.GE.0 )
426: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
427: IF( IINFO.NE.0 ) THEN
428: INFO = N + 3
429: GO TO 10
430: END IF
431: *
432: IF( ILVSL ) THEN
433: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
434: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
435: $ VSL( ILO+1, ILO ), LDVSL )
436: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
437: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
438: $ IINFO )
439: IF( IINFO.GE.0 )
440: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
441: IF( IINFO.NE.0 ) THEN
442: INFO = N + 4
443: GO TO 10
444: END IF
445: END IF
446: *
447: IF( ILVSR )
448: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
449: *
450: * Reduce to generalized Hessenberg form
451: *
452: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
453: $ LDVSL, VSR, LDVSR, IINFO )
454: IF( IINFO.NE.0 ) THEN
455: INFO = N + 5
456: GO TO 10
457: END IF
458: *
459: * Perform QZ algorithm, computing Schur vectors if desired
460: * Workspace layout: ("work..." must have at least 1 word)
461: * left_permutation, right_permutation, work...
462: *
463: IWORK = ITAU
464: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
465: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
466: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
467: IF( IINFO.GE.0 )
468: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
469: IF( IINFO.NE.0 ) THEN
470: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
471: INFO = IINFO
472: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
473: INFO = IINFO - N
474: ELSE
475: INFO = N + 6
476: END IF
477: GO TO 10
478: END IF
479: *
480: * Apply permutation to VSL and VSR
481: *
482: IF( ILVSL ) THEN
483: CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
484: $ WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
485: IF( IINFO.NE.0 ) THEN
486: INFO = N + 7
487: GO TO 10
488: END IF
489: END IF
490: IF( ILVSR ) THEN
491: CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
492: $ WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
493: IF( IINFO.NE.0 ) THEN
494: INFO = N + 8
495: GO TO 10
496: END IF
497: END IF
498: *
499: * Undo scaling
500: *
501: IF( ILASCL ) THEN
502: CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
503: IF( IINFO.NE.0 ) THEN
504: INFO = N + 9
505: RETURN
506: END IF
507: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
508: $ IINFO )
509: IF( IINFO.NE.0 ) THEN
510: INFO = N + 9
511: RETURN
512: END IF
513: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
514: $ IINFO )
515: IF( IINFO.NE.0 ) THEN
516: INFO = N + 9
517: RETURN
518: END IF
519: END IF
520: *
521: IF( ILBSCL ) THEN
522: CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
523: IF( IINFO.NE.0 ) THEN
524: INFO = N + 9
525: RETURN
526: END IF
527: CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
528: IF( IINFO.NE.0 ) THEN
529: INFO = N + 9
530: RETURN
531: END IF
532: END IF
533: *
534: 10 CONTINUE
535: WORK( 1 ) = LWKOPT
536: *
537: RETURN
538: *
539: * End of DGEGS
540: *
541: 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