Program e02rafe

!     E02RAF Example Program Text

!     Mark 28.6 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: c02abf, e02raf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: k = 4, m = 4, nout = 6
      Integer, Parameter               :: ia = k + 1
      Integer, Parameter               :: ib = m + 1
      Integer, Parameter               :: ic = ia + ib - 1
      Integer, Parameter               :: jw = ib*(2*ib+3)
!     .. Local Scalars ..
      Integer                          :: i, ifail, itmax, polish
!     .. Local Arrays ..
      Complex (Kind=nag_wp)            :: z(max(k,m))
      Real (Kind=nag_wp)               :: a(ia), b(ib), berr(max(k,m)), c(ic), &
                                          cond(max(k,m)), dd(ia+ib), w(jw)
      Integer                          :: conv(max(k,m))
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, real
!     .. Executable Statements ..
      Write (nout,*) 'E02RAF Example Program Results'

!     Power series coefficients in C

      c(1) = 1.0E0_nag_wp

      Do i = 1, ic - 1
        c(i+1) = c(i)/real(i,kind=nag_wp)
      End Do

      ifail = 0
      Call e02raf(ia,ib,c,ic,a,b,w,jw,ifail)

      Write (nout,*)
      Write (nout,*) 'The given series coefficients are'
      Write (nout,99999) c(1:ic)
      Write (nout,*)
      Write (nout,*) 'Numerator coefficients'
      Write (nout,99999) a(1:ia)
      Write (nout,*)
      Write (nout,*) 'Denominator coefficients'
      Write (nout,99999) b(1:ib)

!     Calculate zeros of the approximant using C02ABF
!     First need to reverse order of coefficients

      dd(ia:1:-1) = a(1:ia)

      itmax = 30
      polish = 0
      ifail = 0
      Call c02abf(dd,k,itmax,polish,z,berr,cond,conv,ifail)

      Write (nout,*)
      Write (nout,*) 'Zeros of approximant are at'
      Write (nout,*)
      Write (nout,*) '    Real part    Imag part'
      Write (nout,99998)(z(i),i=1,k)

!     Calculate poles of the approximant using C02ABF
!     Reverse order of coefficients

      dd(ib:1:-1) = b(1:ib)

      ifail = 0
      Call c02abf(dd,m,itmax,polish,z,berr,cond,conv,ifail)

      Write (nout,*)
      Write (nout,*) 'Poles of approximant are at'
      Write (nout,*)
      Write (nout,*) '    Real part    Imag part'
      Write (nout,99998)(z(i),i=1,k)

99999 Format (1X,5E13.4)
99998 Format (1X,2E13.4)
    End Program e02rafe