*     F08XEF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, LDH, LDT, LDQ, LDZ, LWORK
      PARAMETER        (NMAX=10,LDH=NMAX,LDT=NMAX,LDQ=1,LDZ=1,
     +                 LWORK=6*NMAX)
*     .. Local Scalars ..
      INTEGER          I, IFAIL, IHI, ILO, INFO, IROWS, J, JWORK, N
      CHARACTER        COMPQ, COMPZ, JOB
*     .. Local Arrays ..
      DOUBLE PRECISION ALPHAI(NMAX), ALPHAR(NMAX), BETA(NMAX),
     +                 H(LDH,NMAX), LSCALE(NMAX), Q(LDQ,LDQ),
     +                 RSCALE(NMAX), T(LDT,NMAX), TAU(NMAX),
     +                 WORK(LWORK), Z(LDZ,LDZ)
*     .. External Subroutines ..
      EXTERNAL         DGEQRF, DGGBAL, DGGHRD, DHGEQZ, DORMQR, X04CAF
*     .. Intrinsic Functions ..
      INTRINSIC        NINT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08XEF Example Program Results'
*
*     Skip heading in data file
*
      READ (NIN,*)
      READ (NIN,*) N
      IF (N.LE.NMAX) THEN
*
*        READ matrix H from data file
*
         READ (NIN,*) ((H(I,J),J=1,N),I=1,N)
*
*        READ matrix T from data file
*
         READ (NIN,*) ((T(I,J),J=1,N),I=1,N)
*
*        Balance matrix pair (H,T)
*
         JOB = 'B'
         CALL DGGBAL(JOB,N,H,LDH,T,LDT,ILO,IHI,LSCALE,RSCALE,WORK,INFO)
*
*        Matrix H after balancing
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,H,LDH,'Matrix H after balancing',
     +               IFAIL)
         WRITE (NOUT,*)
*
*        Matrix T after balancing
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,T,LDT,'Matrix T after balancing',
     +               IFAIL)
         WRITE (NOUT,*)
*
*        Reduce T to triangular form using QR
*
         IROWS = IHI + 1 - ILO
         CALL DGEQRF(IROWS,IROWS,T(ILO,ILO),LDT,TAU,WORK,LWORK,INFO)
*
*        Apply the orthogonal transformation to matrix H
*
         CALL DORMQR('L','T',IROWS,IROWS,IROWS,T(ILO,ILO),LDT,TAU,
     +               H(ILO,ILO),LDH,WORK,LWORK,INFO)
*
*        Compute the generalized Hessenberg form of (H,T)
*
         COMPQ = 'N'
         COMPZ = 'N'
         CALL DGGHRD(COMPQ,COMPZ,IROWS,1,IROWS,H(ILO,ILO),LDH,T(ILO,ILO)
     +               ,LDT,Q,LDQ,Z,LDZ,INFO)
*
*        Matrix H in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,H,LDH,
     +               'Matrix H in Hessenberg form',IFAIL)
         WRITE (NOUT,*)
*
*        Matrix T in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,T,LDT,'Matrix T is triangular',
     +               IFAIL)
*
*        Routine DHGEQZ
*        Workspace query: JWORK = -1
*
         JWORK = -1
         JOB = 'E'
         CALL DHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,H,LDH,T,LDT,ALPHAR,
     +               ALPHAI,BETA,Q,LDQ,Z,LDZ,WORK,JWORK,INFO)
         WRITE (NOUT,*)
         WRITE (NOUT,99999) NINT(WORK(1))
         WRITE (NOUT,99998) LWORK
         WRITE (NOUT,*)
*
*        Compute the generalized Schur form
*        if the workspace LWORK is adequate
*
         IF (NINT(WORK(1)).LE.LWORK) THEN
            CALL DHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,H,LDH,T,LDT,ALPHAR,
     +                  ALPHAI,BETA,Q,LDQ,Z,LDZ,WORK,LWORK,INFO)
*
*           Print the generalized eigenvalues
*
            WRITE (NOUT,99997)
*
            DO 20 I = 1, N
               IF (BETA(I).NE.0.0D0) THEN
                  WRITE (NOUT,99996) I, '(', ALPHAR(I)/BETA(I), ',',
     +              ALPHAI(I)/BETA(I), ')'
               ELSE
                  WRITE (NOUT,99994) I
               END IF
   20       CONTINUE
         ELSE
            WRITE (NOUT,99995)
         END IF
      END IF
*
99999 FORMAT (1X,'Minimal required LWORK = ',I6)
99998 FORMAT (1X,'Actual value of  LWORK = ',I6)
99997 FORMAT (1X,'Generalized eigenvalues')
99996 FORMAT (1X,I4,5X,A,F7.3,A,F7.3,A)
99995 FORMAT (1X,'Insufficient workspace for array WORK',/' in F08XEF/',
     +       'DHGEQZ')
99994 FORMAT (1X,I4,'Eigenvalue is infinite')
      END