Program d02nnfe ! D02NNF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: d02mzf, d02nnf, d02nrf, d02nuf, d02nvf, d02nxf, & d02nyf, nag_wp, x04abf ! .. 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 :: zero = 0.0_nag_wp Real (Kind=nag_wp), Parameter :: gamm2 = 2.0_nag_wp*gamma Integer, Parameter :: iset = 1, itrace = 0, nelts = 8, & neq = 3, nia = 1, nin = 5, nja = 1 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 ! .. Local Scalars .. Real (Kind=nag_wp) :: eta, h, h0, hmax, hmin, hu, hxd, & sens, t, tcrit, tcur, tolsf, tout, u Integer :: i, icall, ifail, igrow, imon, imxer, & indd, indr, inln, iplace, ires, & irevcm, isplit, itask, itol, j, & liwreq, liwusd, lrwreq, lrwusd, & maxord, maxstp, mxhnil, nblock, & nfails, 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(:), ysav(:,:) Real (Kind=nag_wp) :: con(6) Integer :: ia(nia), inform(23), ja(nja) Integer, Allocatable :: jacpvt(:) Logical, Allocatable :: algequ(:) Logical :: lderiv(2) ! .. Executable Statements .. Write (nout,*) 'D02NNF Example Program Results' ! Skip heading in data file Read (nin,*) ! neq: number of differential equations Read (nin,*) maxord, maxstp, mxhnil sdysav = maxord + 1 Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), & ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq)) outchn = nout Write (nout,*) Call x04abf(iset,outchn) ! Integrate towards tout stopping at the first mesh point beyond ! tout (itask=3) using the B.D.F. formulae with a Newton method. ! Employ scalar tolerances and the Jacobian is supplied, but its ! structure is evaluated internally by calls to the Jacobian ! forming part of the program (irevcm=8). Default values for the ! array con are used. Also count the number of step failures ! (irevcm=10). The solution is interpolated using D02MZF to give ! the solution at tout. Read (nin,*) hmin, hmax, h0, tcrit Read (nin,*) eta, sens, u Read (nin,*) lblock, petzld Read (nin,*) t, tout Read (nin,*) itol, isplit Read (nin,*) y(1:neq) Select Case (itol) Case (1) Read (nin,*) rtol(1), atol(1) Case (2) Read (nin,*) rtol(1), atol(1:neq) Case (3) Read (nin,*) rtol(1:neq), atol(1) Case (4) Read (nin,*) rtol(1:neq), atol(1:neq) End Select itask = 3 lderiv(1:2) = .False. con(1:6) = zero nfails = 0 ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, & maxstp,mxhnil,'Average-l2',rwork,ifail) ifail = 0 Call d02nuf(neq,neq,'Analytical',nwkjac,ia,nia,ja,nja,jacpvt,njcpvt, & sens,u,eta,lblock,isplit,rwork,ifail) irevcm = 0 Write (nout,*) ' X Y(1) Y(2) Y(3)' Write (nout,99999) t, (y(i),i=1,neq) Flush (nout) revcm: Do ifail = -1 Call d02nnf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,ysav, & sdysav,wkjac,nwkjac,jacpvt,njcpvt,imon,inln,ires,irevcm,lderiv, & itask,itrace,ifail) Select Case (irevcm) Case (0) ! Final exit. Exit revcm Case (1,3,4,6,7,11) If (irevcm==4 .Or. irevcm==7) Then indr = 50 + neq Else indr = 50 + 2*neq End If ! Return residual in rwork(indr+1:indr+neq) using y' in ydot. rwork(indr+1) = -ydot(1) - ydot(2) - ydot(3) rwork(indr+2) = -ydot(2) rwork(indr+3) = -ydot(3) If (ires==1) Then rwork(indr+1) = rwork(indr+1) + zero rwork(indr+2) = rwork(indr+2) + alpha*y(1) - beta*y(2)*y(3) - & gamma*y(2)*y(2) rwork(indr+3) = rwork(indr+3) + gamma*y(2)*y(2) End If Case (2) ! Return residual in rwork(51+2*neq:) using y' in rwork(51+neq:). indd = 50 + neq indr = 50 + 2*neq rwork(indr+1) = -rwork(indd+1) - rwork(indd+2) - rwork(indd+3) rwork(indr+2) = -rwork(indd+2) rwork(indr+3) = -rwork(indd+3) Case (5) ! Return residual in ydot, using y' in rwork(51+2*neq:). indd = 50 + 2*neq ydot(1) = -rwork(indd+1) - rwork(indd+2) - rwork(indd+3) ydot(2) = -rwork(indd+2) ydot(3) = -rwork(indd+3) ydot(1) = ydot(1) + zero ydot(2) = ydot(2) + alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) ydot(3) = ydot(3) + gamma*y(2)*y(2) Case (8) ! Return Jacobian in rwork(51+neq:) or rwork(51+2*neq:). ! Get index J for Jacoban evaluation. Call d02nrf(j,iplace,inform) hxd = rwork(16)*rwork(20) If (iplace<2) Then ! return Jacobian in rwork(51+2*neq:). indr = 50 + 2*neq Else ! return Jacobian in rwork(51+neq:). indr = 50 + neq End If ! 8 non-zero elements in Jacobian. If (j<2) Then rwork(indr+1) = one - hxd*(zero) rwork(indr+2) = zero - hxd*(alpha) ! rwork(indr+3) = zero - hxd*(zero) Else If (j==2) Then rwork(indr+1) = one - hxd*(zero) rwork(indr+2) = one - hxd*(-beta*y(3)-gamm2*y(2)) rwork(indr+3) = zero - hxd*(gamm2*y(2)) Else If (j>2) Then rwork(indr+1) = one - hxd*(zero) rwork(indr+2) = zero - hxd*(-beta*y(2)) rwork(indr+3) = one - hxd*(zero) End If Case (10) ! Step failure nfails = nfails + 1 End Select End Do revcm ! Print solution and statistics. If (ifail==0) Then Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, & imxer,algequ,inform,ifail) ifail = 0 Call d02mzf(tout,y,neq,ldysav,neq,ysav,sdysav,rwork,ifail) Write (nout,99999) tout, (y(i),i=1,neq) Write (nout,99997) hu, h, tcur Write (nout,99996) nst, nre, nje Write (nout,99995) nqu, nq, niter Write (nout,99994) imxer, nfails icall = 0 Call d02nxf(icall,liwreq,liwusd,lrwreq,lrwusd,nlu,nnz,ngp,isplit, & igrow,lblock,nblock,inform) 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,99993) liwreq, liwusd Write (nout,99992) lrwreq, lrwusd Else Write (nout,99998) ifail, t End If 99999 Format (1X,F8.3,3(F13.5,2X)) 99998 Format (/1X,'Exit D02NNF with IFAIL = ',I5,' and T = ',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,' Max err comp = ',I4,' No. of failed steps = ',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) End Program d02nnfe