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