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

    MODULE d03phfe_mod

!      D03PHF 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       :: one = 1.0_nag_wp
       INTEGER, PARAMETER                  :: itrace = 0, ncode = 1, nin = 5,  &
                                              nout = 6, npde = 1, nxi = 1
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: ts
    CONTAINS
!dec$ attributes offload : mic :: odedef
       SUBROUTINE odedef(npde,t,ncode,v,vdot,nxi,xi,ucp,ucpx,rcp,ucpt,ucptx,f, &
             ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: ncode, npde, nxi
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(ncode)
          REAL (KIND=nag_wp), INTENT (IN)     :: rcp(npde,*), ucp(npde,*),     &
                                                 ucpt(npde,*), ucptx(npde,*),  &
                                                 ucpx(npde,*), v(ncode),       &
                                                 vdot(ncode), xi(nxi)
!         .. Executable Statements ..
          IF (ires==1) THEN
             f(1) = vdot(1) - v(1)*ucp(1,1) - ucpx(1,1) - one - t
          ELSE IF (ires==-1) THEN
             f(1) = vdot(1)
          END IF
          RETURN
!dec$ attributes offload : mic :: odedef
       END SUBROUTINE odedef
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,ux,ncode,v,vdot,p,q,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t, x
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: ncode, npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: p(npde,npde), q(npde), r(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
!         .. Executable Statements ..
          p(1,1) = v(1)*v(1)
          r(1) = ux(1)
          q(1) = -x*ux(1)*v(1)*vdot(1)
          RETURN
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndary
       SUBROUTINE bndary(npde,t,u,ux,ncode,v,vdot,ibnd,beta,gamma,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, ncode, npde
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: beta(npde), gamma(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          beta(1) = one
          IF (ibnd==0) THEN
             gamma(1) = -v(1)*exp(t)
          ELSE
             gamma(1) = -v(1)*vdot(1)
          END IF
          RETURN
!dec$ attributes offload : mic :: bndary
       END SUBROUTINE bndary
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,x,u,ncode,neqn)

!         Routine for PDE initial values

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ncode, neqn, npde, npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(neqn)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          DO i = 1, npts
             u(i) = exp(ts*(one-x(i))) - one
          END DO
          u(neqn) = ts
          RETURN
!dec$ attributes offload : mic :: uvinit
       END SUBROUTINE uvinit
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(time,npts,x,u)
!         Exact solution (for comparison purpose)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: time
          INTEGER, INTENT (IN)                :: npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          DO i = 1, npts
             u(i) = exp(time*(one-x(i))) - one
          END DO
          RETURN
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
    END MODULE d03phfe_mod
    PROGRAM d03phfe

!      D03PHF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03phf
       USE d03phfe_mod, ONLY : bndary, exact, itrace, nag_wp, ncode, nin,      &
                               nout, npde, nxi, odedef, pdedef, ts, uvinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: tout
       INTEGER                             :: i, ifail, ind, it, itask, itol,  &
                                              latol, lenode, lisave, lrsave,   &
                                              lrtol, m, neqn, npts, nwkres
       LOGICAL                             :: theta
       CHARACTER (1)                       :: laopt, norm
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: algopt(30), xi(nxi)
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), exy(:), rsave(:),       &
                                              rtol(:), u(:), x(:)
       INTEGER, ALLOCATABLE                :: isave(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              mod, real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PHF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) m, npts
       neqn = npde*npts + ncode
       nwkres = npde*(npts+6*nxi+3*npde+15) + ncode + nxi + 7*npts + 2
       lenode = 11*neqn + 50
       lrsave = neqn*neqn + neqn + nwkres + lenode
       lisave = 25*neqn + 24

       ALLOCATE (exy(npts),u(neqn),rsave(lrsave),x(npts),isave(lisave))

       READ (nin,*) itol
       latol = 1
       lrtol = 1
       IF (itol>2) latol = neqn
       IF (mod(itol,2)==0) lrtol = neqn
       ALLOCATE (atol(latol),rtol(lrtol))
       READ (nin,*) atol(1:latol), rtol(1:lrtol)

       READ (nin,*) ts

!      Set break-points
       DO i = 1, npts
          x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
       END DO

       READ (nin,*) xi(1:nxi)
       READ (nin,*) norm, laopt
       ind = 0
       itask = 1

!      Set theta to .TRUE. if the Theta integrator is required
       theta = .FALSE.
       algopt(1:30) = 0.0_nag_wp
       IF (theta) THEN
          algopt(1) = 2.0_nag_wp
       END IF

!      Loop over output value of t

       CALL uvinit(npde,npts,x,u,ncode,neqn)
       tout = 0.2_nag_wp
       DO it = 1, 5

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          CALL d03phf(npde,m,ts,tout,pdedef,bndary,u,npts,x,ncode,odedef,nxi, &
             xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
             lisave,itask,itrace,ind,ifail)

          IF (it==1) THEN
             WRITE (nout,99997) atol, npts
             WRITE (nout,99999) x(1:npts-5:4), x(npts)
          END IF

!         Print against the exact solution.
          CALL exact(tout,npts,x,exy)
          WRITE (nout,99998) ts
          WRITE (nout,99995) u(1:npts-5:4), u(npts:neqn)
          WRITE (nout,99994) exy(1:npts-5:4), exy(npts), ts
!         Select next time to solve to for output.
          tout = 2.0_nag_wp*tout
       END DO
       WRITE (nout,99996) isave(1), isave(2), isave(3), isave(5)

99999  FORMAT ('  X        ',5F9.3/)
99998  FORMAT (' T = ',F6.3)
99997  FORMAT (//'  Simple coupled PDE using BDF '/'  Accuracy require', &
          'ment =',E10.3,' Number of points = ',I4/)
99996  FORMAT (' Number of integration steps in time = ',I6/' Number o', &
          'f function evaluations = ',I6/' Number of Jacobian eval', &
          'uations =',I6/' Number of iterations = ',I6)
99995  FORMAT (1X,'App.  sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3)
99994  FORMAT (1X,'Exact sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3/)
    END PROGRAM d03phfe
