!   D02TGF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.

    Module d02tgfe_mod

!     D02TGF 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                           :: bdyc, coeff
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: coeff_tol = 1.0E-5_nag_wp
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Integer, Parameter, Public       :: n = 2, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: x0, x1
      Integer, Public, Save            :: k1
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable, Public, Save :: b(:,:)
      Integer, Public, Save            :: l(n) = (/1,2/)
      Integer, Public, Save            :: m(n) = (/1,2/)
    Contains
      Subroutine coeff(x,i,a,ia,ia1,rhs)

!       .. Use Statements ..
        Use nag_library, Only: e02akf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: rhs
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: i, ia, ia1
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: z1, z2
        Integer                        :: ifail
!       .. Executable Statements ..
        If (i<=1) Then

!         Evaluate z1, z2 at x using previous coeffs b.
          ifail = 0
          Call e02akf(k1,x0,x1,b(1,1),1,k1,x,z1,ifail)
          Call e02akf(k1,x0,x1,b(1,2),1,k1,x,z2,ifail)

          a(1,1) = z2*z2 - one
          a(1,2) = two
          a(2,1) = two*z1*z2 + one
          rhs = two*z1*z2*z2
        Else
          a(1,2) = -one
          a(2,3) = two
        End If
        Return
      End Subroutine coeff
      Subroutine bdyc(x,i,j,a,ia,ia1,rhs)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: rhs
        Real (Kind=nag_wp), Intent (Out) :: x
        Integer, Intent (In)           :: i, ia, ia1, j
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
!       .. Executable Statements ..
        x = -one
        a(i,j) = one
        If (i==2 .And. j==1) Then
          rhs = 3.0_nag_wp
        End If
        Return
      End Subroutine bdyc
    End Module d02tgfe_mod
    Program d02tgfe

!     D02TGF Example Main Program

!     .. Use Statements ..
      Use d02tgfe_mod, Only: b, bdyc, coeff, coeff_tol, k1, l, m, n, nin,      &
                             nout, x0, x1
      Use nag_library, Only: d02tgf, e02akf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: emax, x, xinc
      Integer                          :: i, ia1, ifail, iter, j, k, kp, ldc,  &
                                          liw, lsum, lw, mimax
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:,:), w(:), y(:)
      Integer, Allocatable             :: iw(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, max, real, sum
!     .. Executable Statements ..
      Write (nout,*) 'D02TGF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) mimax, kp
      Read (nin,*) x0, x1

      lsum = sum(l(1:n))
      k1 = mimax + 1
      ldc = k1
      liw = n*k1
      lw = 2*(n*kp+lsum)*(n*k1+1) + 7*n*k1
      Allocate (b(k1,n),c(ldc,n),w(lw),y(n),iw(liw))

!     Initialize coefficients b(:,:) such that z1 = 0 and z2 = 3.
      b(1:k1,1:n) = 0.0_nag_wp
      b(1,2) = 6.0_nag_wp

!     Iterate until coefficients of linearized systems converge.
      iter = 0
      emax = 1.0_nag_wp
iters: Do While (emax>=coeff_tol)
        iter = iter + 1
        Write (nout,*)
        Write (nout,99999) ' Iteration', iter, ' Chebyshev coefficients are'

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d02tgf(n,m,l,x0,x1,k1,kp,c,ldc,coeff,bdyc,w,lw,iw,liw,ifail)

        Write (nout,99998)(c(1:k1,j),j=1,n)
        emax = 0.0_nag_wp
        Do j = 1, n
          Do i = 1, k1
            emax = max(emax,abs(c(i,j)-b(i,j)))
          End Do
        End Do
        b(1:k1,1:n) = c(1:k1,1:n)
      End Do iters

      Deallocate (b)

!     Print solution on uniform mesh.
      k = 9
      ia1 = 1
      Write (nout,*)
      Write (nout,99999) 'Solution evaluated at', k, '  equally spaced points'
      Write (nout,*)
      Write (nout,99997) '      X ', (j,j=1,n)
      xinc = (x1-x0)/real(k-1,kind=nag_wp)
      x = x0
      Do i = 1, k
        Do j = 1, n
          ifail = 0
          Call e02akf(k1,x0,x1,c(1,j),ia1,k1,x,y(j),ifail)
        End Do
        Write (nout,99996) x, (y(j),j=1,n)
        x = x + xinc
      End Do

99999 Format (1X,A,I3,A)
99998 Format (1X,9F8.4)
99997 Format (1X,A,2('      Y(',I1,')'))
99996 Format (1X,3F10.4)
    End Program d02tgfe