! D02NDF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d02ndfe_mod
! D02NDF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: fcn, monitr
! .. Parameters ..
Integer, Parameter, Public :: iset = 1, neq = 3
Integer, Parameter, Public :: nia = neq + 1
Integer, Parameter, Public :: nin = 5
Integer, Parameter, Public :: njcpvt = 20*neq + 100
Integer, Parameter, Public :: nout = 6
Integer, Parameter, Public :: nrw = 50 + 4*neq
Integer, Parameter, Public :: nwkjac = 4*neq + 100
Integer, Parameter, Public :: sdysav = 6
Integer, Parameter, Public :: ldysav = neq
Contains
Subroutine fcn(neq,t,y,f,ires)
! .. 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
! .. 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) :: f(neq)
Real (Kind=nag_wp), Intent (In) :: y(neq)
! .. Executable Statements ..
f(1) = -alpha*y(1) + beta*y(2)*y(3)
f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
f(3) = gamma*y(2)*y(2)
Return
End Subroutine fcn
Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
hmin,hmax,nqu)
! .. Use Statements ..
Use nag_library, Only: d02xkf
! .. 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)
! .. Local Scalars ..
Real (Kind=nag_wp) :: xout
Integer :: i, ifail
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: z(:)
! .. Executable Statements ..
inln = 3
If (imon==1) Then
Allocate (z(neq))
! Interpolate at multiples of 2.0 between t-hlast and t
xout = 2.0E0_nag_wp
Do While (xout<=t-hlast)
xout = xout + 2.0E0_nag_wp
End Do
loop: Do While (t-hlast<xout .And. xout<=t)
! C1 interpolation
ifail = 1
Call d02xkf(xout,z,neq,ysav,ldysav,sdysav,acor(1,2),neq,t,nqu, &
hlast,hnext,ifail)
If (ifail/=0) Then
imon = -2
Exit loop
End If
Write (nout,99999) xout, (z(i),i=1,neq)
xout = xout + 2.0E0_nag_wp
If (xout>=10.0E0_nag_wp) Then
Exit loop
End If
End Do loop
End If
Deallocate (z)
Return
99999 Format (1X,F8.3,3(F13.5,2X))
End Subroutine monitr
End Module d02ndfe_mod
Program d02ndfe
! D02NDF Example Main Program
! .. Use Statements ..
Use d02ndfe_mod, Only: fcn, iset, ldysav, monitr, neq, nia, nin, njcpvt, &
nout, nrw, nwkjac, sdysav
Use nag_library, Only: d02ndf, d02ndz, d02nuf, d02nvf, d02nxf, d02nyf, &
nag_wp, x04abf
! .. 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, &
itrace, liwreq, liwusd, lrwreq, &
lrwusd, maxord, maxstp, mxhnil, &
nblock, ngp, niter, nja, nje, nlu, &
nnz, nq, nqu, nre, nst, outchn
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(:)
! .. Executable Statements ..
Write (nout,*) 'D02NDF Example Program Results'
! Skip heading in data file
Read (nin,*)
! Allocations based on number of differential equations (neq)
Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
yinit(neq),ydot(neq),ysav(ldysav,sdysav),ia(nia),jacpvt(njcpvt), &
algequ(neq))
Read (nin,*) maxord, maxstp, mxhnil
Read (nin,*) ia(1:nia)
nja = ia(nia) - 1
Allocate (ja(nja))
Read (nin,*) ja(1:nja)
! Read algorithmic parameters
Read (nin,*) hmin, hmax, h0, tcrit
Read (nin,*) eta, sens, u
Read (nin,*) lblock, petzld
Read (nin,*) tinit, tout
Read (nin,*) itol, isplt
Read (nin,*) yinit(1:neq)
Read (nin,*) rtol(1), atol(1)
outchn = nout
Call x04abf(iset,outchn)
Do icase = 1, 2
t = tinit
isplit = isplt
y(1:neq) = yinit(1:neq)
! In both cases: integrate to tout by overshooting (itask=1) using BDF
! with Newton; use default con values, scalar tolerances and numerical
! Jacobian; interpolate using MONITR and D02XKF.
con(1:6) = 0.0_nag_wp
itask = 1
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)
! No Jacobian Structure Supplied.
ifail = 0
Call d02nuf(neq,neq,'Numerical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
sens,u,eta,lblock,isplit,rwork,ifail)
Write (nout,*) ' Numerical Jacobian, structure not supplied'
Case (2)
! Jacobian Structure Supplied.
ifail = 0
Call d02nuf(neq,neq,'Structural',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
sens,u,eta,lblock,isplit,rwork,ifail)
Write (nout,*) ' Numerical 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
itrace = 0
ifail = -1
Call d02ndf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, &
ysav,sdysav,d02ndz,wkjac,nwkjac,jacpvt,njcpvt,monitr,itask,itrace, &
ifail)
If (ifail==0) Then
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
Write (nout,*)
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
Write (nout,*)
Write (nout,99998) 'Exit D02NDF with IFAIL = ', ifail, ' and T = ', &
t
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 D02NDF with IFAIL = ', ifail, ' and T = ', &
t
End If
End Do
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 d02ndfe