Program d02nmfe
! D02NMF Example Program Text
! Mark 29.2 Release. NAG Copyright 2023.
! .. 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