*     F08YXF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, LDA, LDB, LDQ, LDZ, LWORK
      PARAMETER        (NMAX=10,LDA=NMAX,LDB=NMAX,LDQ=NMAX,LDZ=NMAX,
     +                 LWORK=6*NMAX)
      COMPLEX *16      CONE, CZERO
      PARAMETER        (CONE=1.0D0,CZERO=0.0D0)
*     .. Local Scalars ..
      COMPLEX *16      E
      INTEGER          I, ICOLS, IFAIL, IHI, ILO, INFO, IROWS, J, JWORK,
     +                 M, N
      LOGICAL          ILEFT, IRIGHT
      CHARACTER        COMPQ, COMPZ, HOWMNY, JOB, SIDE
*     .. Local Arrays ..
      COMPLEX *16      A(LDA,NMAX), ALPHA(NMAX), B(LDB,NMAX),
     +                 BETA(NMAX), Q(LDQ,LDQ), TAU(NMAX), WORK(LWORK),
     +                 Z(LDZ,LDZ)
      DOUBLE PRECISION LSCALE(NMAX), RSCALE(NMAX), RWORK(6*NMAX)
      LOGICAL          SELECT(NMAX)
      CHARACTER        CLABS(1), RLABS(1)
*     .. External Subroutines ..
      EXTERNAL         F06TFF, F06THF, X04DBF, ZGEQRF, ZGGBAK, ZGGBAL,
     +                 ZGGHRD, ZHGEQZ, ZTGEVC, ZUNGQR, ZUNMQR
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE, DIMAG, NINT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08YXF Example Program Results'
*
*     ILEFT  is TRUE if left  eigenvectors are required
*     IRIGHT is TRUE if right eigenvectors are required
*
      ILEFT = .TRUE.
      IRIGHT = .TRUE.
*
*     Skip heading in data file
*
      READ (NIN,*)
      READ (NIN,*) N
      IF (N.LE.NMAX) THEN
*
*        READ matrix A from data file
*
         READ (NIN,*) ((A(I,J),J=1,N),I=1,N)
*
*        READ matrix B from data file
*
         READ (NIN,*) ((B(I,J),J=1,N),I=1,N)
*
*        Balance matrix pair (A,B)
*
         JOB = 'B'
         CALL ZGGBAL(JOB,N,A,LDA,B,LDB,ILO,IHI,LSCALE,RSCALE,RWORK,INFO)
*
*        Matrix A after balancing
*
         IFAIL = 0
         CALL X04DBF('General',' ',N,N,A,LDA,'Bracketed','F7.4',
     +               'Matrix A after balancing','Integer',RLABS,
     +               'Integer',CLABS,80,0,IFAIL)
         WRITE (NOUT,*)
*
*        Matrix B after balancing
*
         IFAIL = 0
         CALL X04DBF('General',' ',N,N,B,LDB,'Bracketed','F7.4',
     +               'Matrix B after balancing','Integer',RLABS,
     +               'Integer',CLABS,80,0,IFAIL)
         WRITE (NOUT,*)
*
*        Reduce B to triangular form using QR
*
         IROWS = IHI + 1 - ILO
         ICOLS = N + 1 - ILO
         CALL ZGEQRF(IROWS,ICOLS,B(ILO,ILO),LDB,TAU,WORK,LWORK,INFO)
*
*        Apply the orthogonal transformation to A
*
         CALL ZUNMQR('L','C',IROWS,ICOLS,IROWS,B(ILO,ILO),LDB,TAU,
     +               A(ILO,ILO),LDA,WORK,LWORK,INFO)
*
*        Initialize Q (for left eigenvectors)
*
         IF (ILEFT) THEN
*
            CALL F06THF('General',N,N,CZERO,CONE,Q,LDQ)
            CALL F06TFF('Lower',IROWS-1,IROWS-1,B(ILO+1,ILO),LDB,
     +                  Q(ILO+1,ILO),LDQ)
            CALL ZUNGQR(IROWS,IROWS,IROWS,Q(ILO,ILO),LDQ,TAU,WORK,LWORK,
     +                  INFO)
         END IF
*
*        Initialize Z for right eigenvectors
*
         IF (IRIGHT) CALL F06THF('General',N,N,CZERO,CONE,Z,LDZ)
*
*        Compute the generalized Hessenberg form of (A,B)
*
         COMPQ = 'V'
         COMPZ = 'V'
         CALL ZGGHRD(COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,Q,LDQ,Z,LDZ,INFO)
*
*        Matrix A in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04DBF('General',' ',N,N,A,LDA,'Bracketed','F7.3',
     +               'Matrix A in Hessenberg form','Integer',RLABS,
     +               'Integer',CLABS,80,0,IFAIL)
         WRITE (NOUT,*)
*
*        Matrix B in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04DBF('General',' ',N,N,B,LDB,'Bracketed','F7.3',
     +               'Matrix B in Hessenberg form','Integer',RLABS,
     +               'Integer',CLABS,80,0,IFAIL)
*
*        Routine ZHGEQZ
*        Workspace query: JWORK = -1
*
         JWORK = -1
         JOB = 'S'
         CALL ZHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,ALPHA,BETA,Q,
     +               LDQ,Z,LDZ,WORK,JWORK,RWORK,INFO)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) NINT(DBLE(WORK(1)))
         WRITE (NOUT,99998) LWORK
         WRITE (NOUT,*)
*
*        Compute the generalized Schur form
*        if the workspace LWORK is adequate
*
         IF (NINT(DBLE(WORK(1))).LE.LWORK) THEN
            CALL ZHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,ALPHA,
     +                  BETA,Q,LDQ,Z,LDZ,WORK,LWORK,RWORK,INFO)
*
*        Print the generalized eigenvalues
*        Note: the actual values of beta are real and non-negative
*
            WRITE (NOUT,99997)
            DO 20 I = 1, N
               IF (DBLE(BETA(I)).NE.0.0D0) THEN
                  E = ALPHA(I)/BETA(I)
                  WRITE (NOUT,99995) I, '(', DBLE(E), ',', DIMAG(E), ')'
               ELSE
                  WRITE (NOUT,99996) I
               END IF
   20       CONTINUE
            WRITE (NOUT,*)
*
*        Compute left and right generalized eigenvectors
*        of the balanced matrix
*
            HOWMNY = 'B'
            IF (ILEFT .AND. IRIGHT) THEN
               SIDE = 'B'
            ELSE IF (ILEFT) THEN
               SIDE = 'L'
            ELSE IF (IRIGHT) THEN
               SIDE = 'R'
            END IF
*
            CALL ZTGEVC(SIDE,HOWMNY,SELECT,N,A,LDA,B,LDB,Q,LDQ,Z,LDZ,N,
     +                  M,WORK,RWORK,INFO)
*
*        Compute right eigenvectors of the original matrix
*
            IF (IRIGHT) THEN
               JOB = 'B'
               SIDE = 'R'
*
               CALL ZGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,Z,LDZ,
     +                     INFO)
*
*        Print the right eigenvectors
*
               IFAIL = 0
               CALL X04DBF('General',' ',N,N,Z,LDZ,'Bracketed','F7.4',
     +                     'Right eigenvectors','Integer',RLABS,
     +                     'Integer',CLABS,80,0,IFAIL)
               WRITE (NOUT,*)
            END IF
*
*        Compute left eigenvectors of the original matrix
*
            IF (IRIGHT) THEN
               JOB = 'B'
               SIDE = 'L'
*
               CALL ZGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,Q,LDQ,
     +                     INFO)
*
*        Print the left eigenvectors
*
               IFAIL = 0
               CALL X04DBF('General',' ',N,N,Q,LDQ,'Bracketed','F7.4',
     +                     'Left eigenvectors','Integer',RLABS,
     +                     'Integer',CLABS,80,0,IFAIL)
            END IF
         ELSE
            WRITE (NOUT,99994)
         END IF
      END IF
      STOP
*
99999 FORMAT (1X,'Minimal required LWORK = ',I6)
99998 FORMAT (1X,'Actual value of  LWORK = ',I6)
99997 FORMAT (1X,'Generalized eigenvalues')
99996 FORMAT (1X,I4,' Infinite eigenvalue')
99995 FORMAT (1X,I4,5X,A,F7.3,A,F7.3,A)
99994 FORMAT (1X,'Insufficient workspace for array WORK',/' in F08XSF/',
     +       'ZHGEQZ')
      END
