Program d02nmfe ! D02NMF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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-hlasttout) Exit interp 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