Program d02nnfe
! D02NNF Example Program Text
! Mark 26.2 Release. NAG Copyright 2017.
! .. 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