! D02MWF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02mwfe_mod ! D02MWF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: iset = 1, neq = 5, nin = 5, & nout = 6 Integer, Parameter :: licom = 50 + neq Contains Subroutine res(neq,t,y,ydot,r,ires,iuser,ruser) ! .. 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 (Inout) :: ruser(*) Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq) Integer, Intent (Inout) :: iuser(*) ! .. Executable Statements .. r(1) = y(3) - ydot(1) r(2) = y(4) - ydot(2) r(3) = -y(5)*y(1) - ydot(3) r(4) = -y(5)*y(2) - 1.0_nag_wp - ydot(4) r(5) = y(3)**2 + y(4)**2 - y(5) - y(2) Return End Subroutine res Subroutine jac(neq,t,y,ydot,pd,cj,iuser,ruser) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: cj, t Integer, Intent (In) :: neq ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Inout) :: pd(*), ruser(*) Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq) Integer, Intent (Inout) :: iuser(*) ! .. Executable Statements .. pd(1) = -cj pd(3) = -y(5) pd(7) = -cj pd(9) = -y(5) pd(10) = -1.0_nag_wp pd(11) = 1.0_nag_wp pd(13) = -cj pd(15) = 2.0_nag_wp*y(3) pd(17) = 1.0_nag_wp pd(19) = -cj pd(20) = 2.0_nag_wp*y(4) pd(23) = -y(1) pd(24) = -y(2) pd(25) = -1.0_nag_wp Return End Subroutine jac End Module d02mwfe_mod Program d02mwfe ! D02MWF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02mwf, d02nef, nag_wp, x04abf Use d02mwfe_mod, Only: iset, jac, licom, neq, nin, nout, res ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: g1, g2, h0, hmax, t, tout Integer :: i, ifail, ijac, itask, itol, & lcom, maxord, nadv Character (8) :: jceval ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: atol(:), com(:), rtol(:), y(:), & ydot(:) Real (Kind=nag_wp) :: ruser(1) Integer, Allocatable :: icom(:) Integer :: iuser(1) ! .. Executable Statements .. Write (nout,*) 'D02MWF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) ! neq: number of differential-algebraic equations Read (nin,*) maxord lcom = 40 + (maxord+4)*neq + neq**2 Allocate (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom)) nadv = nout Call x04abf(iset,nadv) Read (nin,*) t, tout Read (nin,*) itol, itask Read (nin,*) rtol(1:neq) Read (nin,*) atol(1:neq) ! Set initial values Read (nin,*) y(1:neq) Read (nin,*) hmax, h0 Read (nin,*) ijac If (ijac==1) Then jceval = 'Analytic' Else jceval = 'Numeric' End If ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail) Write (nout,99995)(i,i=1,neq) Write (nout,99998) t, y(1:neq) ydot(1:neq) = 0.0_nag_wp loop: Do Call d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res,jac,icom,com,lcom, & iuser,ruser,ifail) Write (nout,99998) t, y(1:neq) Write (nout,99999) itask If ((itask>=0) .And. (itask<=3)) Then If (t>=tout) Then g1 = y(1)**2 + y(2)**2 - 1.0_nag_wp g2 = y(1)*y(3) + y(2)*y(4) Write (nout,99997) g1 Write (nout,99996) g2 Exit loop End If Else Exit loop End If End Do loop 99999 Format (/' D02NEF returned with ITASK = ',I4/) 99998 Format (1X,F7.4,2X,5(F11.6)) 99997 Format (1X,'The position-level constraint G1 = ',E12.4) 99996 Format (1X,'The velocity-level constraint G2 = ',E12.4) 99995 Format (/1X,' t ',3X,5(' y(',I1,')')) End Program d02mwfe