*     F08XEF 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=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 A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
     +                 B(LDB,NMAX), BETA(NMAX), LSCALE(NMAX),
     +                 Q(LDQ,LDQ), RSCALE(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 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 DGGBAL(JOB,N,A,LDA,B,LDB,ILO,IHI,LSCALE,RSCALE,WORK,INFO)
*
*        Matrix A after balancing
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,A,LDA,'Matrix A after balancing',
     +               IFAIL)
         WRITE (NOUT,*)
*
*        Matrix B after balancing
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,B,LDB,'Matrix B after balancing',
     +               IFAIL)
         WRITE (NOUT,*)
*
*         Reduce B to triangular form using QR
*
         IROWS = IHI + 1 - ILO
         CALL DGEQRF(IROWS,IROWS,B(ILO,ILO),LDB,TAU,WORK,LWORK,INFO)
*
*         Apply the orthogonal transformation to matrix A
*
         CALL DORMQR('L','T',IROWS,IROWS,IROWS,B(ILO,ILO),LDB,TAU,
     +               A(ILO,ILO),LDA,WORK,LWORK,INFO)
*
*        Compute the generalized Hessenberg form of (A,B)
*
         COMPQ = 'N'
         COMPZ = 'N'
         CALL DGGHRD(COMPQ,COMPZ,IROWS,1,IROWS,A(ILO,ILO),LDA,B(ILO,ILO)
     +               ,LDB,Q,LDQ,Z,LDZ,INFO)
*
*        Matrix A in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,A,LDA,
     +               'Matrix A in Hessenberg form',IFAIL)
         WRITE (NOUT,*)
*
*        Matrix B in generalized Hessenberg form
*
         IFAIL = 0
         CALL X04CAF('General',' ',N,N,B,LDB,'Matrix B is triangular',
     +               IFAIL)
*
*        Routine DHGEQZ
*        Workspace query: JWORK = -1
*
         JWORK = -1
         JOB = 'E'
         CALL DHGEQZ(JOB,COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,LDB,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,A,LDA,B,LDB,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
      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,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
