!   D05ABF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module d05abfe_mod

!     D05ABF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: g, k
!     .. Parameters ..
      Integer, Parameter, Public       :: nmax = 10, nout = 6
    Contains
      Function k(x,s)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: k
!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: alpha = 1.0_nag_wp
        Real (Kind=nag_wp), Parameter  :: w = alpha**2
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: s, x
!       .. Executable Statements ..
        k = alpha/(w+(x-s)*(x-s))

        Return
      End Function k
      Function g(x)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: g
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Executable Statements ..
        g = 1.0_nag_wp

        Return
      End Function g
    End Module d05abfe_mod
    Program d05abfe

!     D05ABF Example Main Program

!     .. Use Statements ..
      Use d05abfe_mod, Only: g, k, nmax, nout
      Use nag_library, Only: c06dcf, d05abf, nag_wp, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: a, b, lambda, x0
      Integer                          :: i, ifail, ldcm, lx, n, nt2p1, ss
      Logical                          :: ev, odorev
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:), chebr(:), cm(:,:), f(:),       &
                                          f1(:,:), wk(:,:), x(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, int, real
!     .. Executable Statements ..
      Write (nout,*) 'D05ABF Example Program Results'

      odorev = .True.
      ev = .True.
      lambda = -0.3183_nag_wp
      a = -1.0_nag_wp
      b = 1.0_nag_wp

      If (odorev) Then
        Write (nout,*)
        If (ev) Then
          Write (nout,*) 'Solution is even'
          ss = 2
        Else
          Write (nout,*) 'Solution is odd'
          ss = 3
        End If
        x0 = 0.5_nag_wp*(a+b)
      Else
        ss = 1
        x0 = a
      End If

!     Set up uniform grid to evaluate Chebyshev polynomials.
      lx = int(4.000001_nag_wp*(b-x0)) + 1
      Allocate (x(lx),chebr(lx))
      x(1) = x0
      Do i = 2, lx
        x(i) = x(i-1) + 0.25_nag_wp
      End Do

      Do n = 5, nmax, 5
        ldcm = n
        nt2p1 = 2*n + 1
        Allocate (c(n),cm(ldcm,ldcm),f(n),f1(ldcm,1),wk(2,nt2p1))

        ifail = -1
        Call d05abf(k,g,lambda,a,b,odorev,ev,n,cm,f1,wk,ldcm,nt2p1,f,c,ifail)

        If (ifail==0) Then
          Write (nout,*)
          Write (nout,99999) 'Results for N =', n
          Write (nout,*)
          Write (nout,99996) 'Solution on first ', n,                          &
            ' Chebyshev points and Chebyshev coefficients'
          Write (nout,*) '  I           X             F(I)         C(I)'
          Write (nout,99998)(i,cos(x01aaf(a)*real(i,kind=nag_wp)/real(2*n-1,   &
            kind=nag_wp)),f(i),c(i),i=1,n)

!         Evaluate and print solution on uniform grid.
          ifail = 0
          Call c06dcf(x,lx,a,b,c,n,ss,chebr,ifail)

          Write (nout,*)
          Write (nout,*) 'Solution on evenly spaced grid'
          Write (nout,*) '     X           F(X)'
          Write (nout,99997)(x(i),chebr(i),i=1,lx)

        End If

        Deallocate (c,cm,f,f1,wk)
      End Do
      Deallocate (x,chebr)

99999 Format (1X,A,I3)
99998 Format (1X,I3,2F15.5,E15.5)
99997 Format (1X,F8.4,F15.5)
99996 Format (1X,A,I2,A)
    End Program d05abfe