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

    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
!      .. 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                  :: iset = 1, itrace = 0, nelts = 8, &
                                              neq = 3
       INTEGER, PARAMETER                  :: nia = neq + 1
       INTEGER, PARAMETER                  :: nin = 5
       INTEGER, PARAMETER                  :: nja = nelts
       INTEGER, PARAMETER                  :: njcpvt = 20*neq + 12*nelts
       INTEGER, PARAMETER                  :: nout = 6
       INTEGER, PARAMETER                  :: nrw = 50 + 4*neq
       INTEGER, PARAMETER                  :: nwkjac = 4*neq + 12*nelts
       INTEGER, PARAMETER                  :: ldysav = neq
    CONTAINS
!dec$ attributes offload : mic :: resid
       SUBROUTINE resid(neq,t,y,ydot,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: resid
       END SUBROUTINE resid

!dec$ attributes offload : mic :: jac
       SUBROUTINE jac(neq,t,y,ydot,h,d,j,pdj)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 non-zero 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
!dec$ attributes offload : mic :: jac
       END SUBROUTINE jac

!dec$ attributes offload : mic :: monitr
       SUBROUTINE monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon, &
             inln,hmin,hmax,nqu)

!         .. 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)
!         .. Executable Statements ..
          inln = 3
          IF (y(1)<=0.9_nag_wp) imon = -2
          RETURN
!dec$ attributes offload : mic :: monitr
       END SUBROUTINE monitr
    END MODULE d02njfe_mod

    PROGRAM d02njfe

!      D02NJF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02njf, d02nuf, d02nvf, d02nxf, d02nyf, x04abf
       USE d02njfe_mod, ONLY : iset, itrace, jac, ldysav, monitr, nag_wp, 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 is used (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
