1: SUBROUTINE ZGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
2: $ ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI,
3: $ LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV,
4: $ WORK, LWORK, RWORK, IWORK, BWORK, INFO )
5: *
6: * -- LAPACK driver routine (version 3.2) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * November 2006
10: *
11: * .. Scalar Arguments ..
12: CHARACTER BALANC, JOBVL, JOBVR, SENSE
13: INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
14: DOUBLE PRECISION ABNRM, BBNRM
15: * ..
16: * .. Array Arguments ..
17: LOGICAL BWORK( * )
18: INTEGER IWORK( * )
19: DOUBLE PRECISION LSCALE( * ), RCONDE( * ), RCONDV( * ),
20: $ RSCALE( * ), RWORK( * )
21: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
22: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
23: $ WORK( * )
24: * ..
25: *
26: * Purpose
27: * =======
28: *
29: * ZGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
30: * (A,B) the generalized eigenvalues, and optionally, the left and/or
31: * right generalized eigenvectors.
32: *
33: * Optionally, it also computes a balancing transformation to improve
34: * the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
35: * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
36: * the eigenvalues (RCONDE), and reciprocal condition numbers for the
37: * right eigenvectors (RCONDV).
38: *
39: * A generalized eigenvalue for a pair of matrices (A,B) is a scalar
40: * lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
41: * singular. It is usually represented as the pair (alpha,beta), as
42: * there is a reasonable interpretation for beta=0, and even for both
43: * being zero.
44: *
45: * The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
46: * of (A,B) satisfies
47: * A * v(j) = lambda(j) * B * v(j) .
48: * The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
49: * of (A,B) satisfies
50: * u(j)**H * A = lambda(j) * u(j)**H * B.
51: * where u(j)**H is the conjugate-transpose of u(j).
52: *
53: *
54: * Arguments
55: * =========
56: *
57: * BALANC (input) CHARACTER*1
58: * Specifies the balance option to be performed:
59: * = 'N': do not diagonally scale or permute;
60: * = 'P': permute only;
61: * = 'S': scale only;
62: * = 'B': both permute and scale.
63: * Computed reciprocal condition numbers will be for the
64: * matrices after permuting and/or balancing. Permuting does
65: * not change condition numbers (in exact arithmetic), but
66: * balancing does.
67: *
68: * JOBVL (input) CHARACTER*1
69: * = 'N': do not compute the left generalized eigenvectors;
70: * = 'V': compute the left generalized eigenvectors.
71: *
72: * JOBVR (input) CHARACTER*1
73: * = 'N': do not compute the right generalized eigenvectors;
74: * = 'V': compute the right generalized eigenvectors.
75: *
76: * SENSE (input) CHARACTER*1
77: * Determines which reciprocal condition numbers are computed.
78: * = 'N': none are computed;
79: * = 'E': computed for eigenvalues only;
80: * = 'V': computed for eigenvectors only;
81: * = 'B': computed for eigenvalues and eigenvectors.
82: *
83: * N (input) INTEGER
84: * The order of the matrices A, B, VL, and VR. N >= 0.
85: *
86: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
87: * On entry, the matrix A in the pair (A,B).
88: * On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
89: * or both, then A contains the first part of the complex Schur
90: * form of the "balanced" versions of the input A and B.
91: *
92: * LDA (input) INTEGER
93: * The leading dimension of A. LDA >= max(1,N).
94: *
95: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
96: * On entry, the matrix B in the pair (A,B).
97: * On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
98: * or both, then B contains the second part of the complex
99: * Schur form of the "balanced" versions of the input A and B.
100: *
101: * LDB (input) INTEGER
102: * The leading dimension of B. LDB >= max(1,N).
103: *
104: * ALPHA (output) COMPLEX*16 array, dimension (N)
105: * BETA (output) COMPLEX*16 array, dimension (N)
106: * On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the generalized
107: * eigenvalues.
108: *
109: * Note: the quotient ALPHA(j)/BETA(j) ) may easily over- or
110: * underflow, and BETA(j) may even be zero. Thus, the user
111: * should avoid naively computing the ratio ALPHA/BETA.
112: * However, ALPHA will be always less than and usually
113: * comparable with norm(A) in magnitude, and BETA always less
114: * than and usually comparable with norm(B).
115: *
116: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
117: * If JOBVL = 'V', the left generalized eigenvectors u(j) are
118: * stored one after another in the columns of VL, in the same
119: * order as their eigenvalues.
120: * Each eigenvector will be scaled so the largest component
121: * will have abs(real part) + abs(imag. part) = 1.
122: * Not referenced if JOBVL = 'N'.
123: *
124: * LDVL (input) INTEGER
125: * The leading dimension of the matrix VL. LDVL >= 1, and
126: * if JOBVL = 'V', LDVL >= N.
127: *
128: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
129: * If JOBVR = 'V', the right generalized eigenvectors v(j) are
130: * stored one after another in the columns of VR, in the same
131: * order as their eigenvalues.
132: * Each eigenvector will be scaled so the largest component
133: * will have abs(real part) + abs(imag. part) = 1.
134: * Not referenced if JOBVR = 'N'.
135: *
136: * LDVR (input) INTEGER
137: * The leading dimension of the matrix VR. LDVR >= 1, and
138: * if JOBVR = 'V', LDVR >= N.
139: *
140: * ILO (output) INTEGER
141: * IHI (output) INTEGER
142: * ILO and IHI are integer values such that on exit
143: * A(i,j) = 0 and B(i,j) = 0 if i > j and
144: * j = 1,...,ILO-1 or i = IHI+1,...,N.
145: * If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
146: *
147: * LSCALE (output) DOUBLE PRECISION array, dimension (N)
148: * Details of the permutations and scaling factors applied
149: * to the left side of A and B. If PL(j) is the index of the
150: * row interchanged with row j, and DL(j) is the scaling
151: * factor applied to row j, then
152: * LSCALE(j) = PL(j) for j = 1,...,ILO-1
153: * = DL(j) for j = ILO,...,IHI
154: * = PL(j) for j = IHI+1,...,N.
155: * The order in which the interchanges are made is N to IHI+1,
156: * then 1 to ILO-1.
157: *
158: * RSCALE (output) DOUBLE PRECISION array, dimension (N)
159: * Details of the permutations and scaling factors applied
160: * to the right side of A and B. If PR(j) is the index of the
161: * column interchanged with column j, and DR(j) is the scaling
162: * factor applied to column j, then
163: * RSCALE(j) = PR(j) for j = 1,...,ILO-1
164: * = DR(j) for j = ILO,...,IHI
165: * = PR(j) for j = IHI+1,...,N
166: * The order in which the interchanges are made is N to IHI+1,
167: * then 1 to ILO-1.
168: *
169: * ABNRM (output) DOUBLE PRECISION
170: * The one-norm of the balanced matrix A.
171: *
172: * BBNRM (output) DOUBLE PRECISION
173: * The one-norm of the balanced matrix B.
174: *
175: * RCONDE (output) DOUBLE PRECISION array, dimension (N)
176: * If SENSE = 'E' or 'B', the reciprocal condition numbers of
177: * the eigenvalues, stored in consecutive elements of the array.
178: * If SENSE = 'N' or 'V', RCONDE is not referenced.
179: *
180: * RCONDV (output) DOUBLE PRECISION array, dimension (N)
181: * If JOB = 'V' or 'B', the estimated reciprocal condition
182: * numbers of the eigenvectors, stored in consecutive elements
183: * of the array. If the eigenvalues cannot be reordered to
184: * compute RCONDV(j), RCONDV(j) is set to 0; this can only occur
185: * when the true value would be very small anyway.
186: * If SENSE = 'N' or 'E', RCONDV is not referenced.
187: *
188: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
189: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
190: *
191: * LWORK (input) INTEGER
192: * The dimension of the array WORK. LWORK >= max(1,2*N).
193: * If SENSE = 'E', LWORK >= max(1,4*N).
194: * If SENSE = 'V' or 'B', LWORK >= max(1,2*N*N+2*N).
195: *
196: * If LWORK = -1, then a workspace query is assumed; the routine
197: * only calculates the optimal size of the WORK array, returns
198: * this value as the first entry of the WORK array, and no error
199: * message related to LWORK is issued by XERBLA.
200: *
201: * RWORK (workspace) REAL array, dimension (lrwork)
202: * lrwork must be at least max(1,6*N) if BALANC = 'S' or 'B',
203: * and at least max(1,2*N) otherwise.
204: * Real workspace.
205: *
206: * IWORK (workspace) INTEGER array, dimension (N+2)
207: * If SENSE = 'E', IWORK is not referenced.
208: *
209: * BWORK (workspace) LOGICAL array, dimension (N)
210: * If SENSE = 'N', BWORK is not referenced.
211: *
212: * INFO (output) INTEGER
213: * = 0: successful exit
214: * < 0: if INFO = -i, the i-th argument had an illegal value.
215: * = 1,...,N:
216: * The QZ iteration failed. No eigenvectors have been
217: * calculated, but ALPHA(j) and BETA(j) should be correct
218: * for j=INFO+1,...,N.
219: * > N: =N+1: other than QZ iteration failed in ZHGEQZ.
220: * =N+2: error return from ZTGEVC.
221: *
222: * Further Details
223: * ===============
224: *
225: * Balancing a matrix pair (A,B) includes, first, permuting rows and
226: * columns to isolate eigenvalues, second, applying diagonal similarity
227: * transformation to the rows and columns to make the rows and columns
228: * as close in norm as possible. The computed reciprocal condition
229: * numbers correspond to the balanced matrix. Permuting rows and columns
230: * will not change the condition numbers (in exact arithmetic) but
231: * diagonal scaling will. For further explanation of balancing, see
232: * section 4.11.1.2 of LAPACK Users' Guide.
233: *
234: * An approximate error bound on the chordal distance between the i-th
235: * computed generalized eigenvalue w and the corresponding exact
236: * eigenvalue lambda is
237: *
238: * chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
239: *
240: * An approximate error bound for the angle between the i-th computed
241: * eigenvector VL(i) or VR(i) is given by
242: *
243: * EPS * norm(ABNRM, BBNRM) / DIF(i).
244: *
245: * For further explanation of the reciprocal condition numbers RCONDE
246: * and RCONDV, see section 4.11 of LAPACK User's Guide.
247: *
248: * .. Parameters ..
249: DOUBLE PRECISION ZERO, ONE
250: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
251: COMPLEX*16 CZERO, CONE
252: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
253: $ CONE = ( 1.0D+0, 0.0D+0 ) )
254: * ..
255: * .. Local Scalars ..
256: LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
257: $ WANTSB, WANTSE, WANTSN, WANTSV
258: CHARACTER CHTEMP
259: INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
260: $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK, MINWRK
261: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
262: $ SMLNUM, TEMP
263: COMPLEX*16 X
264: * ..
265: * .. Local Arrays ..
266: LOGICAL LDUMMA( 1 )
267: * ..
268: * .. External Subroutines ..
269: EXTERNAL DLABAD, DLASCL, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL,
270: $ ZGGHRD, ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC,
271: $ ZTGSNA, ZUNGQR, ZUNMQR
272: * ..
273: * .. External Functions ..
274: LOGICAL LSAME
275: INTEGER ILAENV
276: DOUBLE PRECISION DLAMCH, ZLANGE
277: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
278: * ..
279: * .. Intrinsic Functions ..
280: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
281: * ..
282: * .. Statement Functions ..
283: DOUBLE PRECISION ABS1
284: * ..
285: * .. Statement Function definitions ..
286: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
287: * ..
288: * .. Executable Statements ..
289: *
290: * Decode the input arguments
291: *
292: IF( LSAME( JOBVL, 'N' ) ) THEN
293: IJOBVL = 1
294: ILVL = .FALSE.
295: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
296: IJOBVL = 2
297: ILVL = .TRUE.
298: ELSE
299: IJOBVL = -1
300: ILVL = .FALSE.
301: END IF
302: *
303: IF( LSAME( JOBVR, 'N' ) ) THEN
304: IJOBVR = 1
305: ILVR = .FALSE.
306: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
307: IJOBVR = 2
308: ILVR = .TRUE.
309: ELSE
310: IJOBVR = -1
311: ILVR = .FALSE.
312: END IF
313: ILV = ILVL .OR. ILVR
314: *
315: NOSCL = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
316: WANTSN = LSAME( SENSE, 'N' )
317: WANTSE = LSAME( SENSE, 'E' )
318: WANTSV = LSAME( SENSE, 'V' )
319: WANTSB = LSAME( SENSE, 'B' )
320: *
321: * Test the input arguments
322: *
323: INFO = 0
324: LQUERY = ( LWORK.EQ.-1 )
325: IF( .NOT.( NOSCL .OR. LSAME( BALANC,'S' ) .OR.
326: $ LSAME( BALANC, 'B' ) ) ) THEN
327: INFO = -1
328: ELSE IF( IJOBVL.LE.0 ) THEN
329: INFO = -2
330: ELSE IF( IJOBVR.LE.0 ) THEN
331: INFO = -3
332: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
333: $ THEN
334: INFO = -4
335: ELSE IF( N.LT.0 ) THEN
336: INFO = -5
337: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
338: INFO = -7
339: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
340: INFO = -9
341: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
342: INFO = -13
343: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
344: INFO = -15
345: END IF
346: *
347: * Compute workspace
348: * (Note: Comments in the code beginning "Workspace:" describe the
349: * minimal amount of workspace needed at that point in the code,
350: * as well as the preferred amount for good performance.
351: * NB refers to the optimal block size for the immediately
352: * following subroutine, as returned by ILAENV. The workspace is
353: * computed assuming ILO = 1 and IHI = N, the worst case.)
354: *
355: IF( INFO.EQ.0 ) THEN
356: IF( N.EQ.0 ) THEN
357: MINWRK = 1
358: MAXWRK = 1
359: ELSE
360: MINWRK = 2*N
361: IF( WANTSE ) THEN
362: MINWRK = 4*N
363: ELSE IF( WANTSV .OR. WANTSB ) THEN
364: MINWRK = 2*N*( N + 1)
365: END IF
366: MAXWRK = MINWRK
367: MAXWRK = MAX( MAXWRK,
368: $ N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
369: MAXWRK = MAX( MAXWRK,
370: $ N + N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
371: IF( ILVL ) THEN
372: MAXWRK = MAX( MAXWRK, N +
373: $ N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, 0 ) )
374: END IF
375: END IF
376: WORK( 1 ) = MAXWRK
377: *
378: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
379: INFO = -25
380: END IF
381: END IF
382: *
383: IF( INFO.NE.0 ) THEN
384: CALL XERBLA( 'ZGGEVX', -INFO )
385: RETURN
386: ELSE IF( LQUERY ) THEN
387: RETURN
388: END IF
389: *
390: * Quick return if possible
391: *
392: IF( N.EQ.0 )
393: $ RETURN
394: *
395: * Get machine constants
396: *
397: EPS = DLAMCH( 'P' )
398: SMLNUM = DLAMCH( 'S' )
399: BIGNUM = ONE / SMLNUM
400: CALL DLABAD( SMLNUM, BIGNUM )
401: SMLNUM = SQRT( SMLNUM ) / EPS
402: BIGNUM = ONE / SMLNUM
403: *
404: * Scale A if max element outside range [SMLNUM,BIGNUM]
405: *
406: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
407: ILASCL = .FALSE.
408: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
409: ANRMTO = SMLNUM
410: ILASCL = .TRUE.
411: ELSE IF( ANRM.GT.BIGNUM ) THEN
412: ANRMTO = BIGNUM
413: ILASCL = .TRUE.
414: END IF
415: IF( ILASCL )
416: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
417: *
418: * Scale B if max element outside range [SMLNUM,BIGNUM]
419: *
420: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
421: ILBSCL = .FALSE.
422: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
423: BNRMTO = SMLNUM
424: ILBSCL = .TRUE.
425: ELSE IF( BNRM.GT.BIGNUM ) THEN
426: BNRMTO = BIGNUM
427: ILBSCL = .TRUE.
428: END IF
429: IF( ILBSCL )
430: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
431: *
432: * Permute and/or balance the matrix pair (A,B)
433: * (Real Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
434: *
435: CALL ZGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
436: $ RWORK, IERR )
437: *
438: * Compute ABNRM and BBNRM
439: *
440: ABNRM = ZLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
441: IF( ILASCL ) THEN
442: RWORK( 1 ) = ABNRM
443: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
444: $ IERR )
445: ABNRM = RWORK( 1 )
446: END IF
447: *
448: BBNRM = ZLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
449: IF( ILBSCL ) THEN
450: RWORK( 1 ) = BBNRM
451: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, RWORK( 1 ), 1,
452: $ IERR )
453: BBNRM = RWORK( 1 )
454: END IF
455: *
456: * Reduce B to triangular form (QR decomposition of B)
457: * (Complex Workspace: need N, prefer N*NB )
458: *
459: IROWS = IHI + 1 - ILO
460: IF( ILV .OR. .NOT.WANTSN ) THEN
461: ICOLS = N + 1 - ILO
462: ELSE
463: ICOLS = IROWS
464: END IF
465: ITAU = 1
466: IWRK = ITAU + IROWS
467: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
468: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
469: *
470: * Apply the unitary transformation to A
471: * (Complex Workspace: need N, prefer N*NB)
472: *
473: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
474: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
475: $ LWORK+1-IWRK, IERR )
476: *
477: * Initialize VL and/or VR
478: * (Workspace: need N, prefer N*NB)
479: *
480: IF( ILVL ) THEN
481: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
482: IF( IROWS.GT.1 ) THEN
483: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
484: $ VL( ILO+1, ILO ), LDVL )
485: END IF
486: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
487: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
488: END IF
489: *
490: IF( ILVR )
491: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
492: *
493: * Reduce to generalized Hessenberg form
494: * (Workspace: none needed)
495: *
496: IF( ILV .OR. .NOT.WANTSN ) THEN
497: *
498: * Eigenvectors requested -- work on whole matrix.
499: *
500: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
501: $ LDVL, VR, LDVR, IERR )
502: ELSE
503: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
504: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
505: END IF
506: *
507: * Perform QZ algorithm (Compute eigenvalues, and optionally, the
508: * Schur forms and Schur vectors)
509: * (Complex Workspace: need N)
510: * (Real Workspace: need N)
511: *
512: IWRK = ITAU
513: IF( ILV .OR. .NOT.WANTSN ) THEN
514: CHTEMP = 'S'
515: ELSE
516: CHTEMP = 'E'
517: END IF
518: *
519: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
520: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
521: $ LWORK+1-IWRK, RWORK, IERR )
522: IF( IERR.NE.0 ) THEN
523: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
524: INFO = IERR
525: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
526: INFO = IERR - N
527: ELSE
528: INFO = N + 1
529: END IF
530: GO TO 90
531: END IF
532: *
533: * Compute Eigenvectors and estimate condition numbers if desired
534: * ZTGEVC: (Complex Workspace: need 2*N )
535: * (Real Workspace: need 2*N )
536: * ZTGSNA: (Complex Workspace: need 2*N*N if SENSE='V' or 'B')
537: * (Integer Workspace: need N+2 )
538: *
539: IF( ILV .OR. .NOT.WANTSN ) THEN
540: IF( ILV ) THEN
541: IF( ILVL ) THEN
542: IF( ILVR ) THEN
543: CHTEMP = 'B'
544: ELSE
545: CHTEMP = 'L'
546: END IF
547: ELSE
548: CHTEMP = 'R'
549: END IF
550: *
551: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
552: $ LDVL, VR, LDVR, N, IN, WORK( IWRK ), RWORK,
553: $ IERR )
554: IF( IERR.NE.0 ) THEN
555: INFO = N + 2
556: GO TO 90
557: END IF
558: END IF
559: *
560: IF( .NOT.WANTSN ) THEN
561: *
562: * compute eigenvectors (DTGEVC) and estimate condition
563: * numbers (DTGSNA). Note that the definition of the condition
564: * number is not invariant under transformation (u,v) to
565: * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
566: * Schur form (S,T), Q and Z are orthogonal matrices. In order
567: * to avoid using extra 2*N*N workspace, we have to
568: * re-calculate eigenvectors and estimate the condition numbers
569: * one at a time.
570: *
571: DO 20 I = 1, N
572: *
573: DO 10 J = 1, N
574: BWORK( J ) = .FALSE.
575: 10 CONTINUE
576: BWORK( I ) = .TRUE.
577: *
578: IWRK = N + 1
579: IWRK1 = IWRK + N
580: *
581: IF( WANTSE .OR. WANTSB ) THEN
582: CALL ZTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
583: $ WORK( 1 ), N, WORK( IWRK ), N, 1, M,
584: $ WORK( IWRK1 ), RWORK, IERR )
585: IF( IERR.NE.0 ) THEN
586: INFO = N + 2
587: GO TO 90
588: END IF
589: END IF
590: *
591: CALL ZTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
592: $ WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
593: $ RCONDV( I ), 1, M, WORK( IWRK1 ),
594: $ LWORK-IWRK1+1, IWORK, IERR )
595: *
596: 20 CONTINUE
597: END IF
598: END IF
599: *
600: * Undo balancing on VL and VR and normalization
601: * (Workspace: none needed)
602: *
603: IF( ILVL ) THEN
604: CALL ZGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
605: $ LDVL, IERR )
606: *
607: DO 50 JC = 1, N
608: TEMP = ZERO
609: DO 30 JR = 1, N
610: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
611: 30 CONTINUE
612: IF( TEMP.LT.SMLNUM )
613: $ GO TO 50
614: TEMP = ONE / TEMP
615: DO 40 JR = 1, N
616: VL( JR, JC ) = VL( JR, JC )*TEMP
617: 40 CONTINUE
618: 50 CONTINUE
619: END IF
620: *
621: IF( ILVR ) THEN
622: CALL ZGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
623: $ LDVR, IERR )
624: DO 80 JC = 1, N
625: TEMP = ZERO
626: DO 60 JR = 1, N
627: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
628: 60 CONTINUE
629: IF( TEMP.LT.SMLNUM )
630: $ GO TO 80
631: TEMP = ONE / TEMP
632: DO 70 JR = 1, N
633: VR( JR, JC ) = VR( JR, JC )*TEMP
634: 70 CONTINUE
635: 80 CONTINUE
636: END IF
637: *
638: * Undo scaling if necessary
639: *
640: IF( ILASCL )
641: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
642: *
643: IF( ILBSCL )
644: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
645: *
646: 90 CONTINUE
647: WORK( 1 ) = MAXWRK
648: *
649: RETURN
650: *
651: * End of ZGGEVX
652: *
653: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggevx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack