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