    PROGRAM d02nnfe

!      D02NNF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : d02mzf, d02nnf, d02nrf, d02nuf, d02nvf, d02nxf, &
                               d02nyf, nag_wp, x04abf
!      .. 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   :: zero = 0.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER   :: gamm2 = 2.0_nag_wp*gamma
       INTEGER, PARAMETER              :: iset = 1, itrace = 0, nelts = 8,     &
                                          neq = 3, nia = 1, nin = 5, nja = 1
       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
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: eta, h, h0, hmax, hmin, hu, hxd,     &
                                          sens, t, tcrit, tcur, tolsf, tout, u
       INTEGER                         :: i, icall, ifail, igrow, imon, imxer, &
                                          indd, indr, inln, iplace, ires,      &
                                          irevcm, isplit, itask, itol, j,      &
                                          liwreq, liwusd, lrwreq, lrwusd,      &
                                          maxord, maxstp, mxhnil, nblock,      &
                                          nfails, 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(:), ysav(:,:)
       REAL (KIND=nag_wp)              :: con(6)
       INTEGER                         :: ia(nia), inform(23), ja(nja)
       INTEGER, ALLOCATABLE            :: jacpvt(:)
       LOGICAL, ALLOCATABLE            :: algequ(:)
       LOGICAL                         :: lderiv(2)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NNF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
!      neq: number of differential equations
       READ (nin,*) maxord, maxstp, mxhnil
       sdysav = maxord + 1
       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
          ydot(neq),ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq))

       outchn = nout
       WRITE (nout,*)
       CALL x04abf(iset,outchn)

!      Integrate towards tout stopping at the first mesh point beyond
!      tout (itask=3) using the B.D.F. formulae with a Newton method.
!      Employ scalar tolerances and the Jacobian is supplied, but its
!      structure is evaluated internally by calls to the Jacobian
!      forming part of the program (irevcm=8). Default values for the
!      array con are used. Also count the number of step failures
!      (irevcm=10). The solution is interpolated using D02MZF to give
!      the solution at tout.

       READ (nin,*) hmin, hmax, h0, tcrit
       READ (nin,*) eta, sens, u
       READ (nin,*) lblock, petzld
       READ (nin,*) t, tout
       READ (nin,*) itol, isplit
       READ (nin,*) y(1:neq)
       SELECT CASE (itol)
       CASE (1)
          READ (nin,*) rtol(1), atol(1)
       CASE (2)
          READ (nin,*) rtol(1), atol(1:neq)
       CASE (3)
          READ (nin,*) rtol(1:neq), atol(1)
       CASE (4)
          READ (nin,*) rtol(1:neq), atol(1:neq)
       END SELECT

       itask = 3
       lderiv(1:2) = .FALSE.
       con(1:6) = zero
       nfails = 0

!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, &
          maxstp,mxhnil,'Average-l2',rwork,ifail)

       ifail = 0
       CALL d02nuf(neq,neq,'Analytical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
          sens,u,eta,lblock,isplit,rwork,ifail)

       irevcm = 0
       WRITE (nout,*) '    X          Y(1)           Y(2)           Y(3)'
       WRITE (nout,99999) t, (y(i),i=1,neq)
       FLUSH (nout)

REVCM: DO
          ifail = -1
          CALL d02nnf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
             ysav,sdysav,wkjac,nwkjac,jacpvt,njcpvt,imon,inln,ires,irevcm, &
             lderiv,itask,itrace,ifail)

          SELECT CASE (irevcm)
          CASE (0)
!            Final exit.
             EXIT REVCM
          CASE (1,3,4,6,7,11)
             IF (irevcm==4 .OR. irevcm==7) THEN
                indr = 50 + neq
             ELSE
                indr = 50 + 2*neq
             END IF
!            Return residual in rwork(indr+1:indr+neq) using y' in ydot.
             rwork(indr+1) = -ydot(1) - ydot(2) - ydot(3)
             rwork(indr+2) = -ydot(2)
             rwork(indr+3) = -ydot(3)
             IF (ires==1) THEN
                rwork(indr+1) = rwork(indr+1) + zero
                rwork(indr+2) = rwork(indr+2) + alpha*y(1) - beta*y(2)*y(3) - &
                   gamma*y(2)*y(2)
                rwork(indr+3) = rwork(indr+3) + gamma*y(2)*y(2)
             END IF
          CASE (2)
!            Return residual in rwork(51+2*neq:) using y' in rwork(51+neq:).
             indd = 50 + neq
             indr = 50 + 2*neq
             rwork(indr+1) = -rwork(indd+1) - rwork(indd+2) - rwork(indd+3)
             rwork(indr+2) = -rwork(indd+2)
             rwork(indr+3) = -rwork(indd+3)
          CASE (5)
!            Return residual in ydot, using y' in rwork(51+2*neq:).
             indd = 50 + 2*neq
             ydot(1) = -rwork(indd+1) - rwork(indd+2) - rwork(indd+3)
             ydot(2) = -rwork(indd+2)
             ydot(3) = -rwork(indd+3)
             ydot(1) = ydot(1) + zero
             ydot(2) = ydot(2) + alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
             ydot(3) = ydot(3) + gamma*y(2)*y(2)
          CASE (8)
!            Return Jacobian in rwork(51+neq:) or rwork(51+2*neq:).

!            Get index J for Jacoban evaluation.
             CALL d02nrf(j,iplace,inform)

             hxd = rwork(16)*rwork(20)
             IF (iplace<2) THEN
!               return Jacobian in rwork(51+2*neq:).
                indr = 50 + 2*neq
             ELSE
!               return Jacobian in rwork(51+neq:).
                indr = 50 + neq
             END IF
!            8 non-zero elements in Jacobian.
             IF (j<2) THEN
                rwork(indr+1) = one - hxd*(zero)
                rwork(indr+2) = zero - hxd*(alpha)
!            rwork(indr+3) = zero - hxd*(zero)
             ELSE IF (j==2) THEN
                rwork(indr+1) = one - hxd*(zero)
                rwork(indr+2) = one - hxd*(-beta*y(3)-gamm2*y(2))
                rwork(indr+3) = zero - hxd*(gamm2*y(2))
             ELSE IF (j>2) THEN
                rwork(indr+1) = one - hxd*(zero)
                rwork(indr+2) = zero - hxd*(-beta*y(2))
                rwork(indr+3) = one - hxd*(zero)
             END IF
          CASE (10)
!            Step failure
             nfails = nfails + 1
          END SELECT
       END DO REVCM

!      Print solution and statistics.
       IF (ifail==0) THEN

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

          ifail = 0
          CALL d02mzf(tout,y,neq,ldysav,neq,ysav,sdysav,rwork,ifail)

          WRITE (nout,99999) tout, (y(i),i=1,neq)
          WRITE (nout,99997) hu, h, tcur
          WRITE (nout,99996) nst, nre, nje
          WRITE (nout,99995) nqu, nq, niter
          WRITE (nout,99994) imxer, nfails
          icall = 0

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

          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,99993) liwreq, liwusd
          WRITE (nout,99992) lrwreq, lrwusd
       ELSE
          WRITE (nout,99998) ifail, t
       END IF


99999  FORMAT (1X,F8.3,3(F13.5,2X))
99998  FORMAT (/1X,'Exit D02NNF with IFAIL = ',I5,'  and T = ',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,' Max err comp = ',I4,'   No. of failed steps = ',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)
    END PROGRAM d02nnfe
