!   D02TZF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module d02tzfe_mod

!     D02TZF 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                           :: ffun, fjac, gafun, gajac, gbfun,     &
                                          gbjac, guess
!     .. Parameters ..
      Integer, Parameter, Public       :: mmax = 2, neq = 1, nin = 5,          &
                                          nlbc = 1, nout = 6, nrbc = 1
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: alpha, beta, eps
!     .. Local Arrays ..
      Integer, Public, Save            :: m(1) = (/2/)
    Contains
      Subroutine ffun(x,y,neq,m,f,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(neq)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        f(1) = (y(1,0)-y(1,0)*y(1,1))/eps
        Return
      End Subroutine ffun
      Subroutine fjac(x,y,neq,m,dfdy,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: x02ajf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dfdy(neq,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: epsh, fac, ptrb
        Integer                        :: i, j, k
!       .. Local Arrays ..
        Real (Kind=nag_wp)             :: f1(1), f2(1), yp(1,0:3)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs, max, sqrt
!       .. Executable Statements ..
        epsh = 100.0_nag_wp*x02ajf()
        fac = sqrt(x02ajf())
        Do i = 1, neq
          Do j = 0, m(i) - 1
            yp(i,j) = y(i,j)
          End Do
        End Do
        Do i = 1, neq
          Do j = 0, m(i) - 1
            ptrb = max(epsh,fac*abs(y(i,j)))
            yp(i,j) = y(i,j) + ptrb
            Call ffun(x,yp,neq,m,f1,iuser,ruser)
            yp(i,j) = y(i,j) - ptrb
            Call ffun(x,yp,neq,m,f2,iuser,ruser)
            Do k = 1, neq
              dfdy(k,i,j) = 0.5_nag_wp*(f1(k)-f2(k))/ptrb
            End Do
            yp(i,j) = y(i,j)
          End Do
        End Do
        Return
      End Subroutine fjac
      Subroutine gafun(ya,neq,m,nlbc,ga,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: ga(nlbc)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        ga(1) = ya(1,0) - alpha
        Return
      End Subroutine gafun
      Subroutine gbfun(yb,neq,m,nrbc,gb,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: gb(nrbc)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        gb(1) = yb(1,0) - beta
        Return
      End Subroutine gbfun
      Subroutine gajac(ya,neq,m,nlbc,dgady,iuser,ruser)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: one = 1.0_nag_wp
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dgady(nlbc,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        dgady(1,1,0) = one
        Return
      End Subroutine gajac
      Subroutine gbjac(yb,neq,m,nrbc,dgbdy,iuser,ruser)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: one = 1.0_nag_wp
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: dgbdy(nrbc,neq,0:*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        dgbdy(1,1,0) = one
        Return
      End Subroutine gbjac
      Subroutine guess(x,neq,m,y,dym,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: dym(neq)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*), y(neq,0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: m(neq)
!       .. Executable Statements ..
        y(1,0) = alpha + (beta-alpha)*x
        y(1,1) = (beta-alpha)
        dym(1) = 0.0_nag_wp
        Return
      End Subroutine guess
    End Module d02tzfe_mod
    Program d02tzfe

!     D02TZF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02tlf, d02tvf, d02txf, d02tyf, d02tzf, nag_wp
      Use d02tzfe_mod, Only: alpha, beta, eps, ffun, fjac, gafun, gajac,       &
                             gbfun, gbjac, guess, m, mmax, neq, nin, nlbc,     &
                             nout, nrbc
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: ermx
      Integer                          :: i, iermx, ifail, ijermx, j, licomm,  &
                                          lrcomm, mxmesh, ncol, nmesh
      Logical                          :: failed
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: mesh(:), rcomm(:), tol(:), y(:,:)
      Real (Kind=nag_wp)               :: ruser(1)
      Integer, Allocatable             :: icomm(:), ipmesh(:)
      Integer                          :: iuser(2)
!     .. Executable Statements ..
      Write (nout,*) 'D02TZF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) ncol, nmesh, mxmesh
      Allocate (mesh(mxmesh),tol(neq),y(neq,0:mmax-1),ipmesh(mxmesh))

      Read (nin,*) alpha, beta, eps
      Read (nin,*) mesh(1:nmesh)
      Read (nin,*) ipmesh(1:nmesh)
      Read (nin,*) tol(1:neq)

!     Workspace query to get size of rcomm and icomm
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,ruser,0,   &
        iuser,2,ifail)
      lrcomm = iuser(1)
      licomm = iuser(2)
      Allocate (rcomm(lrcomm),icomm(licomm))

!     Initialize
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rcomm,     &
        lrcomm,icomm,licomm,ifail)

      eps = 0.1_nag_wp*eps

contn: Do j = 1, 2
        Write (nout,99997) tol(1), eps

!       Solve
        ifail = -1
        Call d02tlf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rcomm,icomm,iuser, &
          ruser,ifail)
        failed = ifail /= 0

!       Extract mesh.
        ifail = -1
        Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rcomm,icomm,    &
          ifail)

!       Print mesh statistics.
        Write (nout,99996) nmesh, ermx, iermx, ijermx

        If (failed) Then
          Exit contn
        End If

!       Print solution at every second point on final mesh.
        Write (nout,99999)
        Do i = 1, nmesh, 2
          ifail = -1
          Call d02tyf(mesh(i),y,neq,mmax,rcomm,icomm,ifail)
          Write (nout,99998) mesh(i), y(1,0), y(1,1)
        End Do

        If (j==1) Then
!         Halve final mesh for new initial mesh and set up for continuation.
          nmesh = (nmesh+1)/2
          ifail = 0
          Call d02txf(mxmesh,nmesh,mesh,ipmesh,rcomm,icomm,ifail)
!         Reduce continuation parameter.
          eps = 0.1_nag_wp*eps
        End If
      End Do contn

99999 Format (/,' Solution and derivative at every second point:',/,           &
        '      x        u          u''')
99998 Format (' ',F8.4,2F11.5)
99997 Format (/,/,' Tolerance = ',E8.1,'  EPS = ',E10.3)
99996 Format (/,' Used a mesh of ',I4,' points',/,' Maximum error = ',E10.2,   &
        '  in interval ',I4,' for component ',I4)
    End Program d02tzfe