*     E02RAF Example Program Text
*     Mark 16 Revised. NAG Copyright 1993.
*     .. Parameters ..
      INTEGER          L, M, IA, IB, IC, JW
      PARAMETER        (L=4,M=4,IA=L+1,IB=M+1,IC=IA+IB-1,JW=IB*(2*IB+3))
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      LOGICAL          SCALE
      PARAMETER        (SCALE=.TRUE.)
*     .. Local Scalars ..
      INTEGER          I, IFAIL
*     .. Local Arrays ..
      DOUBLE PRECISION A(IA), B(IB), C(IC), DD(IA+IB), W(JW),
     +                 WORK(2*(L+M+1)), Z(2,L+M)
*     .. External Subroutines ..
      EXTERNAL         C02AGF, E02RAF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'E02RAF Example Program Results'
*     Power series coefficients in C
      C(1) = 1.0D0
      DO 20 I = 1, IC - 1
         C(I+1) = C(I)/DBLE(I)
   20 CONTINUE
      IFAIL = 1
*
      CALL E02RAF(IA,IB,C,IC,A,B,W,JW,IFAIL)
*
      WRITE (NOUT,*)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99997) ' ** E02RAF returned with IFAIL = ', IFAIL
         GO TO 80
      END IF
*
      WRITE (NOUT,*) 'The given series coefficients are'
      WRITE (NOUT,99999) (C(I),I=1,IC)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Numerator coefficients'
      WRITE (NOUT,99999) (A(I),I=1,IA)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Denominator coefficients'
      WRITE (NOUT,99999) (B(I),I=1,IB)
*     Calculate zeros of the approximant using C02AGF
*     First need to reverse order of coefficients
      DO 40 I = 1, IA
         DD(IA-I+1) = A(I)
   40 CONTINUE
      IFAIL = 0
*
      CALL C02AGF(DD,L,SCALE,Z,WORK,IFAIL)
*
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Zeros of approximant are at'
      WRITE (NOUT,*)
      WRITE (NOUT,*) '    Real part    Imag part'
      WRITE (NOUT,99998) (Z(1,I),Z(2,I),I=1,L)
*     Calculate poles of the approximant using C02AGF
*     Reverse order of coefficients
      DO 60 I = 1, IB
         DD(IB-I+1) = B(I)
   60 CONTINUE
      IFAIL = 0
*
      CALL C02AGF(DD,M,SCALE,Z,WORK,IFAIL)
*
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Poles of approximant are at'
      WRITE (NOUT,*)
      WRITE (NOUT,*) '    Real part    Imag part'
      WRITE (NOUT,99998) (Z(1,I),Z(2,I),I=1,M)
   80 CONTINUE
*
99999 FORMAT (1X,5E13.4)
99998 FORMAT (1X,2E13.4)
99997 FORMAT (1X,A,I5)
      END