!   D02NJF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    Module d02njfe_mod

!     D02NJF 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                           :: jac, monitr, resid
!     .. 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
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: iset = 1, itrace = 0, neq = 3
      Integer, Parameter, Public       :: nia = neq + 1
      Integer, Parameter, Public       :: nin = 5, nout = 6
      Integer, Parameter, Public       :: nrw = 50 + 4*neq
      Integer, Parameter, Public       :: ldysav = neq
      Integer, Parameter               :: nelts = 8
      Integer, Parameter, Public       :: nja = nelts
      Integer, Parameter, Public       :: njcpvt = 20*neq + 12*nelts
      Integer, Parameter, Public       :: nwkjac = 4*neq + 12*nelts
    Contains
      Subroutine resid(neq,t,y,ydot,r,ires)

!       .. 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) :: r(neq)
        Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
!       .. Executable Statements ..
        r(1) = zero
        r(2) = -ydot(2)
        r(3) = -ydot(3)
        If (ires==1) Then
          r(1) = y(1) + y(2) + y(3) - one + r(1)
          r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2)
          r(3) = gamma*y(2)*y(2) + r(3)
        End If
        Return
      End Subroutine resid

      Subroutine jac(neq,t,y,ydot,h,d,j,pdj)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: d, h, t
        Integer, Intent (In)           :: j, neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: pdj(neq)
        Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: hxd
!       .. Executable Statements ..
!       8 nonzero elements in total
        hxd = h*d
        If (j==1) Then
          pdj(1) = zero - hxd*(one)
          pdj(2) = zero - hxd*(alpha)
!         note: pdj(3) is zero
        Else If (j==2) Then
          pdj(1) = zero - hxd*(one)
          pdj(2) = one - hxd*(-beta*y(3)-two*gamma*y(2))
          pdj(3) = zero - hxd*(two*gamma*y(2))
        Else If (j==3) Then
          pdj(1) = zero - hxd*(one)
          pdj(2) = zero - hxd*(-beta*y(2))
          pdj(3) = one - hxd*(zero)
        End If
        Return
      End Subroutine jac

      Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
        hmin,hmax,nqu)

!       .. 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)
!       .. Executable Statements ..
        inln = 3
        If (y(1)<=0.9_nag_wp) Then
          imon = -2
        End If
        Return
      End Subroutine monitr
    End Module d02njfe_mod

    Program d02njfe

!     D02NJF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02njf, d02nuf, d02nvf, d02nxf, d02nyf, nag_wp,   &
                             x04abf
      Use d02njfe_mod, Only: iset, itrace, jac, ldysav, monitr, neq, nia, nin, &
                             nja, njcpvt, nout, nrw, nwkjac, resid
!     .. 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,   &
                                          liwreq, liwusd, lrwreq, lrwusd,      &
                                          maxord, maxstp, mxhnil, nblock, ngp, &
                                          niter, nje, nlu, nnz, nq, nqu, nre,  &
                                          nst, outchn, sdysav
      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(:)
      Logical                          :: lderiv(2)
!     .. Executable Statements ..
      Write (nout,*) 'D02NJF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) maxord, maxstp, mxhnil
      sdysav = maxord + 1
      Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),           &
        yinit(neq),ydot(neq),ysav(ldysav,sdysav),ia(nia),ja(nja),              &
        jacpvt(njcpvt),algequ(neq))
      Read (nin,*) ia(1:nia)
      Read (nin,*) ja(1:nja)

      outchn = nout
      Call x04abf(iset,outchn)

!     Two cases. In both cases:
!          integrate to tout by overshooting (itask=1);
!          use B.D.F formulae with a Newton method;
!          use the Petzold error test (differential algebraic system);
!          use default values for the array con;
!          employ vector relative tolerance and scalar absolute tolerance.
!          the Jacobian is supplied by jac;
!          the monitr routine is used to force a return when y(1) < 0.9.

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

      con(1:6) = 0.0_nag_wp
      itask = 1
      petzld = .True.

cases: Do icase = 1, 2

!       Initialize
        t = tinit
        isplit = isplt
        y(1:neq) = yinit(1:neq)
        lderiv(1:2) = .False.

        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)
!         First case. The Jacobian structure is determined internally by
!                     calls to jac.
          ifail = 0
          Call d02nuf(neq,neq,'Analytical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
            sens,u,eta,lblock,isplit,rwork,ifail)
          Write (nout,*) '  Analytic Jacobian, structure not supplied'
        Case (2)
!         Second case. The Jacobian structure is supplied.
          ifail = 0
          Call d02nuf(neq,neq,'Full info',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt,  &
            sens,u,eta,lblock,isplit,rwork,ifail)
          Write (nout,*) '  Analytic 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

        ifail = 1
        Call d02njf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,      &
          resid,ysav,sdysav,jac,wkjac,nwkjac,jacpvt,njcpvt,monitr,lderiv,      &
          itask,itrace,ifail)

        If (ifail==0 .Or. ifail==12) Then
          Write (nout,99999) t, (y(i),i=1,neq)

          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
          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
          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 D02NJF with IFAIL = ', ifail, '  and T = ', &
            t
        End If
      End Do cases

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 d02njfe