Program e02adfe

!     E02ADF Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

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