001:       SUBROUTINE CGGEVX( 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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
008: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
009: *     November 2006
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
013:       INTEGER            IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
014:       REAL               ABNRM, BBNRM
015: *     ..
016: *     .. Array Arguments ..
017:       LOGICAL            BWORK( * )
018:       INTEGER            IWORK( * )
019:       REAL               LSCALE( * ), RCONDE( * ), RCONDV( * ),
020:      $                   RSCALE( * ), RWORK( * )
021:       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
022:      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
023:      $                   WORK( * )
024: *     ..
025: *
026: *  Purpose
027: *  =======
028: *
029: *  CGGEVX computes for a pair of N-by-N complex nonsymmetric matrices
030: *  (A,B) the generalized eigenvalues, and optionally, the left and/or
031: *  right generalized eigenvectors.
032: *
033: *  Optionally, it also computes a balancing transformation to improve
034: *  the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
035: *  LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
036: *  the eigenvalues (RCONDE), and reciprocal condition numbers for the
037: *  right eigenvectors (RCONDV).
038: *
039: *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar
040: *  lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
041: *  singular. It is usually represented as the pair (alpha,beta), as
042: *  there is a reasonable interpretation for beta=0, and even for both
043: *  being zero.
044: *
045: *  The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
046: *  of (A,B) satisfies
047: *                   A * v(j) = lambda(j) * B * v(j) .
048: *  The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
049: *  of (A,B) satisfies
050: *                   u(j)**H * A  = lambda(j) * u(j)**H * B.
051: *  where u(j)**H is the conjugate-transpose of u(j).
052: *
053: *
054: *  Arguments
055: *  =========
056: *
057: *  BALANC  (input) CHARACTER*1
058: *          Specifies the balance option to be performed:
059: *          = 'N':  do not diagonally scale or permute;
060: *          = 'P':  permute only;
061: *          = 'S':  scale only;
062: *          = 'B':  both permute and scale.
063: *          Computed reciprocal condition numbers will be for the
064: *          matrices after permuting and/or balancing. Permuting does
065: *          not change condition numbers (in exact arithmetic), but
066: *          balancing does.
067: *
068: *  JOBVL   (input) CHARACTER*1
069: *          = 'N':  do not compute the left generalized eigenvectors;
070: *          = 'V':  compute the left generalized eigenvectors.
071: *
072: *  JOBVR   (input) CHARACTER*1
073: *          = 'N':  do not compute the right generalized eigenvectors;
074: *          = 'V':  compute the right generalized eigenvectors.
075: *
076: *  SENSE   (input) CHARACTER*1
077: *          Determines which reciprocal condition numbers are computed.
078: *          = 'N': none are computed;
079: *          = 'E': computed for eigenvalues only;
080: *          = 'V': computed for eigenvectors only;
081: *          = 'B': computed for eigenvalues and eigenvectors.
082: *
083: *  N       (input) INTEGER
084: *          The order of the matrices A, B, VL, and VR.  N >= 0.
085: *
086: *  A       (input/output) COMPLEX array, dimension (LDA, N)
087: *          On entry, the matrix A in the pair (A,B).
088: *          On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
089: *          or both, then A contains the first part of the complex Schur
090: *          form of the "balanced" versions of the input A and B.
091: *
092: *  LDA     (input) INTEGER
093: *          The leading dimension of A.  LDA >= max(1,N).
094: *
095: *  B       (input/output) COMPLEX array, dimension (LDB, N)
096: *          On entry, the matrix B in the pair (A,B).
097: *          On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
098: *          or both, then B contains the second part of the complex
099: *          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 array, dimension (N)
105: *  BETA    (output) COMPLEX 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 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 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) REAL 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) REAL 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) REAL
170: *          The one-norm of the balanced matrix A.
171: *
172: *  BBNRM   (output) REAL
173: *          The one-norm of the balanced matrix B.
174: *
175: *  RCONDE  (output) REAL 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) REAL array, dimension (N)
181: *          If SENSE = '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 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 CHGEQZ.
220: *                =N+2: error return from CTGEVC.
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:       REAL               ZERO, ONE
250:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
251:       COMPLEX            CZERO, CONE
252:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
253:      $                   CONE = ( 1.0E+0, 0.0E+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:       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
262:      $                   SMLNUM, TEMP
263:       COMPLEX            X
264: *     ..
265: *     .. Local Arrays ..
266:       LOGICAL            LDUMMA( 1 )
267: *     ..
268: *     .. External Subroutines ..
269:       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
270:      $                   CLASCL, CLASET, CTGEVC, CTGSNA, CUNGQR, CUNMQR,
271:      $                   SLABAD, SLASCL, XERBLA
272: *     ..
273: *     .. External Functions ..
274:       LOGICAL            LSAME
275:       INTEGER            ILAENV
276:       REAL               CLANGE, SLAMCH
277:       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
278: *     ..
279: *     .. Intrinsic Functions ..
280:       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
281: *     ..
282: *     .. Statement Functions ..
283:       REAL               ABS1
284: *     ..
285: *     .. Statement Function definitions ..
286:       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( 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, 'CGEQRF', ' ', N, 1, N, 0 ) )
369:             MAXWRK = MAX( MAXWRK,
370:      $                    N + N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, 0 ) )
371:             IF( ILVL ) THEN
372:                MAXWRK = MAX( MAXWRK, N +
373:      $                       N*ILAENV( 1, 'CUNGQR', ' ', 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( 'CGGEVX', -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 = SLAMCH( 'P' )
398:       SMLNUM = SLAMCH( 'S' )
399:       BIGNUM = ONE / SMLNUM
400:       CALL SLABAD( 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 = CLANGE( '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 CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
417: *
418: *     Scale B if max element outside range [SMLNUM,BIGNUM]
419: *
420:       BNRM = CLANGE( '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 CLASCL( '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 CGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
436:      $             RWORK, IERR )
437: *
438: *     Compute ABNRM and BBNRM
439: *
440:       ABNRM = CLANGE( '1', N, N, A, LDA, RWORK( 1 ) )
441:       IF( ILASCL ) THEN
442:          RWORK( 1 ) = ABNRM
443:          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, RWORK( 1 ), 1,
444:      $                IERR )
445:          ABNRM = RWORK( 1 )
446:       END IF
447: *
448:       BBNRM = CLANGE( '1', N, N, B, LDB, RWORK( 1 ) )
449:       IF( ILBSCL ) THEN
450:          RWORK( 1 ) = BBNRM
451:          CALL SLASCL( '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 CGEQRF( 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 CUNMQR( '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 CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
482:          IF( IROWS.GT.1 ) THEN
483:             CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
484:      $                   VL( ILO+1, ILO ), LDVL )
485:          END IF
486:          CALL CUNGQR( 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 CLASET( '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 CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
501:      $                LDVL, VR, LDVR, IERR )
502:       ELSE
503:          CALL CGGHRD( '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 CHGEQZ( 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: *     CTGEVC: (Complex Workspace: need 2*N )
535: *             (Real Workspace:    need 2*N )
536: *     CTGSNA: (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 CTGEVC( 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 (STGEVC) and estimate condition
563: *           numbers (STGSNA). 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 CTGEVC( '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 CTGSNA( 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 CGGBAK( 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 CGGBAK( 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 CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
642: *
643:       IF( ILBSCL )
644:      $   CALL CLASCL( '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 CGGEVX
652: *
653:       END
654: