Fix f2c compilation
[scilab.git] / scilab / modules / differential_equations / tests / unit_tests / daskr.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING,
7 // which you should have received as part of this distribution.
8 // The terms are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 // =============================================================================
11
12 // Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
13
14 //C-----------------------------------------------------------------------
15 //C First problem.
16 //C The initial value problem is..
17 //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  1 .LE. T .LE. 6
18 //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
19 //C The two root functions are..
20 //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
21 //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
22 //C-----------------------------------------------------------------------
23 y0=1;t=2:6;t0=1;y0d=3;
24 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
25 atol=1.d-6;rtol=0;ng=2;
26 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
27 if abs(nn(1)-2.47)>0.001 then pause,end
28 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
29 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info,'psol1','pjac1');
30 if abs(nn(1)-2.5)>0.001 then pause,end
31 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
32 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
33 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
34 if abs(nn(1)-2.53)>0.001 then pause,end
35
36 // Same problem, but using macro for the derivative evaluation function 'res1'
37 deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y')
38 deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
39
40 y0=1;t=2:6;t0=1;y0d=3;
41 atol=1.d-6;rtol=0;ng=2;
42 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
43 if abs(nn(1)-2.47)>0.001 then pause,end
44 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
45 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
46 if abs(nn(1)-2.5)>0.001 then pause,end
47 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
48 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
49 if abs(nn(1)-2.53)>0.001 then pause,end
50
51 // Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
52 // pjac uses the macro res1 defined above.
53 function [r, ier] = psol(wp, iwp, b)
54     ier = 0;
55     //Compute the LU factorization of R.
56     sp = sparse(iwp, wp);
57     [h, rk] = lufact(sp);
58     //Solve the system LU*X = b
59     r = lusolve(h, b);
60     ludel(h);
61 endfunction
62 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
63     ires = 0;
64     SQuround = 1.490D-08;
65     tx = t;
66     nrow = 0;
67     e = zeros(1, neq);
68     wp = zeros(neq*neq, 1);
69     iwp = zeros(neq*neq, 2);
70     for i=1:neq
71         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
72         if h*ydot(i) < 0 then del = -del; end
73         ysave = y(i);
74         ypsave = ydot(i);
75         y(i) = y(i) + del;
76         ydot(i) = ydot(i) + cj*del;
77         [e ires]=res1(tx, y, ydot);
78         if ires < 0 then return; end
79         delinv = 1/del;
80         for j=1:neq
81             wp(nrow+j) = delinv*(e(j)-savr(j));
82             iwp(nrow+j,1) = i;
83             iwp(nrow+j,2) = j;
84         end
85         nrow = nrow + neq;
86         y(i) = ysave;
87         ydot(i) = ypsave;
88     end
89 endfunction
90
91 y0=1;t=2:6;t0=1;y0d=3;
92 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
93 atol=1.d-6;rtol=0;ng=2;
94 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
95 if abs(nn(1)-2.47)>0.001 then pause,end
96 y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
97 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
98 if abs(nn(1)-2.5)>0.001 then pause,end
99 y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
100 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
101 if abs(nn(1)-2.53)>0.001 then pause,end
102
103 //C
104 //C-----------------------------------------------------------------------
105 //C Second problem (Van Der Pol oscillator).
106 //C The initial value problem is..
107 //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
108 //C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
109 //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
110 //C The root function is  G = Y1.
111 //C An analytic solution is not known, but the zeros of Y1 are known
112 //C to 15 figures for purposes of checking the accuracy.
113 //C-----------------------------------------------------------------------
114 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
115 rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
116 t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
117 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
118 if abs(nn(1)-81.163512)>0.001 then pause,end
119
120 deff('[delta,ires]=res2(t,y,ydot)',...
121 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
122 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
123 deff('J=jac2(t,y,ydot,c)','y1=y(1);y2=y(2);J=[c,-1;200*y1*y2+1,c-100*(1-y1*y1)]')
124 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
125 deff('s=gr2(t,y,yd)','s=y(1)')
126 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
127
128 // Same problem, with psol and pjac example routines
129
130 info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
131 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,'psol1','pjac1');
132 if abs(nn(1)-81.163512)>0.009 then pause,end
133 deff('s=gr2(t,y,yd)','s=y(1)')
134 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,'psol1','pjac1');
135 if abs(nn(1)-81.163512)>0.009 then pause,end
136
137 // Same problem, with psol and pjac macros
138
139 // Redefine pjac to use res2
140 function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
141     ires = 0;
142     SQuround = 1.490D-08;
143     tx = t;
144     nrow = 0;
145     e = zeros(1, neq);
146     wp = zeros(neq*neq, 1);
147     iwp = zeros(neq*neq, 2);
148     for i=1:neq
149         del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
150         if h*ydot(i) < 0 then del = -del; end
151         ysave = y(i);
152         ypsave = ydot(i);
153         y(i) = y(i) + del;
154         ydot(i) = ydot(i) + cj*del;
155         [e ires]=res2(tx, y, ydot);
156         if ires < 0 then return; end
157         delinv = 1/del;
158         for j=1:neq
159             wp(nrow+j) = delinv*(e(j)-savr(j));
160             iwp(nrow+j,1) = i;
161             iwp(nrow+j,2) = j;
162         end
163         nrow = nrow + neq;
164         y(i) = ysave;
165         ydot(i) = ypsave;
166     end
167 endfunction
168 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,psol,pjac);
169 if abs(nn(1)-81.163512)>0.003 then pause,end
170 deff('s=gr2(t,y,yd)','s=y(1)')
171 [yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,psol,pjac);
172 if abs(nn(1)-81.163512)>0.003 then pause,end
173 info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
174
175 //           Hot Restart
176
177 [yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
178 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
179 [yy,nn,hotd]=daskr([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hotd);
180 if abs(nn(1)-162.57763)>0.004 then pause,end
181
182 //Test of Dynamic link (Require f77!)
183 //         1 making the routines
184 res22=[...
185 '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
186 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
187 '      INTEGER NEQ';
188 '      DIMENSION Y(*), YDOT(*), DELTA(*)';
189 '      NEQ=2';
190 '      CALL F2(NEQ,T,Y,DELTA)';
191 '      DO 10 I = 1,NEQ';
192 '         DELTA(I) = YDOT(I) - DELTA(I)';
193 ' 10   CONTINUE';
194 '      RETURN';
195 '      END';
196 '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
197 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
198 '      INTEGER NEQ';
199 '      DOUBLE PRECISION T, Y, YDOT';
200 '      DIMENSION Y(*), YDOT(*)';
201 '      YDOT(1) = Y(2)';
202 '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
203 '      RETURN';
204 '      END';]
205
206 jac22=[...
207 '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
208 ' ';
209 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
210 '      INTEGER  NROWPD';
211 '      DOUBLE PRECISION T, Y, PD';
212 '      PARAMETER (NROWPD=2)';
213 '      DIMENSION Y(2), PD(NROWPD,2)';
214 '      PD(1,1) = 0.0D0';
215 '      PD(1,2) = 1.0D0';
216 '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
217 '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
218 '      PD(1,1) = CJ - PD(1,1)';
219 '      PD(1,2) =    - PD(1,2)';
220 '      PD(2,1) =    - PD(2,1)';
221 '      PD(2,2) = CJ - PD(2,2)';
222 '      RETURN';
223 '      END';]
224
225
226 gr22=[...
227 '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
228 '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
229 '      INTEGER NEQ, NG';
230 '      DOUBLE PRECISION T, Y, GROOT';
231 '      DIMENSION Y(*), GROOT(*)';
232 '      GROOT(1) = Y(1)';
233 '      RETURN';
234 '      END';]
235
236 //Uncomment lines below: link may be machine dependent if some f77 libs are
237 //needed for linking
238 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
239 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
240 //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
241 //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
242 //unix_g('cd /tmp;make /tmp/jac22.o');
243 //          2  Linking the routines
244 //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
245 //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
246 //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
247 //          3 Calling the routines by dasrt
248 //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
249 // hot restart
250 //[yy,nn,hotd]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
251 //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
252 //[yy,nn,hotd]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hotd);