! D02NJF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02njfe_mod ! D02NJF 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 Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp Integer, Parameter :: iset = 1, itrace = 0, nelts = 8, & neq = 3 Integer, Parameter :: nia = neq + 1 Integer, Parameter :: nin = 5 Integer, Parameter :: nja = nelts 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 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) = zero r(2) = -ydot(2) r(3) = -ydot(3) If (ires==1) Then r(1) = y(1) + y(2) + y(3) - one + 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 Subroutine jac(neq,t,y,ydot,h,d,j,pdj) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: d, h, t Integer, Intent (In) :: j, neq ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Inout) :: pdj(neq) Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq) ! .. Local Scalars .. Real (Kind=nag_wp) :: hxd ! .. Executable Statements .. ! 8 non-zero elements in total hxd = h*d If (j==1) Then pdj(1) = zero - hxd*(one) pdj(2) = zero - hxd*(alpha) ! note: pdj(3) is zero Else If (j==2) Then pdj(1) = zero - hxd*(one) pdj(2) = one - hxd*(-beta*y(3)-two*gamma*y(2)) pdj(3) = zero - hxd*(two*gamma*y(2)) Else If (j==3) Then pdj(1) = zero - hxd*(one) pdj(2) = zero - hxd*(-beta*y(2)) pdj(3) = one - hxd*(zero) End If Return End Subroutine jac Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, & hmin,hmax,nqu) ! .. 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) ! .. Executable Statements .. inln = 3 If (y(1)<=0.9_nag_wp) imon = -2 Return End Subroutine monitr End Module d02njfe_mod Program d02njfe ! D02NJF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02njf, d02nuf, d02nvf, d02nxf, d02nyf, nag_wp, & x04abf Use d02njfe_mod, Only: iset, itrace, jac, ldysav, monitr, neq, nia, nin, & nja, njcpvt, nout, nrw, nwkjac, resid ! .. 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, liwreq, liwusd, lrwreq, & lrwusd, maxord, maxstp, mxhnil, & nblock, 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(:), & yinit(:), ysav(:,:) Real (Kind=nag_wp) :: con(6) Integer, Allocatable :: ia(:), ja(:), jacpvt(:) Integer :: inform(23) Logical, Allocatable :: algequ(:) Logical :: lderiv(2) ! .. Executable Statements .. Write (nout,*) 'D02NJF 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),wkjac(nwkjac),y(neq), & yinit(neq),ydot(neq),ysav(ldysav,sdysav),ia(nia),ja(nja), & jacpvt(njcpvt),algequ(neq)) Read (nin,*) ia(1:nia) Read (nin,*) ja(1:nja) outchn = nout Call x04abf(iset,outchn) ! Two cases. In both cases: ! integrate to tout by overshooting (itask=1); ! use B.D.F formulae with a Newton method; ! use the Petzold error test is used (differential algebraic system); ! use default values for the array con; ! employ vector relative tolerance and scalar absolute tolerance. ! the Jacobian is supplied by jac; ! the monitr routine is used to force a return when y(1) < 0.9. Read (nin,*) hmin, hmax, h0, tcrit Read (nin,*) eta, sens, u Read (nin,*) lblock Read (nin,*) tinit, tout Read (nin,*) itol, isplt Read (nin,*) yinit(1:neq) Read (nin,*) rtol(1:neq) Read (nin,*) atol(1) con(1:6) = 0.0_nag_wp itask = 1 petzld = .True. cases: Do icase = 1, 2 ! Initialize t = tinit isplit = isplt y(1:neq) = yinit(1:neq) lderiv(1:2) = .False. 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) ! First case. The Jacobian structure is determined internally by ! calls to jac. ifail = 0 Call d02nuf(neq,neq,'Analytical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, & sens,u,eta,lblock,isplit,rwork,ifail) Write (nout,*) ' Analytic Jacobian, structure not supplied' Case (2) ! Second case. The Jacobian structure is supplied. ifail = 0 Call d02nuf(neq,neq,'Full info',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, & sens,u,eta,lblock,isplit,rwork,ifail) Write (nout,*) ' Analytic Jacobian, structure supplied' End Select Write (nout,99988)(i,i=1,neq) Write (nout,99999) t, (y(i),i=1,neq) ! Soft fail and error messages only ifail = 1 Call d02njf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, & resid,ysav,sdysav,jac,wkjac,nwkjac,jacpvt,njcpvt,monitr,lderiv, & itask,itrace,ifail) If (ifail==0 .Or. ifail==12) Then Write (nout,99999) t, (y(i),i=1,neq) 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 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 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 D02NJF with IFAIL = ', ifail, ' and T = ', & t End If End Do cases 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 d02njfe