! D02NDF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. Parameters .. Integer, Parameter :: iset = 1, neq = 3 Integer, Parameter :: nia = neq + 1 Integer, Parameter :: nin = 5 Integer, Parameter :: njcpvt = 20*neq + 100 Integer, Parameter :: nout = 6 Integer, Parameter :: nrw = 50 + 4*neq Integer, Parameter :: nwkjac = 4*neq + 100 Integer, Parameter :: sdysav = 6 Integer, Parameter :: ldysav = neq ! .. Local Scalars .. Real (Kind=nag_wp), Save :: xout 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 .. Integer :: i, ifail ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: z(:) ! .. Executable Statements .. inln = 3 If (imon==1) Then Allocate (z(neq)) loop: Do While (t-hlast=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 nag_library, Only: d02ndf, d02ndz, d02nuf, d02nvf, d02nxf, d02nyf, & nag_wp, x04abf Use d02ndfe_mod, Only: fcn, iset, ldysav, monitr, neq, nia, nin, njcpvt, & nout, nrw, nwkjac, sdysav, xout ! .. 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 Structiure 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 Structiure 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) xout = 2.0_nag_wp ! 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