! D02NJF Example Program Text
! Mark 27.3 Release. NAG Copyright 2021.
Module d02njfe_mod
! D02NJF 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 :: jac, monitr, resid
! .. 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 :: two = 2.0_nag_wp
Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp
Integer, Parameter, Public :: iset = 1, itrace = 0, neq = 3
Integer, Parameter, Public :: nia = neq + 1
Integer, Parameter, Public :: nin = 5, nout = 6
Integer, Parameter, Public :: nrw = 50 + 4*neq
Integer, Parameter, Public :: ldysav = neq
Integer, Parameter :: nelts = 8
Integer, Parameter, Public :: nja = nelts
Integer, Parameter, Public :: njcpvt = 20*neq + 12*nelts
Integer, Parameter, Public :: nwkjac = 4*neq + 12*nelts
Contains
Subroutine resid(neq,t,y,ydot,r,ires)
! .. 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) :: r(neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Executable Statements ..
r(1) = zero
r(2) = -ydot(2)
r(3) = -ydot(3)
If (ires==1) Then
r(1) = y(1) + y(2) + y(3) - one + r(1)
r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2)
r(3) = gamma*y(2)*y(2) + r(3)
End If
Return
End Subroutine resid
Subroutine jac(neq,t,y,ydot,h,d,j,pdj)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: d, h, t
Integer, Intent (In) :: j, neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: pdj(neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: hxd
! .. Executable Statements ..
! 8 nonzero elements in total
hxd = h*d
If (j==1) Then
pdj(1) = zero - hxd*(one)
pdj(2) = zero - hxd*(alpha)
! note: pdj(3) is zero
Else If (j==2) Then
pdj(1) = zero - hxd*(one)
pdj(2) = one - hxd*(-beta*y(3)-two*gamma*y(2))
pdj(3) = zero - hxd*(two*gamma*y(2))
Else If (j==3) Then
pdj(1) = zero - hxd*(one)
pdj(2) = zero - hxd*(-beta*y(2))
pdj(3) = one - hxd*(zero)
End If
Return
End Subroutine jac
Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
hmin,hmax,nqu)
! .. 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)
! .. Executable Statements ..
inln = 3
If (y(1)<=0.9_nag_wp) Then
imon = -2
End If
Return
End Subroutine monitr
End Module d02njfe_mod
Program d02njfe
! D02NJF Example Main Program
! .. Use Statements ..
Use d02njfe_mod, Only: iset, itrace, jac, ldysav, monitr, neq, nia, nin, &
nja, njcpvt, nout, nrw, nwkjac, resid
Use nag_library, Only: d02njf, 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, &
liwreq, liwusd, lrwreq, lrwusd, &
maxord, maxstp, mxhnil, nblock, 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(:), yinit(:), &
ysav(:,:)
Real (Kind=nag_wp) :: con(6)
Integer, Allocatable :: ia(:), ja(:), jacpvt(:)
Integer :: inform(23)
Logical, Allocatable :: algequ(:)
Logical :: lderiv(2)
! .. Executable Statements ..
Write (nout,*) 'D02NJF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) maxord, maxstp, mxhnil
sdysav = maxord + 1
Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
yinit(neq),ydot(neq),ysav(ldysav,sdysav),ia(nia),ja(nja), &
jacpvt(njcpvt),algequ(neq))
Read (nin,*) ia(1:nia)
Read (nin,*) ja(1:nja)
outchn = nout
Call x04abf(iset,outchn)
! Two cases. In both cases:
! integrate to tout by overshooting (itask=1);
! use B.D.F formulae with a Newton method;
! use the Petzold error test (differential algebraic system);
! use default values for the array con;
! employ vector relative tolerance and scalar absolute tolerance.
! the Jacobian is supplied by jac;
! the monitr routine is used to force a return when y(1) < 0.9.
Read (nin,*) hmin, hmax, h0, tcrit
Read (nin,*) eta, sens, u
Read (nin,*) lblock
Read (nin,*) tinit, tout
Read (nin,*) itol, isplt
Read (nin,*) yinit(1:neq)
Read (nin,*) rtol(1:neq)
Read (nin,*) atol(1)
con(1:6) = 0.0_nag_wp
itask = 1
petzld = .True.
cases: Do icase = 1, 2
! Initialize
t = tinit
isplit = isplt
y(1:neq) = yinit(1:neq)
lderiv(1:2) = .False.
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)
! First case. The Jacobian structure is determined internally by
! calls to jac.
ifail = 0
Call d02nuf(neq,neq,'Analytical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
sens,u,eta,lblock,isplit,rwork,ifail)
Write (nout,*) ' Analytic Jacobian, structure not supplied'
Case (2)
! Second case. The Jacobian structure is supplied.
ifail = 0
Call d02nuf(neq,neq,'Full info',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, &
sens,u,eta,lblock,isplit,rwork,ifail)
Write (nout,*) ' Analytic 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
ifail = 1
Call d02njf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
resid,ysav,sdysav,jac,wkjac,nwkjac,jacpvt,njcpvt,monitr,lderiv, &
itask,itrace,ifail)
If (ifail==0 .Or. ifail==12) Then
Write (nout,99999) t, (y(i),i=1,neq)
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
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
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 D02NJF with IFAIL = ', ifail, ' and T = ', &
t
End If
End Do cases
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 d02njfe