PROGRAM drddasl2
c>> 2009-10-27 DRDDASL2 Krogh Declared DUM(4,*) for NAG compiler.
c>> 2009-10-19 DRDDASL2 Krogh Changed def. of d(4,4) so 0 doesn't abort.
c>> 2008-10-26 DRDDASL2 Krogh Moved Format statements up for C conv.
c>> 2008-10-24 DRDDASL2 Krogh Removed in INCLUDE statement & cDEC$...
C>> 2008-09-04 DRDDASL2 Hanson added starting computation of y'
C>> 2008-08-26 DRDDASL2 Hanson added row dimensions to evaluators
C>> 2001-10-11 DRDDASL2 R. J. Hanson Example 2 Code, with Download
C Solve Enright and Pryce stiff test problem E5.
C The equation is presented as y'= F(t,y), y_0 given.
C This is solved by defining the residual function
C f(t,y,y')=F(t,y)-y'. Integration is done twice.
C The first time analytic partials are provided.
C The second time, difference quotients are used for
C partials, with the integrator requesting only values
C of f(t,y,y'). The two solutions are equivalent but
C not exactly equal.
C--D replaces "?": DR?DASL2, ?DASLS, ?DASLX, ?DASSF2, ?DASSF3, ?COPY
C-- & ?ROTG
EXTERNAL ddassf2,ddassf3
INTEGER ndig,nepoch,neq,maxcon,ldc,ltd,lrw,liw
C Set number of equations:
PARAMETER (neq=4)
C Set number of constraints.
PARAMETER (maxcon=0)
C Work space sizes:
PARAMETER (liw=30+neq)
PARAMETER (lrw=45+(5+2*maxcon+4)*neq+neq**2)
PARAMETER (ldc=neq,ltd=neq)
DOUBLE PRECISION tol
C++S Default NDIG = 3
C++ Default NDIG = 8
C++ Substitute for NDIG below
PARAMETER (ndig=8,nepoch=10)
PARAMETER (tol=10.d0**(-ndig))
INTEGER i,j,info(16),idid,iwork(liw)
DOUBLE PRECISION t,y(neq),yprime(neq),rtol(neq),atol(neq),
& rwork(lrw),tout(nepoch),y1(neq,nepoch),y2(neq,nepoch)
DOUBLE PRECISION c(ldc,ltd),ftol,rnktol
LOGICAL match
INTEGER kf2,ks2,kf3,ks3
COMMON /counts/ kf2,ks2,kf3,ks3
90 FORMAT (/20x,'Example Results for a Stiff ODE Problem, E5')
100 FORMAT (6x,'Integration Values do not Match')
110 FORMAT (/15x,'No. of Residual Evaluations',6x,'No. of User Solves'
& /' Partials Provided-',8x,i4,24x,i4)
120 FORMAT (' Divided Differences-',6x,i4,24x,i4)
130 FORMAT (/10x,'T',7x,'RErr in y_1',3x,'RErr in y_2',3x,'RErr in y_3
& ',3x,'RErr in y_4')
140 FORMAT (1p,5d14.4)
C Tolerances:
DO 10 i=1,neq
atol(i)=tol
rtol(i)=tol
10 CONTINUE
C Setup options:
DO 20 i=1,16
info(i)=0
20 CONTINUE
C Compute, factor and solve the partial derivative matrix in routine
C DDASSF2. This illustrates how to completely control the linear
C algebra. Using analytic partial derivatives.
info(5)=5
C Compute the initial value of YPRIME(*):
t=0
ftol=tol
rnktol=tol
C Give starting values to y and y' before Newton method:
do 25 i = 1, neq
y(i) = 0.D0
yprime(i) = 0.D0
25 continue
y(1) = 1.76D-3
CALL ddasls (ddassf2, neq, t, y, yprime, info, ftol, rnktol, c,
& ldc, ltd, idid, rwork, lrw, iwork, liw)
C Set output points for integration:
tout(1)=0.d0
tout(2)=1.d-3
DO 30 i=3,nepoch
tout(i)=10*tout(i-1)
30 CONTINUE
DO 40 i=2,nepoch
t=tout(i-1)
CALL ddaslx (ddassf2, neq, t, y, yprime, tout(i), info, rtol,
& atol, idid, rwork, lrw, iwork, liw)
CALL dcopy (neq, y, 1, y1(1,i), 1)
40 CONTINUE
C Compute the initial value of YPRIME(*):
t=0
ftol=tol
rnktol=tol
C Give starting values to y and y' before Newton method:
do 45 i = 1, neq
y(i) = 0.D0
yprime(i) = 0.D0
45 continue
y(1) = 1.76D-3
CALL ddasls (ddassf3, neq, t, y, yprime, info, ftol, rnktol, c,
& ldc, ltd, idid, rwork, lrw, iwork, liw)
C Restart the integration, but do not use analytic derivatives.
info(1)=0
C Use divided differences instead of user-provided partials.
C And the default linear solver provided.
info(5)=1
DO 50 i=2,nepoch
t=tout(i-1)
CALL ddaslx (ddassf3, neq, t, y, yprime, tout(i), info, rtol,
& atol, idid, rwork, lrw, iwork, liw)
CALL dcopy (neq, y, 1, y2(1,i), 1)
50 CONTINUE
C See if the solutions "match" the requested tolerances.
match=.true.
DO 70 j=2,nepoch
DO 60 i=1,neq
y2(i,j)=(y2(i,j)-y1(i,j))/(abs(y1(i,j))*rtol(i)+atol(i))
match=match.and.(abs(y2(i,j)).le.100.d0)
60 CONTINUE
70 CONTINUE
C Output the relative errors and a summary.
WRITE (*,90)
WRITE (*,130)
DO 80 j=2,nepoch
WRITE (*,140) tout(j),(y2(i,j),i=1,neq)
80 CONTINUE
IF (match) THEN
WRITE (*,110) kf2,ks2
WRITE (*,120) kf3,ks3
ELSE
WRITE (*,100)
END IF
END
SUBROUTINE ddassf2 (t, y, yprime, delta, dum, ldd, cj, ires,
& rwork, iwork)
C
C Routine for the Enright and Pryce problem E5.
C DUM(*,*) is not used in this example.
DOUBLE PRECISION t,y(*),yprime(*),delta(*),dum(ldd,*),cj,rwork(*),
& dc(3),ds(3),u,b0,b1,b2,b3,b4
DOUBLE PRECISION d(4,4)
INTEGER i,ires,iwork(*),ldd
SAVE b0,b1,b2,b3,b4,d,dc,ds
INTEGER kf2,ks2,kf3,ks3
COMMON /counts/ kf2,ks2,kf3,ks3
C This is the setup call.
IF (ires.eq.0) THEN
b0=1.76d-03
b1=7.89d-10
b2=1.1d+07
b3=1.13d+09
b4=1.13d+03
C Set initial conditions.
y(1)=0.d0
CALL dcopy (4, y, 0, y, 1)
y(1)=b0
C Count evaluations in KF2 and KS2.
kf2=0
ks2=0
RETURN
END IF
C The sytem residual value:
IF (ires.eq.1) THEN
delta(1)=-b1*y(1)-b2*y(1)*y(3)-yprime(1)
delta(2)=b1*y(1)-b3*y(2)*y(3)-yprime(2)
delta(3)=b1*y(1)-b2*y(1)*y(3)+b4*y(4)-b3*y(2)*y(3)-yprime(3)
delta(4)=b2*y(1)*y(3)-b4*y(4)-yprime(4)
kf2=kf2+1
RETURN
END IF
C The partial derivative matrix:
IF (ires.eq.2) THEN
d(4,1)=b2*y(3)
d(1,1)=-b1-d(4,1)-cj
d(2,1)=b1
d(3,1)=b1-d(4,1)
d(1,2)=0.d0
d(3,2)=-b3*y(3)
d(2,2)=d(3,2)-cj
d(4,2)=0.d0
d(1,3)=-b2*y(1)
d(2,3)=-b3*y(2)
d(3,3)=d(1,3)+d(2,3)-cj
d(4,3)=-d(1,3)
d(1,4)=0.d0
d(2,4)=0.d0
d(3,4)=b4
d(4,4)=-b4-cj
C This matrix is right factored to lower triangular for with three
C plane rotations. Use planes 1 and 3:
C Tell integrator system is non-singular.
ires=0
CALL drotg (d(1,1), d(1,3), dc(1), ds(1))
u=dc(1)*d(2,1)+ds(1)*d(2,3)
d(2,3)=-ds(1)*d(2,1)+dc(1)*d(2,3)
d(2,1)=u
u=dc(1)*d(3,1)+ds(1)*d(3,3)
d(3,3)=-ds(1)*d(3,1)+dc(1)*d(3,3)
d(3,1)=u
u=dc(1)*d(4,1)+ds(1)*d(4,3)
d(4,3)=-ds(1)*d(4,1)+dc(1)*d(4,3)
d(4,1)=u
d(1,1)=1.d0/d(1,1)
C Eliminate in planes 2 and 3.
CALL drotg (d(2,2), d(2,3), dc(2), ds(2))
u=dc(2)*d(3,2)+ds(2)*d(3,3)
d(3,3)=-ds(2)*d(3,2)+dc(2)*d(3,3)
d(3,2)=u
u=dc(2)*d(4,2)+ds(2)*d(4,3)
d(4,3)=-ds(2)*d(4,2)+dc(2)*d(4,3)
d(4,2)=u
d(2,2)=1.d0/d(2,2)
C Eliminate in planes 3 and 4.
CALL drotg (d(3,3), d(3,4), dc(3), ds(3))
u=dc(3)*d(4,3)+ds(3)*d(4,4)
d(4,4)=(-ds(3)*d(4,3)+dc(3)*d(4,4))
if (d(4,4) .eq. 0.D0) then
ires = 1
return
end if
d(4,4) = 1.d0 / d(4,4)
d(4,3)=u
d(3,3)=1.d0/d(3,3)
RETURN
END IF
C Solve the corrector equation.
IF (ires.eq.4) THEN
C Count number of solves.
ks2=ks2+1
C Forward substitute:
delta(1)=delta(1)*d(1,1)
delta(4)=delta(4)-delta(1)*d(4,1)
delta(3)=delta(3)-delta(1)*d(3,1)
delta(2)=delta(2)-delta(1)*d(2,1)
delta(2)=delta(2)*d(2,2)
delta(4)=delta(4)-delta(2)*d(4,2)
delta(3)=delta(3)-delta(2)*d(3,2)
delta(3)=delta(3)*d(3,3)
delta(4)=delta(4)-delta(3)*d(4,3)
delta(4)=delta(4)*d(4,4)
C Apply three plane rotations,(3,4),(2,3),(1,3).
u=dc(3)*delta(3)-ds(3)*delta(4)
delta(4)=ds(3)*delta(3)+dc(3)*delta(4)
delta(3)=u
u=dc(2)*delta(2)-ds(2)*delta(3)
delta(3)=ds(2)*delta(2)+dc(2)*delta(3)
delta(2)=u
u=dc(1)*delta(1)-ds(1)*delta(3)
delta(3)=ds(1)*delta(1)+dc(1)*delta(3)
delta(1)=u
RETURN
END IF
C This is df/dy' for the starting procedure. This
C problem is a linear ODE.
IF (ires.eq.7) THEN
DO 10 i=1,4
d(i,i)=-1.d0
10 CONTINUE
RETURN
END IF
END
SUBROUTINE ddassf3 (t, y, yprime, delta, d, ldd, cj, ires, rwork,
& iwork)
C
C Routine for the Enright and Pryce problem E5.
DOUBLE PRECISION t,y(*),yprime(*),delta(*),d(ldd,4),cj,rwork(*),
& dc(3),ds(3),b1,b2,b3,b4
INTEGER ldd, ires,iwork(*)
PARAMETER(b1=7.89d-10, b2=1.1d+07, b3=1.13d+09, b4=1.13d+03)
INTEGER i,kf2,ks2,kf3,ks3
COMMON /counts/ kf2,ks2,kf3,ks3
C This is the setup call.
IF (ires.eq.0) THEN
C Count evaluations in KF3 and KS3.
kf3=0
ks3=0
END IF
C The sytem residual value:
IF (ires.eq.1) THEN
delta(1)=-b1*y(1)-b2*y(1)*y(3)-yprime(1)
delta(2)=b1*y(1)-b3*y(2)*y(3)-yprime(2)
delta(3)=b1*y(1)-b2*y(1)*y(3)+b4*y(4)-b3*y(2)*y(3)-yprime(3)
delta(4)=b2*y(1)*y(3)-b4*y(4)-yprime(4)
kf3=kf3+1
END IF
C The partial derivative matrix:
IF (ires.eq.2) THEN
d(4,1)=b2*y(3)
d(1,1)=-b1-d(4,1)-cj
d(2,1)=b1
d(3,1)=b1-d(4,1)
d(2,2)=-b3*y(3)-cj
d(3,2)=d(2,2)
d(1,3)=-b2*y(1)
d(2,3)=-b3*y(2)
d(3,3)=d(1,3)+d(2,3)-cj
d(4,3)=-d(1,3)
d(3,4)=b4
d(4,4)=-b4-cj
END IF
IF (ires.eq.3) THEN
WRITE (*,*) 'Should be no call to DDASSF3 with IRES = 3.'
END IF
IF (ires.eq.4) THEN
WRITE (*,*) 'Should be no call to DDASSF3 with IRES = 4.'
END IF
C This is df/dy' for the starting procedure. This
C problem is a linear ODE of index 0. So only IRES=1,6 occur.
IF (ires.eq.7) THEN
DO 10 i=1,4
d(i,i)=-1.d0
10 CONTINUE
RETURN
END IF
END