! D02NGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02ngfe_mod ! D02NGF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. 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 Integer, Parameter :: iset = 1, itrace = 0, neq = 3, & nin = 5, nout = 6 Integer, Parameter :: nrw = 50 + 4*neq Integer, Parameter :: nwkjac = neq*(neq+1) Integer, Parameter :: ldysav = neq 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) = -ydot(1) r(2) = -ydot(2) r(3) = -ydot(3) If (ires==1) Then r(1) = -alpha*y(1) + beta*y(2)*y(3) + 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 End Module d02ngfe_mod Program d02ngfe ! D02NGF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02nby, d02ngf, d02ngz, d02nsf, d02nvf, d02nyf, & d02xjf, nag_wp, x04abf Use d02ngfe_mod, Only: iset, itrace, ldysav, neq, nin, nout, nrw, & nwkjac, resid ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: h, h0, hmax, hmin, hu, t, tcrit, & tcur, tolsf, tout, xout Integer :: i, ifail, imxer, iout, itask, & itol, maxord, maxstp, mxhnil, & niter, nje, nq, nqu, nre, nst, & outchn, sdysav Logical :: petzld ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: atol(:), rtol(:), rwork(:), & sol(:), wkjac(:), y(:), ydot(:), & ysav(:,:) Real (Kind=nag_wp) :: con(6) Integer :: inform(23) Logical, Allocatable :: algequ(:) Logical :: lderiv(2) ! .. Intrinsic Procedures .. Intrinsic :: real ! .. Executable Statements .. Write (nout,*) 'D02NGF 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),sol(neq),wkjac(nwkjac),y(neq), & ydot(neq),ysav(ldysav,sdysav),algequ(neq)) outchn = nout Call x04abf(iset,outchn) ! Integrate to tout by overshooting tout in one step mode (itask=2) ! using B.D.F formulae with a functional iteration method. ! Default values for the array con are used. Employ vector ! tolerances and the Jacobian is evaluated internally, if necessary. ! monitr subroutine replaced by NAG dummy routine D02NBY. ! Interpolation outside D02NGF using D02XJF. itask = 2 Read (nin,*) petzld Read (nin,*) hmin, hmax, h0, tcrit Read (nin,*) t, tout Read (nin,*) y(1:neq) Read (nin,*) lderiv(1:2) Read (nin,*) itol Read (nin,*) rtol(1:neq) Read (nin,*) atol(1:neq) con(1:6) = 0.0_nag_wp ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02nvf(neq,sdysav,maxord,'Functional-iteration',petzld,con,tcrit, & hmin,hmax,h0,maxstp,mxhnil,'Average-L2',rwork,ifail) ! Linear algebra setup required (in case functional iteration ! encounters any difficulty). ifail = 0 Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail) xout = 0.02E0_nag_wp iout = 1 Write (nout,99993)(i,i=1,neq) Write (nout,99999) t, (y(i),i=1,neq) ! Soft fail and error messages only loop1: Do ifail = -1 Call d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, & resid,ysav,sdysav,d02ngz,wkjac,nwkjac,d02nby,lderiv,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) loop2: Do If (tcur-hu=6) Then 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 Exit loop1 End If Else Exit loop2 End If End Do loop2 Else Write (nout,*) Write (nout,99998) 'Exit D02NGF with IFAIL = ', ifail, ' and T = ', & t End If End Do loop1 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) 99993 Format (/1X,' X ',3(' Y(',I1,') ')) End Program d02ngfe