Fix f2c compilation 68/11568/8
Paul BIGNIER [Fri, 24 May 2013 09:28:16 +0000 (11:28 +0200)]
Also fix some numbering issues, making ddaskr suitable for more than one dimension problems

Change-Id: I335c648747ad3144ef657840284f5f123427fd93

scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c
scilab/modules/differential_equations/sci_gateway/fortran/bpjacd.f
scilab/modules/differential_equations/sci_gateway/fortran/bpsold.f
scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref
scilab/modules/differential_equations/tests/unit_tests/daskr.tst
scilab/modules/elementary_functions/includes/unsfdcopy.h
scilab/modules/elementary_functions/src/c/unsfdcopy.c

index fe292e2..01ae0fd 100644 (file)
@@ -46,8 +46,8 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
         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;
index 1634dd5..03a45f5 100644 (file)
@@ -24,8 +24,8 @@ c
       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
@@ -153,14 +153,26 @@ c+
 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
index 074599f..cbe8bd1 100644 (file)
@@ -24,8 +24,7 @@ c
       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
@@ -66,14 +65,31 @@ c
 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
index a8c62ef..edee5f2 100644 (file)
@@ -31,7 +31,7 @@ info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 [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;
@@ -93,7 +93,6 @@ if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
 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).
@@ -105,6 +104,7 @@ info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 //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);
@@ -116,6 +116,51 @@ 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)]'
 [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);
index 686ce31..2e91c81 100644 (file)
@@ -34,7 +34,7 @@ info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 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;
@@ -99,7 +99,6 @@ if abs(nn(1)-2.5)>0.001 then pause,end
 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-----------------------------------------------------------------------
@@ -112,6 +111,7 @@ info=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
 //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);
@@ -125,6 +125,53 @@ 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)]'
 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);
index f6058a2..d739520 100644 (file)
@@ -2,11 +2,11 @@
  * 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);
index 20649e8..ea80575 100644 (file)
 #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;
+}
 /*--------------------------------------------------------------------------*/