NAG Library Manual, Mark 30.2
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program d02nnfe

!     D02NNF Example Program Text

!     Mark 30.2 Release. NAG Copyright 2024.

!     .. 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 nonzero 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