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

NAG FL Interface Introduction
Example description
    Program d02nmfe

!     D02NMF Example Program Text

!     Mark 30.1 Release. NAG Copyright 2024.

!     .. Use Statements ..
      Use nag_library, Only: d02nmf, d02nsf, d02nvf, d02nyf, d02xkf, nag_wp,   &
                             x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: tstep = 2.0E0_nag_wp
      Integer, Parameter               :: iset = 1, neq = 3, nin = 5,          &
                                          njcpvt = 1, nout = 6
      Integer, Parameter               :: nrw = 50 + 4*neq
      Integer, Parameter               :: nwkjac = neq*(neq+1)
      Integer, Parameter               :: sdysav = 6
      Integer, Parameter               :: ldysav = neq
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h, h0, hlast, hmax, hmin, hnext, hu, &
                                          t, tc, tcrit, tcur, tolsf, tout,     &
                                          xout
      Integer                          :: i, ifail, iflag, imon, imxer, inln,  &
                                          ires, irevcm, itask, itol, itrace,   &
                                          lacor1, lacor2, lacor3, lacorb,      &
                                          lsavr1, lsavr2, lsavr3, lsavrb,      &
                                          maxord, maxstp, mxhnil, niter, nje,  &
                                          nq, nqu, nre, nst, outchn
      Logical                          :: petzld
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: atol(:), rtol(:), rwork(:),          &
                                          wkjac(:), y(:), ydot(:), ysav(:,:)
      Real (Kind=nag_wp)               :: con(6)
      Integer                          :: inform(23)
      Integer, Allocatable             :: jacpvt(:)
      Logical, Allocatable             :: algequ(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: int
!     .. Executable Statements ..
      Write (nout,*) 'D02NMF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     neq: number of differential equations
      Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), &
        ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq))

!     Integrate to tout by overshooting tout (itask=1) using B.D.F.
!     formulae with a Newton method. Default values for the array con
!     are used. Employ scalar tolerances and the Jacobian is evaluated
!     internally. On the reverse communication call equivalent to the
!     monitr call in forward communication routines carry out
!     interpolation using D02XKF.

      Read (nin,*) maxord, maxstp, mxhnil
      Read (nin,*) hmin, hmax, h0, tcrit
      Read (nin,*) petzld
      Read (nin,*) t, tout
      Read (nin,*) itol
      Read (nin,*) y(1:neq)
      Read (nin,*) rtol(1), atol(1)

      outchn = nout
      Call x04abf(iset,outchn)
      itask = 1
      xout = tstep
      con(1:6) = 0.0E0_nag_wp

      ifail = 0
      Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0,    &
        maxstp,mxhnil,'Average-L2',rwork,ifail)

      ifail = 0
      Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)

      lacorb = 50 + neq
      lacor1 = lacorb + 1
      lacor2 = lacorb + 2
      lacor3 = lacorb + 3
      lsavrb = lacorb + neq
      lsavr1 = lsavrb + 1
      lsavr2 = lsavrb + 2
      lsavr3 = lsavrb + 3
      Write (nout,*) '    X          Y(1)           Y(2)           Y(3)'
      Write (nout,99999) t, (y(i),i=1,neq)

!     Soft fail and error messages only
      irevcm = 0
      itrace = 0

steps: Do
        ifail = 1
        Call d02nmf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,ysav, &
          sdysav,wkjac,nwkjac,jacpvt,njcpvt,imon,inln,ires,irevcm,itask,       &
          itrace,ifail)

        Select Case (irevcm)
        Case (0)
          Exit steps
        Case (1,3)
!         Equivalent to fcn evaluation in forward communication
!         routines
          rwork(lsavr1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
          rwork(lsavr2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) -        &
            3.0E7_nag_wp*y(2)*y(2)
          rwork(lsavr3) = 3.0E7_nag_wp*y(2)*y(2)
        Case (4)
!         Equivalent to fcn evaluation in forward communication
!         routines
          rwork(lacor1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
          rwork(lacor2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) -        &
            3.0E7_nag_wp*y(2)*y(2)
          rwork(lacor3) = 3.0E7_nag_wp*y(2)*y(2)
        Case (5)
!         Equivalent to fcn evaluation in forward communication
!         routines
          ydot(1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
          ydot(2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) -              &
            3.0E7_nag_wp*y(2)*y(2)
          ydot(3) = 3.0E7_nag_wp*y(2)*y(2)
        Case (9)
!         Equivalent to monitr call in forward communication routines
          If (imon==1) Then
            tc = rwork(19)
            hlast = rwork(15)
            hnext = rwork(16)
            nqu = int(rwork(10))
interp:     Do
              If (tc-hlast<xout .And. xout<=tc) Then
                iflag = 1

                Call d02xkf(xout,rwork(lsavr1),neq,ysav,ldysav,sdysav,         &
                  rwork(lacor1),neq,tc,nqu,hlast,hnext,iflag)

                If (iflag/=0) Then
                  imon = -2
                  Exit interp
                Else
                  Write (nout,99999) xout, (rwork(lsavrb+i),i=1,neq)
                  xout = xout + tstep
                  If (xout>tout) Then
                    Exit interp
                  End If
                End If
              Else
                Exit interp
              End If
            End Do interp
          End If
        Case (2,6:8)
          Write (nout,*)
          Write (nout,99994) 'Illegal value of IREVCM = ', irevcm
          Exit steps
        End Select
      End Do steps
      If (ifail==0) Then

        iflag = 0
        Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter,    &
          imxer,algequ,inform,iflag)

        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,*)
      Else
        Write (nout,*)
        Write (nout,99998) 'Exit D02NMF with IFAIL = ', ifail, '  and T = ', t
      End If

99999 Format (1X,F8.3,3(F13.5,2X))
99998 Format (1X,A,I2,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)
    End Program d02nmfe