PROGRAM e01aefe

!      E01AEF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       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