Annotation of rpl/lapack/lapack/zgegs.f, revision 1.2
1.1 bertrand 1: SUBROUTINE ZGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
2: $ VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
3: $ INFO )
4: *
5: * -- LAPACK driver routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER JOBVSL, JOBVSR
12: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION RWORK( * )
16: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
17: $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
18: $ WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * This routine is deprecated and has been replaced by routine ZGGES.
25: *
26: * ZGEGS computes the eigenvalues, Schur form, and, optionally, the
27: * left and or/right Schur vectors of a complex matrix pair (A,B).
28: * Given two square matrices A and B, the generalized Schur
29: * factorization has the form
30: *
31: * A = Q*S*Z**H, B = Q*T*Z**H
32: *
33: * where Q and Z are unitary matrices and S and T are upper triangular.
34: * The columns of Q are the left Schur vectors
35: * and the columns of Z are the right Schur vectors.
36: *
37: * If only the eigenvalues of (A,B) are needed, the driver routine
38: * ZGEGV should be used instead. See ZGEGV for a description of the
39: * eigenvalues of the generalized nonsymmetric eigenvalue problem
40: * (GNEP).
41: *
42: * Arguments
43: * =========
44: *
45: * JOBVSL (input) CHARACTER*1
46: * = 'N': do not compute the left Schur vectors;
47: * = 'V': compute the left Schur vectors (returned in VSL).
48: *
49: * JOBVSR (input) CHARACTER*1
50: * = 'N': do not compute the right Schur vectors;
51: * = 'V': compute the right Schur vectors (returned in VSR).
52: *
53: * N (input) INTEGER
54: * The order of the matrices A, B, VSL, and VSR. N >= 0.
55: *
56: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
57: * On entry, the matrix A.
58: * On exit, the upper triangular matrix S from the generalized
59: * Schur factorization.
60: *
61: * LDA (input) INTEGER
62: * The leading dimension of A. LDA >= max(1,N).
63: *
64: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
65: * On entry, the matrix B.
66: * On exit, the upper triangular matrix T from the generalized
67: * Schur factorization.
68: *
69: * LDB (input) INTEGER
70: * The leading dimension of B. LDB >= max(1,N).
71: *
72: * ALPHA (output) COMPLEX*16 array, dimension (N)
73: * The complex scalars alpha that define the eigenvalues of
74: * GNEP. ALPHA(j) = S(j,j), the diagonal element of the Schur
75: * form of A.
76: *
77: * BETA (output) COMPLEX*16 array, dimension (N)
78: * The non-negative real scalars beta that define the
79: * eigenvalues of GNEP. BETA(j) = T(j,j), the diagonal element
80: * of the triangular factor T.
81: *
82: * Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
83: * represent the j-th eigenvalue of the matrix pair (A,B), in
84: * one of the forms lambda = alpha/beta or mu = beta/alpha.
85: * Since either lambda or mu may overflow, they should not,
86: * in general, be computed.
87: *
88: *
89: * VSL (output) COMPLEX*16 array, dimension (LDVSL,N)
90: * If JOBVSL = 'V', the matrix of left Schur vectors Q.
91: * Not referenced if JOBVSL = 'N'.
92: *
93: * LDVSL (input) INTEGER
94: * The leading dimension of the matrix VSL. LDVSL >= 1, and
95: * if JOBVSL = 'V', LDVSL >= N.
96: *
97: * VSR (output) COMPLEX*16 array, dimension (LDVSR,N)
98: * If JOBVSR = 'V', the matrix of right Schur vectors Z.
99: * Not referenced if JOBVSR = 'N'.
100: *
101: * LDVSR (input) INTEGER
102: * The leading dimension of the matrix VSR. LDVSR >= 1, and
103: * if JOBVSR = 'V', LDVSR >= N.
104: *
105: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
106: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
107: *
108: * LWORK (input) INTEGER
109: * The dimension of the array WORK. LWORK >= max(1,2*N).
110: * For good performance, LWORK must generally be larger.
111: * To compute the optimal value of LWORK, call ILAENV to get
112: * blocksizes (for ZGEQRF, ZUNMQR, and CUNGQR.) Then compute:
113: * NB -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and CUNGQR;
114: * the optimal LWORK is N*(NB+1).
115: *
116: * If LWORK = -1, then a workspace query is assumed; the routine
117: * only calculates the optimal size of the WORK array, returns
118: * this value as the first entry of the WORK array, and no error
119: * message related to LWORK is issued by XERBLA.
120: *
121: * RWORK (workspace) DOUBLE PRECISION array, dimension (3*N)
122: *
123: * INFO (output) INTEGER
124: * = 0: successful exit
125: * < 0: if INFO = -i, the i-th argument had an illegal value.
126: * =1,...,N:
127: * The QZ iteration failed. (A,B) are not in Schur
128: * form, but ALPHA(j) and BETA(j) should be correct for
129: * j=INFO+1,...,N.
130: * > N: errors that usually indicate LAPACK problems:
131: * =N+1: error return from ZGGBAL
132: * =N+2: error return from ZGEQRF
133: * =N+3: error return from ZUNMQR
134: * =N+4: error return from ZUNGQR
135: * =N+5: error return from ZGGHRD
136: * =N+6: error return from ZHGEQZ (other than failed
137: * iteration)
138: * =N+7: error return from ZGGBAK (computing VSL)
139: * =N+8: error return from ZGGBAK (computing VSR)
140: * =N+9: error return from ZLASCL (various places)
141: *
142: * =====================================================================
143: *
144: * .. Parameters ..
145: DOUBLE PRECISION ZERO, ONE
146: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
147: COMPLEX*16 CZERO, CONE
148: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
149: $ CONE = ( 1.0D0, 0.0D0 ) )
150: * ..
151: * .. Local Scalars ..
152: LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
153: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
154: $ IRIGHT, IROWS, IRWORK, ITAU, IWORK, LOPT,
155: $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
156: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
157: $ SAFMIN, SMLNUM
158: * ..
159: * .. External Subroutines ..
160: EXTERNAL XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
161: $ ZLACPY, ZLASCL, ZLASET, ZUNGQR, ZUNMQR
162: * ..
163: * .. External Functions ..
164: LOGICAL LSAME
165: INTEGER ILAENV
166: DOUBLE PRECISION DLAMCH, ZLANGE
167: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
168: * ..
169: * .. Intrinsic Functions ..
170: INTRINSIC INT, MAX
171: * ..
172: * .. Executable Statements ..
173: *
174: * Decode the input arguments
175: *
176: IF( LSAME( JOBVSL, 'N' ) ) THEN
177: IJOBVL = 1
178: ILVSL = .FALSE.
179: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
180: IJOBVL = 2
181: ILVSL = .TRUE.
182: ELSE
183: IJOBVL = -1
184: ILVSL = .FALSE.
185: END IF
186: *
187: IF( LSAME( JOBVSR, 'N' ) ) THEN
188: IJOBVR = 1
189: ILVSR = .FALSE.
190: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
191: IJOBVR = 2
192: ILVSR = .TRUE.
193: ELSE
194: IJOBVR = -1
195: ILVSR = .FALSE.
196: END IF
197: *
198: * Test the input arguments
199: *
200: LWKMIN = MAX( 2*N, 1 )
201: LWKOPT = LWKMIN
202: WORK( 1 ) = LWKOPT
203: LQUERY = ( LWORK.EQ.-1 )
204: INFO = 0
205: IF( IJOBVL.LE.0 ) THEN
206: INFO = -1
207: ELSE IF( IJOBVR.LE.0 ) THEN
208: INFO = -2
209: ELSE IF( N.LT.0 ) THEN
210: INFO = -3
211: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
212: INFO = -5
213: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
214: INFO = -7
215: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
216: INFO = -11
217: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
218: INFO = -13
219: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
220: INFO = -15
221: END IF
222: *
223: IF( INFO.EQ.0 ) THEN
224: NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
225: NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
226: NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
227: NB = MAX( NB1, NB2, NB3 )
228: LOPT = N*( NB+1 )
229: WORK( 1 ) = LOPT
230: END IF
231: *
232: IF( INFO.NE.0 ) THEN
233: CALL XERBLA( 'ZGEGS ', -INFO )
234: RETURN
235: ELSE IF( LQUERY ) THEN
236: RETURN
237: END IF
238: *
239: * Quick return if possible
240: *
241: IF( N.EQ.0 )
242: $ RETURN
243: *
244: * Get machine constants
245: *
246: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
247: SAFMIN = DLAMCH( 'S' )
248: SMLNUM = N*SAFMIN / EPS
249: BIGNUM = ONE / SMLNUM
250: *
251: * Scale A if max element outside range [SMLNUM,BIGNUM]
252: *
253: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
254: ILASCL = .FALSE.
255: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
256: ANRMTO = SMLNUM
257: ILASCL = .TRUE.
258: ELSE IF( ANRM.GT.BIGNUM ) THEN
259: ANRMTO = BIGNUM
260: ILASCL = .TRUE.
261: END IF
262: *
263: IF( ILASCL ) THEN
264: CALL ZLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
265: IF( IINFO.NE.0 ) THEN
266: INFO = N + 9
267: RETURN
268: END IF
269: END IF
270: *
271: * Scale B if max element outside range [SMLNUM,BIGNUM]
272: *
273: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
274: ILBSCL = .FALSE.
275: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
276: BNRMTO = SMLNUM
277: ILBSCL = .TRUE.
278: ELSE IF( BNRM.GT.BIGNUM ) THEN
279: BNRMTO = BIGNUM
280: ILBSCL = .TRUE.
281: END IF
282: *
283: IF( ILBSCL ) THEN
284: CALL ZLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
285: IF( IINFO.NE.0 ) THEN
286: INFO = N + 9
287: RETURN
288: END IF
289: END IF
290: *
291: * Permute the matrix to make it more nearly triangular
292: *
293: ILEFT = 1
294: IRIGHT = N + 1
295: IRWORK = IRIGHT + N
296: IWORK = 1
297: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
298: $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
299: IF( IINFO.NE.0 ) THEN
300: INFO = N + 1
301: GO TO 10
302: END IF
303: *
304: * Reduce B to triangular form, and initialize VSL and/or VSR
305: *
306: IROWS = IHI + 1 - ILO
307: ICOLS = N + 1 - ILO
308: ITAU = IWORK
309: IWORK = ITAU + IROWS
310: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
311: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
312: IF( IINFO.GE.0 )
313: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
314: IF( IINFO.NE.0 ) THEN
315: INFO = N + 2
316: GO TO 10
317: END IF
318: *
319: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
320: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
321: $ LWORK+1-IWORK, IINFO )
322: IF( IINFO.GE.0 )
323: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
324: IF( IINFO.NE.0 ) THEN
325: INFO = N + 3
326: GO TO 10
327: END IF
328: *
329: IF( ILVSL ) THEN
330: CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
331: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
332: $ VSL( ILO+1, ILO ), LDVSL )
333: CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
334: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
335: $ IINFO )
336: IF( IINFO.GE.0 )
337: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
338: IF( IINFO.NE.0 ) THEN
339: INFO = N + 4
340: GO TO 10
341: END IF
342: END IF
343: *
344: IF( ILVSR )
345: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
346: *
347: * Reduce to generalized Hessenberg form
348: *
349: CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
350: $ LDVSL, VSR, LDVSR, IINFO )
351: IF( IINFO.NE.0 ) THEN
352: INFO = N + 5
353: GO TO 10
354: END IF
355: *
356: * Perform QZ algorithm, computing Schur vectors if desired
357: *
358: IWORK = ITAU
359: CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
360: $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWORK ),
361: $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
362: IF( IINFO.GE.0 )
363: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
364: IF( IINFO.NE.0 ) THEN
365: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
366: INFO = IINFO
367: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
368: INFO = IINFO - N
369: ELSE
370: INFO = N + 6
371: END IF
372: GO TO 10
373: END IF
374: *
375: * Apply permutation to VSL and VSR
376: *
377: IF( ILVSL ) THEN
378: CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
379: $ RWORK( IRIGHT ), N, VSL, LDVSL, IINFO )
380: IF( IINFO.NE.0 ) THEN
381: INFO = N + 7
382: GO TO 10
383: END IF
384: END IF
385: IF( ILVSR ) THEN
386: CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
387: $ RWORK( IRIGHT ), N, VSR, LDVSR, IINFO )
388: IF( IINFO.NE.0 ) THEN
389: INFO = N + 8
390: GO TO 10
391: END IF
392: END IF
393: *
394: * Undo scaling
395: *
396: IF( ILASCL ) THEN
397: CALL ZLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
398: IF( IINFO.NE.0 ) THEN
399: INFO = N + 9
400: RETURN
401: END IF
402: CALL ZLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
403: IF( IINFO.NE.0 ) THEN
404: INFO = N + 9
405: RETURN
406: END IF
407: END IF
408: *
409: IF( ILBSCL ) THEN
410: CALL ZLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
411: IF( IINFO.NE.0 ) THEN
412: INFO = N + 9
413: RETURN
414: END IF
415: CALL ZLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
416: IF( IINFO.NE.0 ) THEN
417: INFO = N + 9
418: RETURN
419: END IF
420: END IF
421: *
422: 10 CONTINUE
423: WORK( 1 ) = LWKOPT
424: *
425: RETURN
426: *
427: * End of ZGEGS
428: *
429: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgegs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack