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

    MODULE d02txfe_mod

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

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: half = 0.5_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: three = 3.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: xsplit = 30.0_nag_wp
       INTEGER, PARAMETER                  :: m1 = 3, m2 = 2, mmax = 3,        &
                                              neq = 2, nin = 5, nlbc = 3,      &
                                              nleft = 15, nout = 6, nrbc = 2,  &
                                              nright = 10
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: el, en, s
    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)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: t1, y11, y20
!         .. Executable Statements ..
          t1 = half*(three-en)*y(1,0)
          y11 = y(1,1)
          y20 = y(2,0)
          f(1) = (el**3)*(one-y20**2) + (el**2)*s*y11 - &
             el*(t1*y(1,2)+en*y11**2)
          f(2) = (el**2)*s*(y20-one) - el*(t1*y(2,1)+(en-one)*y11*y20)
          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)
!         .. Executable Statements ..
          dfdy(1,2,0) = -two*el**3*y(2,0)
          dfdy(1,1,0) = -el*half*(three-en)*y(1,2)
          dfdy(1,1,1) = el**2*s - el*two*en*y(1,1)
          dfdy(1,1,2) = -el*half*(three-en)*y(1,0)
          dfdy(2,2,0) = el**2*s - el*(en-one)*y(1,1)
          dfdy(2,2,1) = -el*half*(three-en)*y(1,0)
          dfdy(2,1,0) = -el*half*(three-en)*y(2,1)
          dfdy(2,1,1) = -el*(en-one)*y(2,0)
          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)
          ga(2) = ya(1,1)
          ga(3) = ya(2,0)
          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,1)
          gb(2) = yb(2,0) - one
          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) = one
          dgady(2,1,1) = one
          dgady(3,2,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
!         .. 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,1) = one
          dgbdy(2,2,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)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: ex, expmx
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          ex = x*el
          expmx = exp(-ex)
          y(1,0) = -ex**2*expmx
          y(1,1) = (-two*ex+ex**2)*expmx
          y(1,2) = (-two+4.0E0_nag_wp*ex-ex**2)*expmx
          y(2,0) = one - expmx
          y(2,1) = expmx
          dym(1) = (6.0E0_nag_wp-6.0E0_nag_wp*ex+ex**2)*expmx
          dym(2) = -expmx
          RETURN
!dec$ attributes offload : mic :: guess
       END SUBROUTINE guess
    END MODULE d02txfe_mod
    PROGRAM d02txfe

!      D02TXF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf
       USE d02txfe_mod, ONLY : el, en, ffun, fjac, gafun, gajac, gbfun, gbjac, &
                               guess, m1, m2, mmax, nag_wp, neq, nin, nlbc,    &
                               nleft, nout, nrbc, nright, one, s, two, xsplit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: dx, el_init, ermx, s_init, xx
       INTEGER                             :: i, iermx, ifail, ijermx, j,      &
                                              liwork, lrwork, mxmesh, ncol,    &
                                              ncont, nmesh
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: mesh(:), rwork(:)
       REAL (KIND=nag_wp)                  :: tol(neq), y(neq,0:mmax-1)
       INTEGER, ALLOCATABLE                :: ipmesh(:), iwork(:)
       INTEGER                             :: m(neq)
!      .. Intrinsic Functions ..
       INTRINSIC                              min, real
!      .. Executable Statements ..
       WRITE (nout,*) 'D02TXF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
!      Read method parameters
       READ (nin,*) ncol, nmesh, mxmesh
       READ (nin,*) tol(1:neq)

       liwork = mxmesh*(11*neq+6)
       lrwork = mxmesh*(109*neq**2+78*neq+7)
       ALLOCATE (mesh(mxmesh),rwork(lrwork),ipmesh(mxmesh),iwork(liwork))
!      Read problem (initial) parameters
       READ (nin,*) en, el_init, s_init
!      Initialise data
       el = el_init
       s = s_init
       m(1) = m1
       m(2) = m2

       dx = one/real(nmesh-1,kind=nag_wp)
       mesh(1) = 0.0_nag_wp
       DO i = 2, nmesh - 1
          mesh(i) = mesh(i-1) + dx
       END DO
       mesh(nmesh) = 1.0_nag_wp

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

!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, &
          lrwork,iwork,liwork,ifail)

!      Initialize number of continuation steps in el and s
       READ (nin,*) ncont
CONT:  DO j = 1, ncont
          WRITE (nout,99997) tol(1), el, s

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

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

          WRITE (nout,99996) nmesh, ermx, iermx, ijermx
!         Print solution components on mesh
          WRITE (nout,99999)

!         Left side domain [0,xsplit], evaluate at nleft+1 uniform grid points.
          dx = xsplit/real(nleft,kind=nag_wp)/el
          xx = 0.0_nag_wp
          DO i = 0, nleft
             ifail = 0
             CALL d02tyf(xx,y,neq,mmax,rwork,iwork,ifail)
             WRITE (nout,99998) xx*el, y(1,0), y(2,0)
             xx = xx + dx
          END DO

!         Right side domain (xsplit,L], evaluate at nright uniform grid points.
          dx = (el-xsplit)/real(nright,kind=nag_wp)/el
          xx = xsplit/el
          DO i = 1, nright
             xx = min(1.0_nag_wp,xx+dx)
             ifail = 0
             CALL d02tyf(xx,y,neq,mmax,rwork,iwork,ifail)
             WRITE (nout,99998) xx*el, y(1,0), y(2,0)
          END DO

!         Select mesh for continuation and update continuation parameters.
          IF (j<ncont) THEN
             el = two*el
             s = 0.6_nag_wp*s
             nmesh = (nmesh+1)/2
             ifail = 0
             CALL d02txf(mxmesh,nmesh,mesh,ipmesh,rwork,iwork,ifail)
          END IF
       END DO CONT

99999  FORMAT (/' Solution on original interval:'/6X,'x',8X,'f',10X,'g')
99998  FORMAT (1X,F8.2,2(1X,F10.4))
99997  FORMAT (//' Tolerance = ',E8.1,'  L = ',F8.3,'  S =',F7.4)
99996  FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
          '  in interval ',I4,' for component ',I4)
    END PROGRAM d02txfe
