! D02NEF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02nefe_mod ! D02NEF 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 INTEGER, PARAMETER :: ml = 1, mu = 2, neq1 = 3, & neq2 = 1, nin = 5, nout = 6 CONTAINS SUBROUTINE myjac1(neq,ml,mu,t,y,ydot,pd,cj) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: cj, t INTEGER, INTENT (IN) :: ml, mu, neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: pd(2*ml+mu+1,neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq), ydot(neq) ! .. Local Scalars .. INTEGER :: md, ms ! .. Executable Statements .. ! Main diagonal pdfull(i,i), i=1,neq md = mu + ml + 1 pd(md,1) = -alpha - cj pd(md,2) = -beta*y(3) - two*gamma*y(2) - cj pd(md,3) = -cj ! 1 Sub-diagonal pdfull(i+1:i), i=1,neq-1 ms = md + 1 pd(ms,1) = alpha pd(ms,2) = two*gamma*y(2) ! First super-diagonal pdfull(i-1,i), i=2, neq ms = md - 1 pd(ms,2) = beta*y(3) pd(ms,3) = -beta*y(2) ! Second super-diagonal pdfull(i-2,i), i=3, neq ms = md - 2 pd(ms,3) = beta*y(2) RETURN END SUBROUTINE myjac1 SUBROUTINE myjac2(neq,t,y,ydot,pd,cj) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: cj, t INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: pd(neq*neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq), ydot(neq) ! .. Intrinsic Functions .. INTRINSIC exp ! .. Executable Statements .. pd(1) = -two*y(1) + 0.1E0_nag_wp*t*y(1)*exp(y(1)) RETURN END SUBROUTINE myjac2 SUBROUTINE res1(neq,t,y,ydot,r,ires,iuser,ruser) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) = -alpha*y(1) + beta*y(2)*y(3) - ydot(1) r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) - ydot(2) r(3) = gamma*y(2)*y(2) - ydot(3) RETURN END SUBROUTINE res1 SUBROUTINE jac1(neq,t,y,ydot,pd,cj,iuser,ruser) ! .. Use Statements .. USE nag_library, ONLY : d02nez ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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(*) ! .. Local Scalars .. INTEGER :: ijac, ml, mu ! .. Executable Statements .. ml = iuser(1) mu = iuser(2) ijac = iuser(3) IF (ijac==1) THEN CALL myjac1(neq,ml,mu,t,y,ydot,pd,cj) ELSE CALL d02nez(neq,t,y,ydot,pd,cj,iuser,ruser) END IF RETURN END SUBROUTINE jac1 SUBROUTINE res2(neq,t,y,ydot,r,ires,iuser,ruser) ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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(*) ! .. Intrinsic Functions .. INTRINSIC exp ! .. Executable Statements .. r(1) = 4.0_nag_wp - y(1)**2 + t*0.1E0_nag_wp*exp(y(1)) RETURN END SUBROUTINE res2 SUBROUTINE jac2(neq,t,y,ydot,pd,cj,iuser,ruser) ! .. Use Statements .. USE nag_library, ONLY : d02nez ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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(*) ! .. Local Scalars .. INTEGER :: ijac ! .. Executable Statements .. ijac = iuser(1) IF (ijac==1) THEN CALL myjac2(neq,t,y,ydot,pd,cj) ELSE CALL d02nez(neq,t,y,ydot,pd,cj,iuser,ruser) END IF RETURN END SUBROUTINE jac2 END MODULE d02nefe_mod PROGRAM d02nefe ! D02NEF Example Main Program ! .. Use Statements .. USE d02nefe_mod, ONLY : nout ! .. Implicit None Statement .. IMPLICIT NONE ! .. Executable Statements .. WRITE (nout,*) 'D02NEF Example Program Results' CALL ex1 CALL ex2 CONTAINS SUBROUTINE ex1 ! .. Use Statements .. USE nag_library, ONLY : d02mcf, d02mwf, d02nef, d02npf, nag_wp USE d02nefe_mod, ONLY : jac1, ml, mu, neq1, nin, res1 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h0, hmax, t, tout INTEGER :: i, ifail, ijac, itask, itol, & j, lcom, licom, maxord, neq 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(3) ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) 'D02NEF Example 1' ! Skip heading in data file READ (nin,*) READ (nin,*) maxord neq = neq1 lcom = 40 + (maxord+4)*neq + (2*ml+mu+1)*neq + 2*(neq/(ml+mu+1)+1) licom = 50 + neq ALLOCATE (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq), & icom(licom)) READ (nin,*) ijac, itol READ (nin,*) rtol(1:neq) READ (nin,*) atol(1:neq) READ (nin,*) ydot(1:neq) IF (ijac==1) THEN jceval = 'Analytic' ELSE jceval = 'Numeric' END IF ! Set initial values READ (nin,*) y(1:neq) ! Initialize the problem, specifying that the Jacobian is to be ! evaluated analytically using the provided routine jac. READ (nin,*) hmax, h0 READ (nin,*) t, tout ! 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) ! Specify that the Jacobian is banded. ifail = 0 CALL d02npf(neq,ml,mu,icom,licom,ifail) ! Use the iuser array to pass the band dimensions through to jac. ! An alternative would be to hard code values for ml and mu in jac. iuser(1) = ml iuser(2) = mu iuser(3) = ijac WRITE (nout,99999) (i,i=1,neq) WRITE (nout,99998) t, (y(i),i=1,neq) itask = 0 ! Obtain the solution at 5 equally spaced values of T. LOOP: DO j = 1, 5 ifail = -1 CALL d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res1,jac1,icom,com, & lcom,iuser,ruser,ifail) WRITE (nout,99998) t, (y(i),i=1,neq) IF (ifail/=0) THEN WRITE (nout,99997) ifail EXIT LOOP END IF tout = tout + 0.02_nag_wp CALL d02mcf(icom) END DO LOOP WRITE (nout,*) WRITE (nout,99996) itask 99999 FORMAT (/1X,' t ',5X,3(' Y(',I1,') ')) 99998 FORMAT (1X,F8.4,3X,3(F12.6)) 99997 FORMAT (1X,' ** D02NEF returned with IFAIL = ',I5) 99996 FORMAT (1X,'The integrator completed task, ITASK = ',I4) END SUBROUTINE ex1 SUBROUTINE ex2 ! .. Use Statements .. USE nag_library, ONLY : d02mcf, d02mwf, d02nef, nag_wp USE d02nefe_mod, ONLY : jac2, neq2, nin, res2 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h0, hmax, t, tout INTEGER :: i, ifail, ijac, itask, itol, & j, lcom, licom, maxord, neq 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,*) WRITE (nout,*) 'D02NEF Example 2' WRITE (nout,*) READ (nin,*) READ (nin,*) maxord neq = neq2 lcom = 40 + (maxord+4)*neq + neq*neq licom = 50 + neq ALLOCATE (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq), & icom(licom)) READ (nin,*) ijac, itol READ (nin,*) rtol(1:neq) READ (nin,*) atol(1:neq) READ (nin,*) ydot(1:neq) IF (ijac==1) THEN jceval = 'Analytic' ELSE jceval = 'Numeric' END IF ! Initialize the problem, specifying that the Jacobian is to be ! evaluated analytically using the provided routine jac. READ (nin,*) y(1:neq) READ (nin,*) hmax, h0 READ (nin,*) t, tout ifail = 0 CALL d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom, & ifail) ! Use the iuser array to pass whether numerical or analytic Jacobian ! is to be used. iuser(1) = ijac WRITE (nout,99999) (i,i=1,neq) WRITE (nout,99998) t, y(1:neq) itask = 0 ! Obtain the solution at 5 equally spaced values of t. LOOP: DO j = 1, 5 ifail = -1 CALL d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res2,jac2,icom,com, & lcom,iuser,ruser,ifail) WRITE (nout,99998) t, y(1:neq) IF (ifail/=0) THEN WRITE (nout,99997) ifail EXIT LOOP END IF tout = tout + 0.2_nag_wp CALL d02mcf(icom) END DO LOOP WRITE (nout,*) WRITE (nout,99996) itask 99999 FORMAT (/1X,' t y(',I1,')') 99998 FORMAT (1X,F8.4,3X,3(F12.6)) 99997 FORMAT (1X,' ** D02NEF returned with IFAIL = ',I5) 99996 FORMAT (1X,'The integrator completed task, ITASK = ',I4) END SUBROUTINE ex2 END PROGRAM d02nefe