PROGRAM d02nmfe ! D02NMF Example Program Text ! Mark 23 Release. NAG Copyright 2011. ! .. 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 Functions .. 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