ODEPACK solvers (daskr): adding Scilab macro example for pjac and psol 72/11372/14
Paul BIGNIER [Tue, 23 Apr 2013 15:25:20 +0000 (17:25 +0200)]
Change-Id: Iada46d9fc3025c75e6dd17ef8fb3e1ada536552b

scilab/modules/differential_equations/help/en_US/daskr.xml
scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c
scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref
scilab/modules/differential_equations/tests/unit_tests/daskr.tst

index e762e2a..a745f4b 100644 (file)
@@ -17,7 +17,7 @@
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, jac]] [, hd])</synopsis>
+        <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
@@ -475,7 +475,7 @@ integer ny,ng
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>[r,ier] = psol(neq,R,iR,b)</literal> and the
+                                <literal>[r, ier] = psol(wp, iwp, b)</literal> and the
                                 <literal>psol</literal> function must return the solution of that system in
                                 <literal>r</literal> and an error flag <literal>ier</literal>.
                             </para>
@@ -490,7 +490,7 @@ list(psol,x1,x2,...)
                                 <literal>psol</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-psol(R,iR,b,x1,x2,...)
+psol(wp,iwp,b,x1,x2,...)
  ]]></programlisting>
                             <para>
                                 <literal>psol</literal> still returns the solution in <literal>r</literal>.
@@ -532,8 +532,8 @@ integer neq,iwp(*),ier,ipar(*)
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>[R,iR,ier] = pjac(neq,t,y,ydot,h,cj)</literal> and in return,
-                                the arrays R and iR must contain all preconditioner information in compressed sparse row format.
+                                <literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> and in return,
+                                the arrays wp and iwp must contain all preconditioner information in compressed sparse row format.
                             </para>
                         </listitem>
                         <listitem>
@@ -546,12 +546,12 @@ list(pjac,x1,x2,...)
                                 <literal>pjac</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-pjac(neq,t,y,ydot,h,cj,x1,x2,...)
+pjac(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)
  ]]></programlisting>
                             <para>
                                 <literal>pjac</literal> still returns factorised
                                 <literal>dg/dy + cj*dg/dydot</literal> as a function of
-                                <literal>(neq,t,y,ydot,h,cj,x1,x2,...)</literal>.
+                                <literal>(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)</literal>.
                             </para>
                         </listitem>
                         <listitem>
index ddfc066..fe292e2 100644 (file)
@@ -20,14 +20,14 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
     double ysave = 0;
     double ypsave = 0;
     double* e = NULL;
-
+    int neq = *nequations;
     double SQuround = sqrt(C2F(dlamch)("P"));
 
     tx = *tOld;
 
-    e = (double*)calloc(*nequations, sizeof(double));
+    e = (double*)calloc(neq, sizeof(double));
 
-    for (i = 0 ; i < *nequations ; ++i)
+    for (i = 0 ; i < neq ; ++i)
     {
         del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
         del *= (*h * actualP[i] >= 0) ? 1 : -1;
@@ -43,15 +43,15 @@ int C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *ac
         }
 
         delinv = 1. / del;
-        for (j = 0 ; j < *nequations ; j++)
+        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;
         }
-        nrow       += *nequations;
+        nrow       += neq;
         actual[i]  =  ysave;
         actualP[i] =  ypsave;
-        iwp[i] = i;
-        iwp[i + *nequations] = i;
     }
 
     free(e);
@@ -67,7 +67,7 @@ int C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
     int info = 0;
     int *ipiv = NULL;
 
-    ipiv = (int*)malloc(sizeof(int) * *nequations);
+    ipiv = (int *) malloc(*nequations * sizeof(int));
 
     C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
 
index 6bb8266..a8c62ef 100644 (file)
@@ -8,8 +8,7 @@
 // The terms are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 // =============================================================================
-// Execute with exec("SCI/modules/differential_equations/tests/unit_tests/daskr.tst");
-//  or test_run('differential_equations', 'daskr', ['no_check_error_output']);
+// Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
 //C-----------------------------------------------------------------------
 //C First problem.
 //C The initial value problem is..
@@ -31,6 +30,7 @@ y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
 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('[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;
@@ -43,6 +43,57 @@ 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);
 if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
+// Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
+// pjac uses the macro res1 defined above.
+function [r, ier] = psol(wp, iwp, b)
+    ier = 0;
+    //Compute the LU factorization of R.
+    sp = sparse(iwp, wp);
+    [h, rk] = lufact(sp);
+    //Solve the system LU*X = b
+    r = lusolve(h, b);
+    ludel(h);
+endfunction
+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]=res1(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
+y0=1;t=2:6;t0=1;y0d=3;
+info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
+if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
+y0=yy(2,2);y0d=yy(3,2);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.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).
@@ -213,7 +264,7 @@ gr22=[...
 !      RETURN                                              !
 !                                                          !
 !      END                                                 !
-//Uncomment lines below: link may be machine dependent if some f77 libs are 
+//Uncomment lines below: link may be machine dependent if some f77 libs are
 //needed for linking
 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
 //unix_g('cd /tmp;rm -f /tmp/jac22.f');
index dbbe1cd..686ce31 100644 (file)
@@ -1,4 +1,3 @@
-
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
@@ -10,8 +9,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 // =============================================================================
 
-// Execute with exec("SCI/modules/differential_equations/tests/unit_tests/daskr.tst");
-//  or test_run('differential_equations', 'daskr', ['no_check_error_output']);
+// Run with test_run('differential_equations', 'daskr', ['no_check_error_output']);
 
 //C-----------------------------------------------------------------------
 //C First problem.
@@ -35,6 +33,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 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('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
 
@@ -49,6 +48,59 @@ 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);
 if abs(nn(1)-2.53)>0.001 then pause,end
 
+// Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
+// pjac uses the macro res1 defined above.
+function [r, ier] = psol(wp, iwp, b)
+    ier = 0;
+    //Compute the LU factorization of R.
+    sp = sparse(iwp, wp);
+    [h, rk] = lufact(sp);
+    //Solve the system LU*X = b
+    r = lusolve(h, b);
+    ludel(h);
+endfunction
+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]=res1(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
+
+y0=1;t=2:6;t0=1;y0d=3;
+info=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=daskr([y0,y0d],t0,t,atol,rtol,res1,ng,'gr1',info,psol,pjac);
+if abs(nn(1)-2.47)>0.001 then pause,end
+y0=yy(2,2);y0d=yy(3,2);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.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-----------------------------------------------------------------------
 //C Second problem (Van Der Pol oscillator).
@@ -134,7 +186,7 @@ gr22=[...
 '      RETURN';
 '      END';]
 
-//Uncomment lines below: link may be machine dependent if some f77 libs are 
+//Uncomment lines below: link may be machine dependent if some f77 libs are
 //needed for linking
 //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
 //unix_g('cd /tmp;rm -f /tmp/jac22.f');