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

    MODULE d02ndfe_mod

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

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iset = 1, neq = 3
       INTEGER, PARAMETER                  :: nia = neq + 1
       INTEGER, PARAMETER                  :: nin = 5
       INTEGER, PARAMETER                  :: njcpvt = 20*neq + 100
       INTEGER, PARAMETER                  :: nout = 6
       INTEGER, PARAMETER                  :: nrw = 50 + 4*neq
       INTEGER, PARAMETER                  :: nwkjac = 4*neq + 100
       INTEGER, PARAMETER                  :: sdysav = 6
       INTEGER, PARAMETER                  :: ldysav = neq
!      .. Local Scalars ..
       REAL (KIND=nag_wp), SAVE            :: xout
    CONTAINS
!dec$ attributes offload : mic :: fcn
       SUBROUTINE fcn(neq,t,y,f,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER       :: alpha = 0.04_nag_wp
          REAL (KIND=nag_wp), PARAMETER       :: beta = 1.0E4_nag_wp
          REAL (KIND=nag_wp), PARAMETER       :: gamma = 3.0E7_nag_wp
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq)
!         .. Executable Statements ..
          f(1) = -alpha*y(1) + beta*y(2)*y(3)
          f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
          f(3) = gamma*y(2)*y(2)
          RETURN
!dec$ attributes offload : mic :: fcn
       END SUBROUTINE fcn
!dec$ attributes offload : mic :: monitr
       SUBROUTINE monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon, &
             inln,hmin,hmax,nqu)

!         .. Use Statements ..
          USE nag_library, ONLY : d02xkf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: hlast, t
          REAL (KIND=nag_wp), INTENT (INOUT)  :: hmax, hmin, hnext
          INTEGER, INTENT (INOUT)             :: imon
          INTEGER, INTENT (OUT)               :: inln
          INTEGER, INTENT (IN)                :: ldysav, neq, nqu
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: acor(neq,2), r(neq),          &
                                                 ydot(neq), ysav(ldysav,*)
          REAL (KIND=nag_wp), INTENT (INOUT)  :: y(neq)
!         .. Local Scalars ..
          INTEGER                             :: i, ifail
!         .. Local Arrays ..
          REAL (KIND=nag_wp), ALLOCATABLE     :: z(:)
!         .. Executable Statements ..
          inln = 3

          IF (imon==1) THEN
             ALLOCATE (z(neq))

LOOP:        DO WHILE (t-hlast<xout .AND. xout<=t)

!               C1 interpolation

                ifail = 1
                CALL d02xkf(xout,z,neq,ysav,ldysav,sdysav,acor(1,2),neq,t,nqu, &
                   hlast,hnext,ifail)

                IF (ifail/=0) THEN
                   imon = -2
                   EXIT LOOP
                END IF

                WRITE (nout,99999) xout, (z(i),i=1,neq)
                xout = xout + 2.0E0_nag_wp

                IF (xout>=10.0E0_nag_wp) THEN
                   EXIT LOOP
                END IF

             END DO LOOP

          END IF

          DEALLOCATE (z)
          RETURN

99999     FORMAT (1X,F8.3,3(F13.5,2X))
!dec$ attributes offload : mic :: monitr
       END SUBROUTINE monitr
    END MODULE d02ndfe_mod

    PROGRAM d02ndfe

!      D02NDF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02ndf, d02ndz, d02nuf, d02nvf, d02nxf, d02nyf, &
                               x04abf
       USE d02ndfe_mod, ONLY : fcn, iset, ldysav, monitr, nag_wp, neq, nia,    &
                               nin, njcpvt, nout, nrw, nwkjac, sdysav, xout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: eta, h, h0, hmax, hmin, hu,      &
                                              sens, t, tcrit, tcur, tinit,     &
                                              tolsf, tout, u
       INTEGER                             :: i, icall, icase, ifail, igrow,   &
                                              imxer, isplit, isplt, itask,     &
                                              itol, itrace, liwreq, liwusd,    &
                                              lrwreq, lrwusd, maxord, maxstp,  &
                                              mxhnil, nblock, ngp, niter, nja, &
                                              nje, nlu, nnz, nq, nqu, nre,     &
                                              nst, outchn
       LOGICAL                             :: lblock, petzld
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), rtol(:), rwork(:),      &
                                              wkjac(:), y(:), ydot(:),         &
                                              yinit(:), ysav(:,:)
       REAL (KIND=nag_wp)                  :: con(6)
       INTEGER, ALLOCATABLE                :: ia(:), ja(:), jacpvt(:)
       INTEGER                             :: inform(23)
       LOGICAL, ALLOCATABLE                :: algequ(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NDF Example Program Results'
!      Skip heading in data file
       READ (nin,*)

!      Allocations based on number of differential equations (neq)
       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
          yinit(neq),ydot(neq),ysav(ldysav,sdysav),ia(nia),jacpvt(njcpvt), &
          algequ(neq))

       READ (nin,*) maxord, maxstp, mxhnil
       READ (nin,*) ia(1:nia)

       nja = ia(nia) - 1
       ALLOCATE (ja(nja))
       READ (nin,*) ja(1:nja)

!      Read algorithmic parameters
       READ (nin,*) hmin, hmax, h0, tcrit
       READ (nin,*) eta, sens, u
       READ (nin,*) lblock, petzld
       READ (nin,*) tinit, tout
       READ (nin,*) itol, isplt
       READ (nin,*) yinit(1:neq)
       READ (nin,*) rtol(1), atol(1)

       outchn = nout
       CALL x04abf(iset,outchn)

       DO icase = 1, 2
          t = tinit
          isplit = isplt
          y(1:neq) = yinit(1:neq)

!         In both cases: integrate to tout by overshooting (itask=1) using BDF
!         with Newton; use default con values, scalar tolerances and numerical
!         Jacobian; interpolate using MONITR and D02XKF.

          con(1:6) = 0.0_nag_wp
          itask = 1
          ifail = 0
          CALL d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax, &
             h0,maxstp,mxhnil,'Average-L2',rwork,ifail)
          WRITE (nout,*)

          SELECT CASE (icase)
          CASE (1)
!            No Jacobian Structiure Supplied.
             ifail = 0
             CALL d02nuf(neq,neq,'Numerical',nwkjac,ia,nia,ja,nja,jacpvt, &
                njcpvt,sens,u,eta,lblock,isplit,rwork,ifail)
             WRITE (nout,*) '  Numerical Jacobian, structure not supplied'
          CASE (2)
!            Jacobian Structiure Supplied.
             ifail = 0
             CALL d02nuf(neq,neq,'Structural',nwkjac,ia,nia,ja,nja,jacpvt, &
                njcpvt,sens,u,eta,lblock,isplit,rwork,ifail)
             WRITE (nout,*) '  Numerical Jacobian, structure supplied'
          END SELECT

          WRITE (nout,99988) (i,i=1,neq)
          WRITE (nout,99999) t, (y(i),i=1,neq)
          xout = 2.0_nag_wp

!         Soft fail and error messages only
          itrace = 0

          ifail = -1
          CALL d02ndf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
             fcn,ysav,sdysav,d02ndz,wkjac,nwkjac,jacpvt,njcpvt,monitr,itask, &
             itrace,ifail)

          IF (ifail==0) THEN

             ifail = 0
             CALL d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq, &
                niter,imxer,algequ,inform,ifail)

             WRITE (nout,*)
             WRITE (nout,99997) hu, h, tcur
             WRITE (nout,99996) nst, nre, nje
             WRITE (nout,99995) nqu, nq, niter
             WRITE (nout,99994) ' Max err comp = ', imxer
             WRITE (nout,*)

             icall = 0
             CALL d02nxf(icall,liwreq,liwusd,lrwreq,lrwusd,nlu,nnz,ngp,isplit, &
                igrow,lblock,nblock,inform)

             WRITE (nout,*)
             WRITE (nout,99993) liwreq, liwusd
             WRITE (nout,99992) lrwreq, lrwusd
             WRITE (nout,99991) nlu, nnz
             WRITE (nout,99990) ngp, isplit
             WRITE (nout,99989) igrow, nblock
          ELSE IF (ifail==10) THEN
             WRITE (nout,*)
             WRITE (nout,99998) 'Exit D02NDF with IFAIL = ', ifail, &
                '  and T = ', t

             icall = 1
             CALL d02nxf(icall,liwreq,liwusd,lrwreq,lrwusd,nlu,nnz,ngp,isplit, &
                igrow,lblock,nblock,inform)

             WRITE (nout,*)
             WRITE (nout,99993) liwreq, liwusd
             WRITE (nout,99992) lrwreq, lrwusd
          ELSE
             WRITE (nout,*)
             WRITE (nout,99998) 'Exit D02NDF with IFAIL = ', ifail, &
                '  and T = ', t
          END IF
       END DO

99999  FORMAT (1X,F8.3,3(F13.5,2X))
99998  FORMAT (1X,A,I5,A,E12.5)
99997  FORMAT (1X,' HUSED = ',E12.5,'  HNEXT = ',E12.5,'  TCUR = ',E12.5)
99996  FORMAT (1X,' NST = ',I6,'    NRE = ',I6,'    NJE = ',I6)
99995  FORMAT (1X,' NQU = ',I6,'    NQ  = ',I6,'  NITER = ',I6)
99994  FORMAT (1X,A,I4)
99993  FORMAT (1X,' NJCPVT (required ',I4,'  used ',I8,')')
99992  FORMAT (1X,' NWKJAC (required ',I4,'  used ',I8,')')
99991  FORMAT (1X,' No. of LU-decomps ',I4,'  No. of nonzeros ',I8)
99990  FORMAT (1X,' No. of FCN calls to form Jacobian ',I4,'  Try ISPLIT ',I4)
99989  FORMAT (1X,' Growth est ',I8,'  No. of blocks on diagonal ',I4)
99988  FORMAT (/1X,'    X ',3('         Y(',I1,')  '))
    END PROGRAM d02ndfe
