*     F08YKF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, LDA, LDB, LWORK
      PARAMETER        (NMAX=10,LDA=NMAX,LDB=NMAX,LWORK=6*NMAX)
      INTEGER          LDQ, LDZ
      PARAMETER        (LDQ=NMAX,LDZ=NMAX)
      DOUBLE PRECISION ONE, ZERO
      PARAMETER        (ONE=1.0D0,ZERO=0.0D0)
*     .. Local Scalars ..
      INTEGER          I, ICOLS, IFAIL, IHI, ILO, INFO, IROWS, J, JWORK,
     +                 M, N
      LOGICAL          ILEFT, IRIGHT
      CHARACTER        COMPQ, COMPZ, HOWMNY, JOB, SIDE
*     .. 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)
      LOGICAL          SELECT(NMAX)
*     .. External Subroutines ..
      EXTERNAL         DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DORGQR,
     +                 DORMQR, DTGEVC, F06QFF, F06QHF, X04CAF
*     .. Intrinsic Functions ..
      INTRINSIC        NINT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F08YKF 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 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
         ICOLS = N + 1 - ILO
         CALL DGEQRF(IROWS,ICOLS,B(ILO,ILO),LDB,TAU,WORK,LWORK,INFO)
*
*        Apply the orthogonal transformation to matrix A
*
         CALL DORMQR('L','T',IROWS,ICOLS,IROWS,B(ILO,ILO),LDB,TAU,
     +               A(ILO,ILO),LDA,WORK,LWORK,INFO)
*
*        Initialize Q (if left eigenvectors are required)
*
         IF (ILEFT) THEN
            CALL F06QHF('General',N,N,ZERO,ONE,Q,LDQ)
            CALL F06QFF('Lower',IROWS-1,IROWS-1,B(ILO+1,ILO),LDB,
     +                  Q(ILO+1,ILO),LDQ)
*
            CALL DORGQR(IROWS,IROWS,IROWS,Q(ILO,ILO),LDQ,TAU,WORK,LWORK,
     +                  INFO)
         END IF
*
*        Initialize Z (if right eigenvectors are required)
*
         IF (IRIGHT) CALL F06QHF('General',N,N,ZERO,ONE,Z,LDZ)
*
*        Compute the generalized Hessenberg form of (A,B)
*
         COMPQ = 'V'
         COMPZ = 'V'
         CALL DGGHRD(COMPQ,COMPZ,N,ILO,IHI,A,LDA,B,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 in Hessenberg form',IFAIL)
*
*        Routine DHGEQZ
*        Workspace query: JWORK = -1
*
         JWORK = -1
         JOB = 'S'
         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,*)
         WRITE (NOUT,99997)
         WRITE (NOUT,99996)
*
*        Compute the generalized Schur form
*        if the workspace LWORK is adequate
*        The Schur form also gives parameters
*        required to compute generalized eigenvalues
*
         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
*
            DO 20 I = 1, N
               IF (BETA(I).NE.0.0D0) THEN
                  WRITE (NOUT,99995) I, '(', ALPHAR(I)/BETA(I), ',',
     +              ALPHAI(I)/BETA(I), ')'
               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 DTGEVC(SIDE,HOWMNY,SELECT,N,A,LDA,B,LDB,Q,LDQ,Z,LDZ,N,
     +                  M,WORK,INFO)
            IF (IRIGHT) THEN
*
*        Compute right eigenvectors of the original matrix
*
               JOB = 'B'
               SIDE = 'R'
*
               CALL DGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,Z,LDZ,
     +                     INFO)
*
*        Print the right eigenvectors
*
               IFAIL = 0
               CALL X04CAF('General',' ',N,N,Z,LDZ,'Right eigenvectors',
     +                     IFAIL)
               WRITE (NOUT,*)
            END IF
*
*        Compute left eigenvectors of the original matrix
*
            IF (ILEFT) THEN
               JOB = 'B'
               SIDE = 'L'
*
               CALL DGGBAK(JOB,SIDE,N,ILO,IHI,LSCALE,RSCALE,N,Q,LDQ,
     +                     INFO)
*
*        Print the left eigenvectors
*
               IFAIL = 0
               CALL X04CAF('General',' ',N,N,Q,LDQ,'Left eigenvectors',
     +                     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,5X,'Infinite eigenvalue')
99995 FORMAT (1X,I4,5X,A,F7.3,A,F7.3,A)
99994 FORMAT (1X,'Insufficient workspace for array WORK',/' in F08XEF/',
     +       'DHGEQZ')
      END