!   D02NDF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn, monitr
!     .. Parameters ..
      Integer, Parameter, Public       :: iset = 1, neq = 3
      Integer, Parameter, Public       :: nia = neq + 1
      Integer, Parameter, Public       :: nin = 5
      Integer, Parameter, Public       :: njcpvt = 20*neq + 100
      Integer, Parameter, Public       :: nout = 6
      Integer, Parameter, Public       :: nrw = 50 + 4*neq
      Integer, Parameter, Public       :: nwkjac = 4*neq + 100
      Integer, Parameter, Public       :: sdysav = 6
      Integer, Parameter, Public       :: ldysav = neq
    Contains
      Subroutine fcn(neq,t,y,f,ires)

!       .. 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
      End Subroutine fcn
      Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
        hmin,hmax,nqu)

!       .. Use Statements ..
        Use nag_library, Only: d02xkf
!       .. 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 ..
        Real (Kind=nag_wp)             :: xout
        Integer                        :: i, ifail
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: z(:)
!       .. Executable Statements ..
        inln = 3

        If (imon==1) Then
          Allocate (z(neq))

!         Interpolate at multiples of 2.0 between t-hlast and t
          xout = 2.0E0_nag_wp
          Do While (xout<=t-hlast)
            xout = xout + 2.0E0_nag_wp
          End Do
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))
      End Subroutine monitr
    End Module d02ndfe_mod

    Program d02ndfe

!     D02NDF Example Main Program

!     .. Use Statements ..
      Use d02ndfe_mod, Only: fcn, iset, ldysav, monitr, neq, nia, nin, njcpvt, &
                             nout, nrw, nwkjac, sdysav
      Use nag_library, Only: d02ndf, d02ndz, d02nuf, d02nvf, d02nxf, d02nyf,   &
                             nag_wp, x04abf
!     .. 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 Structure 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 Structure 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)

!       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