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

    MODULE d02tyfe_mod

!      D02TYF 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, nmesh_out = 11,        &
                                              nout = 6, nrbc = 1
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: a
!      .. 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)
!         .. Intrinsic Functions ..
          INTRINSIC                              sqrt
!         .. Executable Statements ..
          IF (y(1,0)<=0.0E0_nag_wp) THEN
             f(1) = 0.0_nag_wp
          ELSE
             f(1) = (y(1,0))**1.5_nag_wp/sqrt(x)
          END IF
          RETURN
!dec$ attributes offload : mic :: ffun
       END SUBROUTINE ffun
!dec$ attributes offload : mic :: fjac
       SUBROUTINE fjac(x,y,neq,m,dfdy)

!         .. 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)
!         .. Intrinsic Functions ..
          INTRINSIC                              sqrt
!         .. Executable Statements ..
          IF (y(1,0)<=0.0E0_nag_wp) THEN
             dfdy(1,1,0) = 0.0_nag_wp
          ELSE
             dfdy(1,1,0) = 1.5_nag_wp*sqrt(y(1,0))/sqrt(x)
          END IF
          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) - 1.0_nag_wp
          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)
          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
!         .. 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) = 1.0_nag_wp
          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
!         .. 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) = 1.0_nag_wp
          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) = 1.0_nag_wp - x/a
          y(1,1) = -1.0_nag_wp/a
          dym(1) = 0.0_nag_wp
          RETURN
!dec$ attributes offload : mic :: guess
       END SUBROUTINE guess
    END MODULE d02tyfe_mod
    PROGRAM d02tyfe

!      D02TYF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02tkf, d02tvf, d02tyf, d02tzf
       USE d02tyfe_mod, ONLY : a, ffun, fjac, gafun, gajac, gbfun, gbjac,      &
                               guess, m, mmax, nag_wp, neq, nin, nlbc,         &
                               nmesh_out, nout, nrbc
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: ainc, ermx, x
       INTEGER                             :: i, iermx, ifail, ijermx, liwork, &
                                              lrwork, mxmesh, ncol, nmesh
       LOGICAL                             :: failed
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: mesh(:), rwork(:), tol(:), y(:,:)
       INTEGER, ALLOCATABLE                :: ipmesh(:), iwork(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D02TYF 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,*) a
       READ (nin,*) tol(1:neq)

       ainc = a/real(nmesh-1,kind=nag_wp)
       mesh(1) = 0.0E0_nag_wp
       DO i = 2, nmesh - 1
          mesh(i) = mesh(i-1) + ainc
       END DO
       mesh(nmesh) = a

       ipmesh(1) = 1
       ipmesh(2:nmesh-1) = 2
       ipmesh(nmesh) = 1

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

       WRITE (nout,99999) tol(1), a

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

       IF (ifail/=1) THEN
!         Print mesh statistics
          WRITE (nout,99998) nmesh, ermx, iermx, ijermx
          WRITE (nout,99997) (i,ipmesh(i),mesh(i),i=1,nmesh)
       END IF
       IF ( .NOT. failed) THEN
!         Print solution on output mesh.
          WRITE (nout,99996)
          x = 0.0_nag_wp
          ainc = a/real(nmesh_out-1,kind=nag_wp)
          DO i = 1, nmesh_out
             ifail = 0
             CALL d02tyf(x,y,neq,mmax,rwork,iwork,ifail)
             WRITE (nout,99995) x, y(1,0:1)
             x = x + ainc
          END DO
       END IF

99999  FORMAT (//' Tolerance = ',E8.1,' A = ',F8.2)
99998  FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
          '  in interval ',I4,' for component ',I4/)
99997  FORMAT (/' Mesh points:'/4(I4,'(',I1,')',E11.4))
99996  FORMAT (/' Computed solution'/'       x     solution   derivative')
99995  FORMAT (' ',F8.2,2F11.5)

    END PROGRAM d02tyfe
