!   D02TZF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: mmax = 2, neq = 1, nin = 5,      &
                                              nlbc = 1, nout = 6, nrbc = 1
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: alpha, beta, eps
!      .. Local Arrays ..
       INTEGER                             :: m(1) = (/ 2/)
    CONTAINS
!dec$ attributes offload : mic :: ffun
       SUBROUTINE ffun(x,y,neq,m,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 (IN)     :: y(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          f(1) = (y(1,0)-y(1,0)*y(1,1))/eps
          RETURN
!dec$ attributes offload : mic :: ffun
       END SUBROUTINE ffun
!dec$ attributes offload : mic :: fjac
       SUBROUTINE fjac(x,y,neq,m,dfdy)

!         .. Use Statements ..
          USE nag_library, ONLY : x02ajf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq,0:*)
          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 Functions ..
          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)
                yp(i,j) = y(i,j) - ptrb
                CALL ffun(x,yp,neq,m,f2)
                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
!dec$ attributes offload : mic :: fjac
       END SUBROUTINE fjac
!dec$ attributes offload : mic :: gafun
       SUBROUTINE gafun(ya,neq,m,nlbc,ga)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nlbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: ga(nlbc)
          REAL (KIND=nag_wp), INTENT (IN)     :: ya(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          ga(1) = ya(1,0) - alpha
          RETURN
!dec$ attributes offload : mic :: gafun
       END SUBROUTINE gafun
!dec$ attributes offload : mic :: gbfun
       SUBROUTINE gbfun(yb,neq,m,nrbc,gb)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: neq, nrbc
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: gb(nrbc)
          REAL (KIND=nag_wp), INTENT (IN)     :: yb(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          gb(1) = yb(1,0) - beta
          RETURN
!dec$ attributes offload : mic :: gbfun
       END SUBROUTINE gbfun
!dec$ attributes offload : mic :: gajac
       SUBROUTINE gajac(ya,neq,m,nlbc,dgady)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: ya(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          dgady(1,1,0) = one
          RETURN
!dec$ attributes offload : mic :: gajac
       END SUBROUTINE gajac
!dec$ attributes offload : mic :: gbjac
       SUBROUTINE gbjac(yb,neq,m,nrbc,dgbdy)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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:*)
          REAL (KIND=nag_wp), INTENT (IN)     :: yb(neq,0:*)
          INTEGER, INTENT (IN)                :: m(neq)
!         .. Executable Statements ..
          dgbdy(1,1,0) = one
          RETURN
!dec$ attributes offload : mic :: gbjac
       END SUBROUTINE gbjac
!dec$ attributes offload : mic :: guess
       SUBROUTINE guess(x,neq,m,y,dym)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)  :: y(neq,0:*)
          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
!dec$ attributes offload : mic :: guess
       END SUBROUTINE guess
    END MODULE d02tzfe_mod
    PROGRAM d02tzfe

!      D02TZF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf
       USE d02tzfe_mod, ONLY : alpha, beta, eps, ffun, fjac, gafun, gajac,     &
                               gbfun, gbjac, guess, m, mmax, nag_wp, neq, nin, &
                               nlbc, nout, nrbc
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: ermx
       INTEGER                             :: i, iermx, ifail, ijermx, j,      &
                                              liwork, lrwork, mxmesh, ncol,    &
                                              nmesh
       LOGICAL                             :: failed
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: mesh(:), rwork(:), tol(:), y(:,:)
       INTEGER, ALLOCATABLE                :: ipmesh(:), iwork(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02TZF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) ncol, nmesh, mxmesh
       liwork = mxmesh*(11*neq+6)
       lrwork = mxmesh*(109*neq**2+78*neq+7)
       ALLOCATE (mesh(mxmesh),tol(neq),rwork(lrwork),y(neq,0:mmax-1), &
          ipmesh(mxmesh),iwork(liwork))

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

!      Initialize
       ifail = 0
       CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, &
          lrwork,iwork,liwork,ifail)

       eps = 0.1_nag_wp*eps

CONTN: DO j = 1, 2
          WRITE (nout,99997) tol(1), eps

!         Solve
          ifail = -1
          CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork, &
             ifail)
          failed = ifail /= 0

!         Extract mesh.
          ifail = -1
          CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, &
             ifail)

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

          IF (failed) EXIT CONTN

!         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,rwork,iwork,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,rwork,iwork,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
