PROGRAM e02adfe

!      E02ADF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : e02adf, e02aef, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: fit, x1, xarg, xcapr, xm
       INTEGER                         :: i, ifail, iwght, j, k, kplus1, lda,  &
                                          m, r
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), ak(:), s(:), w(:), work1(:), &
                                          work2(:), x(:), y(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'E02ADF Example Program Results'

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

       READ (nin,*) m
       READ (nin,*) k, iwght
       kplus1 = k + 1
       lda = kplus1
       ALLOCATE (a(lda,kplus1),s(kplus1),w(m),work1(3*m),work2(2*kplus1),x(m), &
          y(m))

       DO r = 1, m

          IF (iwght/=1) THEN
             READ (nin,*) x(r), y(r), w(r)
          ELSE
             READ (nin,*) x(r), y(r)
             w(r) = 1.0E0_nag_wp
          END IF

       END DO

       ifail = 0
       CALL e02adf(m,kplus1,lda,x,y,w,work1,work2,a,s,ifail)

       DO i = 0, k
          WRITE (nout,*)
          WRITE (nout,99998) 'Degree', i, '   R.M.S. residual =', s(i+1)
          WRITE (nout,*)
          WRITE (nout,*) '  J  Chebyshev coeff A(J)'
          WRITE (nout,99997) (j,a(i+1,j),j=1,i+1)
       END DO

       ALLOCATE (ak(kplus1))

       ak(1:kplus1) = a(kplus1,1:kplus1)
       x1 = x(1)
       xm = x(m)

       WRITE (nout,*)
       WRITE (nout,99996) 'Polynomial approximation and residuals for degree', &
          k
       WRITE (nout,*)
       WRITE (nout,*) &
          '  R   Abscissa     Weight   Ordinate  Polynomial  Residual'

       DO r = 1, m
          xcapr = ((x(r)-x1)-(xm-x(r)))/(xm-x1)

          ifail = 0
          CALL e02aef(kplus1,ak,xcapr,fit,ifail)

          WRITE (nout,99999) r, x(r), w(r), y(r), fit, fit - y(r)

          IF (r<m) THEN
             xarg = 0.5E0_nag_wp*(x(r)+x(r+1))
             xcapr = ((xarg-x1)-(xm-xarg))/(xm-x1)

             ifail = 0
             CALL e02aef(kplus1,ak,xcapr,fit,ifail)

             WRITE (nout,99995) xarg, fit
          END IF

       END DO

99999  FORMAT (1X,I3,4F11.4,E11.2)
99998  FORMAT (1X,A,I4,A,E12.2)
99997  FORMAT (1X,I3,F15.4)
99996  FORMAT (1X,A,I4)
99995  FORMAT (4X,F11.4,22X,F11.4)
    END PROGRAM e02adfe