Program e01aefe

!     E01AEF Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: e01aef, f16dnf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: xmax, xmin
      Integer                          :: i, ifail, ip1, ipmax, ires, itmax,   &
                                          itmin, iy, j, k, liwrk, lwrk, m, n
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), wrk(:), x(:), y(:)
      Integer, Allocatable             :: ip(:), iwrk(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: sum
!     .. Executable Statements ..
      Write (nout,*) 'E01AEF Example Program Results'

!     Skip heading in data file
      Read (nin,*)

      itmin = -1
      itmax = -1

      Read (nin,*) m, xmin, xmax
      liwrk = 2*(m+1)
      Allocate (ip(m),x(m),iwrk(liwrk))

      Read (nin,*)(ip(i),i=1,m)
      Read (nin,*)(x(i),i=1,m)

      n = m + sum(ip(1:m))

!     Get the maximum value of IP

      Call f16dnf(m,ip,1,k,ipmax)

      lwrk = 7*n + 5*ipmax + m + 7
      Allocate (a(n),y(n),wrk(lwrk))

      j = 0

      Do i = 1, m
        Read (nin,*)(y(k),k=j+1,j+ip(i)+1)
        j = j + ip(i) + 1
      End Do

      ifail = -1
      Call e01aef(m,xmin,xmax,x,y,ip,n,itmin,itmax,a,wrk,lwrk,iwrk,liwrk, &
        ifail)

      Write (nout,*)

      Select Case (ifail)
      Case (0,4:)
        Write (nout,99999) 'Total number of interpolating conditions =', n
        Write (nout,*)
        Write (nout,*) 'Interpolating polynomial'
        Write (nout,*)
        Write (nout,*) '   I    Chebyshev Coefficient A(I+1)'

        Do i = 1, n
          Write (nout,99998) i - 1, a(i)
        End Do

        Write (nout,*)
        Write (nout,*) '  X    R   Rth derivative    Residual'

        iy = 0
        ires = ipmax + 1

        Do i = 1, m
          ip1 = ip(i) + 1

          Do j = 1, ip1
            iy = iy + 1
            ires = ires + 1

            If (j-1/=0) Then
              Write (nout,99997) j - 1, y(iy), wrk(ires)
            Else
              Write (nout,99996) x(i), '   0', y(iy), wrk(ires)
            End If

          End Do

        End Do

      End Select

99999 Format (1X,A,I4)
99998 Format (1X,I4,F20.3)
99997 Format (5X,I4,F12.1,F17.6)
99996 Format (1X,F4.1,A,F12.1,F17.6)
    End Program e01aefe