Program e01aefe

!     E01AEF Example Program Text

!     Mark 28.5 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: e01aef, f16dnf, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: restol, xmax, xmin
      Integer                          :: i, ifail, 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                        :: abs, 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                a_i'

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

        Write (nout,*)
        Write (nout,*) '  i   x_i    order(k)   y^(k)    Residual'

!       Residuals less than 100 times machine precision are not printed.
        restol = x02ajf()*100.0_nag_wp
        iy = 0
        ires = ipmax + 1

        Do i = 1, m

          iy = iy + 1
          ires = ires + 1
          If (abs(wrk(ires))>restol) Then
            Write (nout,99996) i, x(i), 0, y(iy), wrk(ires)
          Else
            Write (nout,99995) i, x(i), 0, y(iy), ' - '
          End If
          Do j = 1, ip(i)
            iy = iy + 1
            ires = ires + 1
            If (abs(wrk(ires))>restol) Then
              Write (nout,99997) j, y(iy), wrk(ires)
            Else
              Write (nout,99994) j, y(iy), ' - '
            End If
          End Do

        End Do

      End Select

99999 Format (1X,A,I4)
99998 Format (1X,I4,F20.4)
99997 Format (10X,I7,F11.1,E19.10)
99996 Format (1X,I3,F6.1,I7,F11.1,E19.10)
99995 Format (1X,I3,F6.1,I7,F11.1,A11)
99994 Format (10X,I7,F11.1,A11)
    End Program e01aefe