Also fix some numbering issues, making ddaskr suitable for more than one dimension problems
Change-Id: I335c648747ad3144ef657840284f5f123427fd93
for (j = 0 ; j < neq ; j++)
{
wp[nrow + j] = (e[j] - savr[j]) * delinv;
- iwp[nrow + j] = nrow + j + 1;
- iwp[nrow + j + neq * neq] = nrow + j + 1;
+ iwp[nrow + j] = i + 1;
+ iwp[nrow + j + neq * neq] = j + 1;
}
nrow += neq;
actual[i] = ysave;
double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
* wk(*), h, cj, wp(*), rpar(*)
integer ires, neq, iwp(*), ier, ipar(*)
- double precision dneq, diwp(2*neq*neq)
- integer vol,tops,nordre
+ double precision dneq
+ integer vol,tops,nordre,hsize
data nordre/5/,mlhs/3/
c
iadr(l)=l+l-1
c Transferring the output to Fortran
call btof(ier,1)
if(err.gt.0.or.err1.gt.0) return
- call btof(diwp,2*neq*neq)
+c call btof(iwp,2*neq*neq)
+ il2=iadr(lstk(top))
+ hsize=4
+ n=istk(il2+1)*istk(il2+2)*(istk(il2+3)+1)
+ liwp = 2*neq*neq
+c Test if the variable on the stack has same type and size as the theoretical iwp
+ if (istk(il2).ne.1.or.n.ne.liwp) then
+ call error(98)
+ return
+ endif
+ lx=sadr(il2+hsize)
+ do 900 i=1,liwp
+ iwp(i) = stk(lx+i-1)
+900 continue
+ top = top-1
if(err.gt.0.or.err1.gt.0) return
- call btof(wp,neq*neq)
+ lwp = neq*neq
+ call btof(wp,lwp)
if(err.gt.0.or.err1.gt.0) return
c+
- do 100 i=1, 2*neq*neq
- iwp(i) = diwp(i)
- 100 continue
c Normal return iero set to 0
iero=0
double precision t,y(*),ydot(*),savr(*),wk(*),cj,wght(*),wp(*),
* b(*),eplin,rpar(*)
integer neq,iwp(*),ier,ipar(*)
- integer vol,tops,nordre
- double precision diwp(2*neq*neq)
+ integer vol,tops,nordre,hsize
data nordre/4/,mlhs/2/
c
iadr(l)=l+l-1
c Minimum entry arguments for the simulator. The value of these arguments
c comes from the Fortran context (call list)
c+
- isize = 2*neq*neq
- do 100 i=1, isize
- diwp(i) = iwp(i)
- 100 continue
c
- call ftob(wp,neq*neq,istk(il+1))
+ lwp = neq*neq
+ call ftob(wp,lwp,istk(il+1))
if(err.gt.0.or.err1.gt.0) return
- call ftob(diwp,isize,istk(il+2))
+c call ftob(iwp,isize,istk(il+2))
+ ilx=iadr(lstk(istk(il+2)))
+ hsize=4
+ liwp = 2*neq*neq
+ if(top.ge.bot) then
+ call error(18)
+ return
+ endif
+ top=top+1
+ il2=iadr(lstk(top))
+ err=lstk(top)+sadr(hsize)+liwp-lstk(bot)
+ if(err.gt.0) then
+ call error(17)
+ return
+ endif
+ call icopy(hsize,istk(ilx),1,istk(il2),1)
+ l=sadr(il2+hsize)
+ do 900 i=1,liwp
+ stk(l+i-1) = iwp(i)
+900 continue
+ lstk(top+1)=l+liwp
if(err.gt.0.or.err1.gt.0) return
call ftob(b,neq,istk(il+3))
if(err.gt.0.or.err1.gt.0) return
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
// Same problem, but using macro for the derivative evaluation function 'res1'
-deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
+deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y')
deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
y0=1;t=2:6;t0=1;y0d=3;
atol=1.d-6;rtol=0;ng=2;
y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
-info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
//C
//C-----------------------------------------------------------------------
//C Second problem (Van Der Pol oscillator).
//C An analytic solution is not known, but the zeros of Y1 are known
//C to 15 figures for purposes of checking the accuracy.
//C-----------------------------------------------------------------------
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
deff('s=gr2(t,y,yd)','s=y(1)')
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
+// Same problem, with psol and pjac example routines
+info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,'psol1','pjac1');
+if abs(nn(1)-81.163512)>0.009 then bugmes();quit;end
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,'psol1','pjac1');
+if abs(nn(1)-81.163512)>0.009 then bugmes();quit;end
+// Same problem, with psol and pjac macros
+// Redefine pjac to use res2
+function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
+ ires = 0;
+ SQuround = 1.490D-08;
+ tx = t;
+ nrow = 0;
+ e = zeros(1, neq);
+ wp = zeros(neq*neq, 1);
+ iwp = zeros(neq*neq, 2);
+ for i=1:neq
+ del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
+ if h*ydot(i) < 0 then del = -del; end
+ ysave = y(i);
+ ypsave = ydot(i);
+ y(i) = y(i) + del;
+ ydot(i) = ydot(i) + cj*del;
+ [e ires]=res2(tx, y, ydot);
+ if ires < 0 then return; end
+ delinv = 1/del;
+ for j=1:neq
+ wp(nrow+j) = delinv*(e(j)-savr(j));
+ iwp(nrow+j,1) = i;
+ iwp(nrow+j,2) = j;
+ end
+ nrow = nrow + neq;
+ y(i) = ysave;
+ ydot(i) = ypsave;
+ end
+endfunction
+Warning : redefining function: pjac . Use funcprot(0) to avoid this message
+
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,psol,pjac);
+if abs(nn(1)-81.163512)>0.003 then bugmes();quit;end
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,psol,pjac);
+if abs(nn(1)-81.163512)>0.003 then bugmes();quit;end
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
// Hot Restart
[yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
if abs(nn(1)-2.53)>0.001 then pause,end
// Same problem, but using macro for the derivative evaluation function 'res1'
-deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
+deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2.*log(y)+8)./t-5).*y')
deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
y0=1;t=2:6;t0=1;y0d=3;
y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
if abs(nn(1)-2.53)>0.001 then pause,end
-info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
//C
//C-----------------------------------------------------------------------
//C An analytic solution is not known, but the zeros of Y1 are known
//C to 15 figures for purposes of checking the accuracy.
//C-----------------------------------------------------------------------
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
deff('s=gr2(t,y,yd)','s=y(1)')
[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
+// Same problem, with psol and pjac example routines
+
+info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,'psol1','pjac1');
+if abs(nn(1)-81.163512)>0.009 then pause,end
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,'psol1','pjac1');
+if abs(nn(1)-81.163512)>0.009 then pause,end
+
+// Same problem, with psol and pjac macros
+
+// Redefine pjac to use res2
+function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
+ ires = 0;
+ SQuround = 1.490D-08;
+ tx = t;
+ nrow = 0;
+ e = zeros(1, neq);
+ wp = zeros(neq*neq, 1);
+ iwp = zeros(neq*neq, 2);
+ for i=1:neq
+ del = max(SQuround*max(abs(y(i)), abs(h*ydot(i))), 1/rewt(i))
+ if h*ydot(i) < 0 then del = -del; end
+ ysave = y(i);
+ ypsave = ydot(i);
+ y(i) = y(i) + del;
+ ydot(i) = ydot(i) + cj*del;
+ [e ires]=res2(tx, y, ydot);
+ if ires < 0 then return; end
+ delinv = 1/del;
+ for j=1:neq
+ wp(nrow+j) = delinv*(e(j)-savr(j));
+ iwp(nrow+j,1) = i;
+ iwp(nrow+j,2) = j;
+ end
+ nrow = nrow + neq;
+ y(i) = ysave;
+ ydot(i) = ypsave;
+ end
+endfunction
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info,psol,pjac);
+if abs(nn(1)-81.163512)>0.003 then pause,end
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info,psol,pjac);
+if abs(nn(1)-81.163512)>0.003 then pause,end
+info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+
// Hot Restart
[yy,nn,hotd]=daskr([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2007 - INRIA - Sylvestre LEDRU
* Copyright (C) 2010 - DIGITEO - Allan CORNET
- *
+ *
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
- * are also available at
+ * are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
*
*/
/**
* <long-description>
*
- * @param n
- * @param dx
- * @param incx
- * @param dy
- * @param incy
+ * @param n
+ * @param dx
+ * @param incx
+ * @param dy
+ * @param incy
* @return <ReturnValue>
*/
ELEMENTARY_FUNCTIONS_IMPEXP int C2F(unsfdcopy)(int *n, long long *dx, int *incx, long long *dy, int *incy);
#include "machine.h"
#include "unsfdcopy.h"
/*--------------------------------------------------------------------------*/
- int C2F(unsfdcopy)(int *n, long long *dx, int *incx, long long *dy, int *incy)
+int C2F(unsfdcopy)(int *n, long long *dx, int *incx, long long *dy, int *incy)
{
- if (*n <= 0) return 0;
+ if (*n <= 0)
+ {
+ return 0;
+ }
- if ( (*incx == 1) && (*incy == 1) )
- {
- /* code for both increments equal to 1 */
- /* clean-up loop */
- memmove(dy , dx , (*n *sizeof(double)) );
- }
- else
- {
- int i = 0;
+ if ( (*incx == 1) && (*incy == 1) )
+ {
+ /* code for both increments equal to 1 */
+ /* clean-up loop */
+ memmove(dy , dx , (*n * sizeof(double)) );
+ }
+ else
+ {
+ int i = 0;
- /* code for unequal increments or equal increments */
- /* not equal to 1 */
- int ix = *incx >= 0 ? 0 : (1 - *n) * *incx;
- int iy = *incy >= 0 ? 0 : (1 - *n) * *incy ;
+ /* code for unequal increments or equal increments */
+ /* not equal to 1 */
+ int ix = *incx >= 0 ? 0 : (1 - *n) * *incx;
+ int iy = *incy >= 0 ? 0 : (1 - *n) * *incy ;
- for (i = 0; i < *n; i++)
- {
- dy[iy] = dx[ix];
- ix += *incx;
- iy += *incy;
- }
- }
- return 0;
-}
+ for (i = 0; i < *n; i++)
+ {
+ dy[iy] = dx[ix];
+ ix += *incx;
+ iy += *incy;
+ }
+ }
+ return 0;
+}
/*--------------------------------------------------------------------------*/
projects / scilab.git / commitdiff