Differential_equations: integrate daskr into dae() 87/12387/7
Paul BIGNIER [Fri, 30 Aug 2013 09:19:54 +0000 (11:19 +0200)]
Step 3: integrate to dae()

Change-Id: I43dad4a97c340a80f7f78080cbe6723a62fad7b2

scilab/CHANGES_5.5.X
scilab/modules/differential_equations/help/en_US/dae.xml
scilab/modules/differential_equations/help/en_US/daeoptions.xml
scilab/modules/differential_equations/macros/dae.sci
scilab/modules/differential_equations/macros/daeoptions.sci
scilab/modules/differential_equations/sci_gateway/fortran/sci_f_daskr.f
scilab/modules/differential_equations/tests/unit_tests/dae.dia.ref
scilab/modules/differential_equations/tests/unit_tests/dae.tst
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/dae_2.png [new file with mode: 0644]

index 0aeb0bc..ae208fa 100644 (file)
@@ -122,6 +122,7 @@ Xcos
 ====
 
 * New DAE solver: DDaskr, using BDF methods with direct Newton and preconditioned Krylov linear solvers, which includes rootfinding.
+  It is available from the dae() function.
 
 * In Modelica initialization GUI, inputs (eg. sensor) were not handled.
 
index 54d6ba1..a6deec6 100644 (file)
@@ -1,6 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier : added "roots2" (daskr)
  * Copyright (C) 2008 - INRIA
  * ...
  *
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis> y=dae(initial,t0,t,res)
-            [y [,hd]]=dae(initial,t0,t [,rtol, [atol]],res [,jac] [,hd])
-            [y,rd]=dae("root",initial,t0,t,res,ng,surface)
-            [y ,rd [,hd]]=dae("root",initial,t0,t [,rtol, [atol]],res [,jac], ng, surface [,hd])
+        <synopsis> y = dae(initial, t0, t, res)
+            [y [,hd]] = dae(initial, t0, t [,rtol, [atol]], res [,jac] [,hd])
+            [y, rd] = dae("root", initial, t0, t, res, ng, surface)
+            [y, rd [,hd]] = dae("root", initial, t0, t [,rtol, [atol]], res [,jac], ng, surface [,hd])
+            [y, rd] = dae("root2", initial, t0, t, res, ng, surface)
+            [y, rd [,hd]] = dae("root2", initial, t0, t [,rtol, [atol]], res [,jac], ng, surface [, psol, pjac] [, hd])
         </synopsis>
     </refsynopsisdiv>
     <refsection>
@@ -398,6 +401,22 @@ r=surface(t,x,p1,p2,...)
                 </listitem>
             </varlistentry>
             <varlistentry>
+                <term>pjac</term>
+                <listitem>
+                    <para>
+                        external (function, list or string). See <link linkend="daskr">daskr</link>
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>psol</term>
+                <listitem>
+                    <para>
+                        external (function, list or string). See <link linkend="daskr">daskr</link>
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
                 <term>hd</term>
                 <listitem>
                     <para>a real vector, as an output it stores the
@@ -428,12 +447,17 @@ r=surface(t,x,p1,p2,...)
         <title>Description</title>
         <para>
             The <literal>dae</literal> function is a gateway built above the
-            <link linkend="dassl">dassl</link> and <link linkend="dasrt">dasrt</link>
-            function designed for implicit differential equations integration.
+            <link linkend="dassl">dassl</link>, <link linkend="dasrt">dasrt</link>
+            and <link linkend="daskr">daskr</link>
+            functions designed for implicit differential equations integration.
+        </para>
+        <para>
+            Option <literal>"root"</literal> calls the <link linkend="dasrt">dasrt</link> routine,
+            and <literal>"root2"</literal> calls <link linkend="dasrt">daskr</link>.
         </para>
         <programlisting role="no-scilab-exec"><![CDATA[
-g(t,x,xdot)=0
-x(t0)=x0  and   xdot(t0)=xdot0
+g(t, x, xdot) = 0
+x(t0) = x0 and xdot(t0) = xdot0
  ]]></programlisting>
         <para>
             If <literal>xdot0</literal> is not given in the <emphasis>initial</emphasis>
@@ -447,57 +471,162 @@ x(t0)=x0  and   xdot(t0)=xdot0
             <link linkend="daeoptions">%DAEOPTIONS</link>(7) must be set to 1.
         </para>
         <para>Detailed examples using Scilab and C coded externals are given in
-            <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal>
+            <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal> and
+            <literal>modules/differential_equations/tests/unit_tests/daskr.tst</literal>
         </para>
     </refsection>
     <refsection>
         <title>Examples</title>
+        <para>
+            Example #1: dassl (no root finding)
+        </para>
         <programlisting role="example"><![CDATA[
-//Example with Scilab code
-function [r,ires]=chemres(t,y,yd)
+// Example with Scilab code
+--------------------------------------------------
+function [r, ires] = chemres(t, y, yd)
     r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
     r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
     r(3) =       y(1) +     y(2)      + y(3)-1;
     ires =  0;
 endfunction
 
-function pd=chemjac(x,y,yd,cj)
-    pd=[-0.04-cj , 1d4*y(3)               , 1d4*y(2);
-         0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
-         1       , 1                      , 1       ]
+function pd = chemjac(x, y, yd, cj)
+    pd = [-0.04-cj , 1d4*y(3)               , 1d4*y(2);
+           0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
+           1       , 1                      , 1       ]
 endfunction
 
-x0=[1; 0; 0];
-xd0=[-0.04; 0.04; 0];
-t=[1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
+x0 = [1; 0; 0];
+xd0 = [-0.04; 0.04; 0];
+t = [1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
 
-y=dae([x0,xd0],0,t,chemres);// returns requested observation time points
+y = dae([x0, xd0], 0, t, chemres); // Returns requested observation time points
 
-%DAEOPTIONS=list([],1,[],[],[],0,0); // ask  dae mesh points to be returned
-y=dae([x0,xd0],0,4d10,chemres); // without jacobian
-y=dae([x0,xd0],0,4d10,chemres,chemjac); // with jacobian
+%DAEOPTIONS = list([], 1, [], [], [], 0, 0); // Ask  dae mesh points to be returned
+y = dae([x0, xd0], 0, 4d10, chemres); // Without jacobian
+y = dae([x0, xd0], 0, 4d10, chemres, chemjac); // With jacobian
  ]]></programlisting>
+        <para>
+            Example #2: dasrt ("root")
+        </para>
         <programlisting role="example"><![CDATA[
+// Example with C code (C compiler needed)
+--------------------------------------------------
+bOK = haveacompiler();
+if bOK <> %t
+    [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
+    return
+end
+
+//-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
+code = ['#include <math.h>'
+      'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
+      '{res[0] = yd[0] - y[1];'
+      ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
+      ' '
+      'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
+      '{pd[0] = *cj - 0.0;'
+      ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
+      ' pd[2] =     - 1.0;'
+      ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
+      ' '
+      'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
+      '{ groot[0] = y[0];}']
+previous_dir = pwd();
+cd TMPDIR;
+mputl(code, 't22.c')
+
+//-2- Compile and load them
+ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
+exec('t22loader.sce')
+
+//-3- Run
+rtol = [1.d-6; 1.d-6];
+atol = [1.d-6; 1.d-4];
+t0 = 0; t = [20:20:200];
+y0 = [2; 0]; y0d = [0; -2];
+ng = 1;
 
-//example with C code (C compiler needed)
+// Simple simulation
+t = 0:0.003:300;
+yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
+clf(); plot(yy(1, :), yy(2, :))
+// Find first point where yy(1) = 0
+[yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
+
+// Hot restart for next point
+t01 = nn(1);
+[pp, qq] = size(yy);
+y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
+[yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
+cd(previous_dir);
+ ]]></programlisting>
+        <scilab:image><![CDATA[
+code = ['#include <math.h>'
+      'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
+      '{res[0] = yd[0] - y[1];'
+      ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
+      ' '
+      'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
+      '{pd[0] = *cj - 0.0;'
+      ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
+      ' pd[2] =     - 1.0;'
+      ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
+      ' '
+      'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
+      '{ groot[0] = y[0];}']
+previous_dir = pwd();
+cd TMPDIR;
+mputl(code, 't22.c')
+ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
+exec('t22loader.sce')
+rtol = [1.d-6; 1.d-6];
+atol = [1.d-6; 1.d-4];
+t0 = 0; t = [20:20:200];
+y0 = [2; 0]; y0d = [0; -2];
+ng = 1;
+t = 0:0.003:300;
+yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
+clf(); plot(yy(1, :), yy(2, :))
+[yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
+t01 = nn(1);
+[pp, qq] = size(yy);
+y01 = yy(2:3, qq);
+y0d1 = yy(3:4, qq);
+[yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
+cd(previous_dir);
+ ]]></scilab:image>
+        <para>
+            Example #3: daskr ("root2"), using default 'psol' and 'pjac' routines
+        </para>
+        <programlisting role="example"><![CDATA[
+// Example with C code (C compiler needed)
 --------------------------------------------------
-bOK=haveacompiler();
-if bOK<>%t
-    [btn] = messagebox(["You need a C compiler for this example.";"Execution of this example is canceled."], "Software problem", 'info');
+bOK = haveacompiler();
+if bOK <> %t
+    [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
     return
 end
 
-//-1- create the C codes in TMPDIR - Vanderpol equation, implicit form
-code=['#include <math.h>'
-      'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
+//-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
+code = ['#include <math.h>'
+      'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
       '{res[0] = yd[0] - y[1];'
       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
       ' '
-      'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
-      '{pd[0]=*cj - 0.0;'
-      ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
-      ' pd[2]=    - 1.0;'
-      ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
+      'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
+      '{pd[0] = *cj - 0.0;'
+      ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
+      ' pd[2] =     - 1.0;'
+      ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
       ' '
       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
       '{ groot[0] = y[0];}']
@@ -505,75 +634,74 @@ previous_dir = pwd();
 cd TMPDIR;
 mputl(code, 't22.c')
 
-//-2- compile and load them
-ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',[],'t22loader.sce');
+//-2- Compile and load them
+ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
 exec('t22loader.sce')
 
-//-3- run
-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;
+//-3- Run
+rtol = [1.d-6; 1.d-6];
+atol = [1.d-6; 1.d-4];
+t0 = 0; t = [20:20:200];
+y0 = [2; 0]; y0d = [0; -2];
+ng = 1;
 
-//simple simulation
-t=0:0.003:300;
-yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
-clf();plot(yy(1,:),yy(2,:))
-//find first point where yy(1)=0
-[yy,nn,hotd]=dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
-plot(yy(1,1),yy(2,1),'r+')
-xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
+// Simple simulation
+t = 0:0.003:300;
+yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
+clf(); plot(yy(1, :), yy(2, :))
+// Find first point where yy(1) = 0
+%DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
+[yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
 
-//hot restart for next point
-t01=nn(1);
-[pp,qq]=size(yy);
-y01=yy(2:3,qq);
-y0d1=yy(3:4,qq);
-[yy,nn,hotd]=dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
-plot(yy(1,1),yy(2,1),'r+')
-xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
+// Hot restart for next point
+t01 = nn(1);
+[pp, qq] = size(yy);
+y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
+[yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
 cd(previous_dir);
  ]]></programlisting>
         <scilab:image><![CDATA[
-code=['#include <math.h>'
-      'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
+code = ['#include <math.h>'
+      'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
       '{res[0] = yd[0] - y[1];'
       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
       ' '
-      'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
-      '{pd[0]=*cj - 0.0;'
-      ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
-      ' pd[2]=    - 1.0;'
-      ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
+      'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
+      '{pd[0] = *cj - 0.0;'
+      ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
+      ' pd[2] =     - 1.0;'
+      ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
       ' '
       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
       '{ groot[0] = y[0];}']
 previous_dir = pwd();
 cd TMPDIR;
 mputl(code, 't22.c')
-ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',[],'t22loader.sce');
+ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
 exec('t22loader.sce')
-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;
-t=0:0.003:300;
-yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
-clf();plot(yy(1,:),yy(2,:))
-[yy,nn,hotd]=dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
-plot(yy(1,1),yy(2,1),'r+')
-xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
-t01=nn(1);
-[pp,qq]=size(yy);
-y01=yy(2:3,qq);
-y0d1=yy(3:4,qq);
-[yy,nn,hotd]=dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
-plot(yy(1,1),yy(2,1),'r+')
-xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)));
+rtol = [1.d-6; 1.d-6];
+atol = [1.d-6; 1.d-4];
+t0 = 0; t = [20:20:200];
+y0 = [2; 0]; y0d = [0; -2];
+ng = 1;
+t = 0:0.003:300;
+yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
+clf(); plot(yy(1, :), yy(2, :))
+%DAEOPTIONS = list([], 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
+[yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
+t01 = nn(1);
+[pp, qq] = size(yy);
+y01 = yy(2:3, qq);
+y0d1 = yy(3:4, qq);
+[yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
+plot(yy(1, 1), yy(2, 1), 'r+')
+xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
 cd(previous_dir);
  ]]></scilab:image>
     </refsection>
@@ -590,6 +718,9 @@ cd(previous_dir);
                 <link linkend="dassl">dassl</link>
             </member>
             <member>
+                <link linkend="dasrt">dasrt</link>
+            </member>
+            <member>
                 <link linkend="daskr">daskr</link>
             </member>
             <member>
index 740204a..d59d9d9 100644 (file)
@@ -1,13 +1,14 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier : added daskr options
  * Copyright (C) 2008 - INRIA
  * ...
- * 
+ *
  * 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.1-en.txt
  *
  -->
     </refsynopsisdiv>
     <refsection>
         <title>Description</title>
-        <para>If it exists in the dae function calling context the variable
-            <literal>%DAEOPTIONS</literal> the dae function use it to sets
-            options.
+        <para>
+            Sets a variable <literal>%DAEOPTIONS</literal> that is used in function <literal>dae</literal>
+            to set solver options.
         </para>
-        <para>This daeoptions function interactively displays a command which
-            should be executed to set various options of the <link linkend="dae">dae</link> solver.
+        <para>This daeoptions function interactively displays a window which
+            sets various options of the <link linkend="dae">dae</link> solver.
         </para>
         <para>
             <caution>
         <para>
             The variable <literal>%DAEOPTIONS</literal> is a <link linkend="list">list</link> with the following elements:
         </para>
-        <programlisting role="no-scilab-exec"><![CDATA[  
+        <programlisting role="no-scilab-exec"><![CDATA[
 list(tstop,imode,band,maxstep,stepin,nonneg,isest)
  ]]></programlisting>
         <para>The default value is:</para>
-        <programlisting role="no-scilab-exec"><![CDATA[  
+        <programlisting role="no-scilab-exec"><![CDATA[
 list([],0,[],[],[],0,0)
  ]]></programlisting>
         <para>The meaning of the elements is described below.</para>
         <variablelist>
             <varlistentry>
+                <term>solver</term>
+                <listitem>
+                    <para>
+                        if it is 0, <literal>dae</literal> will use the dassl/dasrt solver,
+                        and if it is 1, <literal>dae</literal> will use daskr.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
                 <term>tstop</term>
                 <listitem>
                     <para>a real scalar or an empty matrix, gives the maximum time for
@@ -63,8 +73,8 @@ list([],0,[],[],[],0,0)
                 <term>imode</term>
                 <listitem>
                     <para>
-                        if it is 0 <literal>dae</literal> returns only the user specified time point
-                        values, and if it is 1 <literal>dae</literal> returns its intermediate computed
+                        if it is 0, <literal>dae</literal> returns only the user specified time point
+                        values, and if it is 1, <literal>dae</literal> returns its intermediate computed
                         values.
                     </para>
                 </listitem>
@@ -116,8 +126,91 @@ list([],0,[],[],[],0,0)
                 <term>isest</term>
                 <listitem>
                     <para>a scalar, must be set to 0 if the given initial condition is
-                        compatible: <literal>g(t0,x0,xdot0)=0</literal>, and set to 1 if
-                        <literal>xdot0</literal> is just an estimation.
+                        compatible: <literal>g(t0,x0,xdot0)=0</literal>. If
+                        <literal>xdot0</literal> is just an estimation, then set to
+                        1 if you are using dassl, or set to
+                        <literal>[+-1, ..., +-1]</literal>, with:
+                        <literal>1</literal> if y(i) is a differential variable and
+                        <literal>-1</literal> if y(i) is an algebraic variable
+                        (if its derivatives do not appear explicitly in the function g(t, y, ydot)).
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>method</term>
+                <listitem>
+                    <para>
+                        a scalar, direct / Krylov method. Set to <literal>1</literal> and provide a routine
+                        <literal>psol</literal> in <literal>dae()</literal>
+                        if you want the solver to use Krylov iterations, else (daskr's direct method) set to
+                        <literal>0</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>Kry_params</term>
+                <listitem>
+                    <para>
+                        a vector, treat as dummy argument if you have set method=0. Otherwise, set to <literal>[]</literal>
+                        or to <literal>[maxl kmp nrmax epli]</literal>, where:
+                        - <literal>maxl </literal> = maximum number of iterations in the GMRes algorithm (default min(5, neq)),
+                        - <literal>kmp  </literal> = number of vectors on which orthogonalization is done in the GMRes algorithm (default maxl),
+                        - <literal>nrmax</literal> = maximum number of restarts of the GMRes algorithm per nonlinear iteration (default 5),
+                        - <literal>epli </literal> = convergence test constant in GMRes algorithm (default 0.05).
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>init</term>
+                <listitem>
+                    <para>
+                        a scalar, treat as dummy argument if you have set isest=0. Set to 1 if the solver should stop right after
+                        computation of the initial conditions, else set to 0.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>precond</term>
+                <listitem>
+                    <para>
+                        a scalar, preconditioner computation and LU-factorization routine for <literal>psol</literal>.
+                        Treat as dummy argument if method=0. Set to 1 and provide a <literal>pjac</literal>
+                        routine in <literal>dae</literal> if the external <literal>psol</literal>
+                        should use a specific routine, else set to 0.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>control</term>
+                <listitem>
+                    <para>
+                        a vector, treat as dummy argument if you have set isest=0if you wish to control errors locally on all the variables then set to [].
+                        Otherwise, set to [+-1, ..., +-1], 1 if <literal>y(i)</literal> is a differential variable and
+                        -1 if <literal>y(i)</literal> is an algebraic variable
+                        (if its derivatives do not appear explicitly in the function <literal>g(t, y, ydot)</literal>).
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>heuristic</term>
+                <listitem>
+                    <para>
+                        a vector, heuristic parameters. Treat as dummy argument if isest=0.
+                        Otherwise, set to [] or to <literal>[mxnit mxnj mxnh lsoff stptol epinit]</literal>, where:
+                        - <literal>mxnit</literal>  = maximum number of Newton iterations per Jacobian or preconditioner evaluation (default 5 if info(8)=0, 15 otherwise),
+                        - <literal>mxnj</literal>   = maximum number of Jacobian or preconditioner evaluations (default 6 if info(8)=0, 2 otherwise),
+                        - <literal>mxnh</literal>   = maximum number of values of the artificial stepsize parameter h to be tried if info(7) ≠ [] (default 5),
+                        - <literal>lsoff</literal>  = flag to turn off the linesearch algorithm (lsoff = 0 means linesearch is on, lsoff = 1 means it is turned off) (default 0),
+                        - <literal>stptol</literal> = minimum scaled step in linesearch algorithm (default (unit roundoff)^(2/3)),
+                        - <literal>epinit</literal> = swing factor in the Newton iteration convergence test (default 0.01).
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>verbosity</term>
+                <listitem>
+                    <para>
+                        a scalar, verbosity. Set to 1 for minimal extra printing, 2 for full printing, else set to 0.
                     </para>
                 </listitem>
             </varlistentry>
index ecb3932..4e90552 100644 (file)
@@ -1,4 +1,5 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
 // Copyright (C) 2008 - INRIA - Sabine GAUZERE
 //
 // This file must be used under the terms of the CeCILL.
 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 
 
-function [varargout]=dae(varargin)
-    [lhs,rhs]=argn();
-    if (isdef("%DAEOPTIONS")==%F),
-        %DAEOPTIONS=list ([],0,[],[],[],0,0);
+function [varargout] = dae(varargin)
+
+    [lhs, rhs] = argn();
+
+    if isdef("%DAEOPTIONS") == %F then
+        if varargin(1) <> "root2" then
+            %DAEOPTIONS = list ([], 0, [], [], [], 0, 0);
+        else
+            %DAEOPTIONS = list ([], 0, [], [], [], 0, [], 0, [], 0, [], [], [], 1);
+        end
     end
-    if (size(varargin) == 0),
-        error(sprintf(gettext("%s: Wrong number of input argument(s): %d or %d expected.\n"), "dae",3,8));
+    if size(varargin) == 0 then
+        error(sprintf(gettext("%s: Wrong number of input argument(s): %d or %d expected.\n"), "dae", 3, 8));
     end
-    if (type(varargin(1))==1), //standard case (dassl)
-
-        if (rhs==4), //call without optional arguments
-            [x0,t0,t,res]=varargin(:)
-            if (lhs==2),
-                [y,hd]=dassl(x0,t0,t,res,%DAEOPTIONS);
-            elseif (lhs==1),
-                [y]=dassl(x0,t0,t,res,%DAEOPTIONS);
-            else,
-                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+    if type(varargin(1)) == 1 then //standard case (dassl)
+
+        if rhs == 4 then //call without optional arguments
+            [x0, t0, t, res] = varargin(:)
+            if lhs == 2 then
+                [y, hd] = dassl(x0, t0, t, res, %DAEOPTIONS);
+            elseif lhs == 1 then
+                [y] = dassl(x0, t0, t, res, %DAEOPTIONS);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
             end
 
-        elseif (rhs==8), // call with all the optional arguments
-            [x0,t0,t,rtol,atol,res,jac,hd]=varargin(:)
-            if (lhs==2),
-                [y,hd]=dassl(x0,t0,t,atol,rtol,res,jac,%DAEOPTIONS,hd);
-            elseif (lhs==1),
-                [y]=dassl(x0,t0,t,atol,rtol,res,jac,%DAEOPTIONS,hd);
-            else,
-                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+        elseif rhs == 8 then // call with all the optional arguments
+            [x0, t0, t, rtol, atol, res, jac, hd] = varargin(:)
+            if lhs == 2 then
+                [y, hd] = dassl(x0,t0, t, atol, rtol, res, jac, %DAEOPTIONS, hd);
+            elseif lhs == 1 then
+                [y] = dassl(x0, t0, t, atol, rtol, res, jac, %DAEOPTIONS, hd);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
             end
 
-        elseif (rhs==5),
+        elseif rhs == 5 then
 
-            if (type(varargin(4))==1),
-                [x0,t0,t,atol,res]=varargin(:)
-                if (lhs==2),
-                    [y,hd]=dassl(x0,t0,t,atol,res,%DAEOPTIONS);
-                elseif (lhs==1),
-                    [y]=dassl(x0,t0,t,atol,res,%DAEOPTIONS);
-                else,
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+            if type(varargin(4)) == 1 then
+                [x0, t0, t, atol, res] = varargin(:)
+                if lhs == 2 then
+                    [y, hd] = dassl(x0, t0, t, atol, res, %DAEOPTIONS);
+                elseif lhs == 1 then
+                    [y] = dassl(x0, t0, t, atol, res, %DAEOPTIONS);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                 end
 
-            else,
+            else
 
-                if (type(varargin(5))==1),
-                    [x0,t0,t,res,hd]=varargin(:)
-                    if (lhs==2),
-                        [y,hd]=dassl(x0,t0,t,res,%DAEOPTIONS,hd);
-                    elseif (lhs==1),
-                        [y]=dassl(x0,t0,t,res,%DAEOPTIONS,hd);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                if type(varargin(5)) == 1 then
+                    [x0, t0, t, res, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, hd] = dassl(x0, t0, t, res, %DAEOPTIONS, hd);
+                    elseif lhs == 1 then
+                        [y] = dassl(x0, t0, t, res, %DAEOPTIONS, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                     end
 
-                else,
-                    [x0,t0,t,res,jac]=varargin(:)
-                    if (lhs==2),
-                        [y,hd]=dassl(x0,t0,t,res,jac,%DAEOPTIONS);
-                    elseif (lhs==1),
-                        [y]=dassl(x0,t0,t,res,jac,%DAEOPTIONS);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                else
+                    [x0, t0, t, res, jac] = varargin(:)
+                    if lhs == 2 then
+                        [y, hd] = dassl(x0, t0, t, res, jac, %DAEOPTIONS);
+                    elseif lhs == 1 then
+                        [y] = dassl(x0, t0, t, res, jac, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                     end
 
                 end
 
             end
 
-        elseif(rhs==6),
+        elseif rhs == 6 then
 
-            if (type(varargin(4))==1),
+            if type(varargin(4)) == 1 then
 
-                if(type(varargin(5))==1),
-                    [x0,t0,t,rtol,atol,res]=varargin(:)
-                    if (lhs==2),
-                        [y,hd]=dassl(x0,t0,t,atol,rtol,res,%DAEOPTIONS);
-                    elseif (lhs==1),
-                        [y]=dassl(x0,t0,t,atol,rtol,res,%DAEOPTIONS);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                if type(varargin(5)) == 1 then
+                    [x0, t0, t, rtol, atol, res] = varargin(:)
+                    if lhs == 2 then
+                        [y, hd] = dassl(x0, t0, t, atol, rtol, res, %DAEOPTIONS);
+                    elseif lhs == 1 then
+                        [y] = dassl(x0, t0, t, atol, rtol, res, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                     end
 
-                else,
+                else
 
-                    if (type(varargin(6))==1),
-                        [x0,t0,t,atol,res,hd]=varargin(:)
-                        if (lhs==2),
-                            [y,hd]=dassl(x0,t0,t,atol,res,%DAEOPTIONS,hd);
-                        elseif (lhs==1),
-                            [y]=dassl(x0,t0,t,atol,res,%DAEOPTIONS,hd);
-                        else,
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                    if type(varargin(6)) == 1 then
+                        [x0, t0, t, atol, res, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, hd] = dassl(x0, t0, t, atol, res, %DAEOPTIONS, hd);
+                        elseif lhs == 1 then
+                            [y] = dassl(x0, t0, t, atol, res, %DAEOPTIONS, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                         end
 
-                    else,
-                        [x0,t0,t,atol,res,jac]=varargin(:)
-                        if (lhs==2),
-                            [y,hd]=dassl(x0,t0,t,atol,res,jac,%DAEOPTIONS);
-                        elseif (lhs==1),
-                            [y]=dassl(x0,t0,t,atol,res,jac,%DAEOPTIONS);
-                        else,
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                    else
+                        [x0, t0, t, atol, res, jac] = varargin(:)
+                        if lhs == 2 then
+                            [y, hd] = dassl(x0, t0, t, atol, res, jac, %DAEOPTIONS);
+                        elseif lhs == 1 then
+                            [y] = dassl(x0, t0, t, atol, res, jac, %DAEOPTIONS);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                         end
 
                     end
                 end
 
-            else,
-                [x0,t0,t,res,jac,hd]=varargin(:)
-                if (lhs==2),
-                    [y,hd]=dassl(x0,t0,t,res,jac,%DAEOPTIONS,hd);
-                elseif (lhs==1),
-                    [y]=dassl(x0,t0,t,res,jac,%DAEOPTIONS,hd);
-                else,
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+            else
+                [x0, t0, t, res, jac, hd] = varargin(:)
+                if lhs == 2 then
+                    [y, hd] = dassl(x0, t0, t, res, jac, %DAEOPTIONS, hd);
+                elseif lhs == 1 then
+                    [y] = dassl(x0, t0, t, res, jac, %DAEOPTIONS, hd);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                 end
 
             end
 
-        elseif (rhs==7),
+        elseif rhs == 7 then
 
-            if(type(varargin(5))==1),
+            if type(varargin(5)) == 1 then
 
-                if (type(varargin(7))==1),
-                    [x0,t0,t,rtol,atol,res,hd]=varargin(:)
-                    if (lhs==2),
-                        [y,hd]=dassl(x0,t0,t,atol,rtol,res,%DAEOPTIONS,hd);
-                    elseif (lhs==1),
-                        [y]=dassl(x0,t0,t,atol,rtol,res,%DAEOPTIONS,hd);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                if type(varargin(7)) == 1 then
+                    [x0, t0, t, rtol, atol, res, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, hd] = dassl(x0, t0, t, atol, rtol, res, %DAEOPTIONS, hd);
+                    elseif lhs == 1 then
+                        [y] = dassl(x0, t0, t, atol, rtol, res, %DAEOPTIONS, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                     end
 
-                else,
-                    [x0,t0,t,rtol,atol,res,jac]=varargin(:)
-                    if (lhs==2),
-                        [y,hd]=dassl(x0,t0,t,atol,rtol,res,jac,%DAEOPTIONS);
-                    elseif (lhs==1),
-                        [y]=dassl(x0,t0,t,atol,rtol,res,jac,%DAEOPTIONS);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+                else
+                    [x0, t0, t, rtol, atol, res, jac] = varargin(:)
+                    if lhs == 2 then
+                        [y, hd] = dassl(x0, t0, t, atol, rtol, res, jac, %DAEOPTIONS);
+                    elseif lhs == 1 then
+                        [y] = dassl(x0, t0, t, atol, rtol, res, jac, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                     end
 
                 end
 
-            else,
-                [x0,t0,t,atol,res,jac,hd]=varargin(:)
-                if (lhs==2),
-                    [y,hd]=dassl(x0,t0,t,atol,res,jac,%DAEOPTIONS,hd);
-                elseif (lhs==1),
-                    [y]=dassl(x0,t0,t,atol,res,jac,%DAEOPTIONS,hd);
-                else,
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+            else
+                [x0, t0, t, atol, res, jac, hd] = varargin(:)
+                if lhs == 2 then
+                    [y, hd] = dassl(x0, t0, t, atol, res, jac, %DAEOPTIONS, hd);
+                elseif lhs == 1 then
+                    [y] = dassl(x0, t0, t, atol, res, jac, %DAEOPTIONS, hd);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
                 end
             end
 
-        else,
-            error(sprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dae",4,8));
+        else
+            error(sprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dae", 4, 8));
         end
 
-        if (%DAEOPTIONS(2)==0),
-            [r,c]=size(y);
-            y=y([2:1:r],:);
+        if %DAEOPTIONS(2) == 0 then
+            [r, c] = size(y);
+            y = y([2:1:r], :);
         end
-        if (lhs==2),
-            varargout=list(y,hd);
-        elseif (lhs==1),
-            varargout=list(y);
-        else,
-            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",1,2));
+        if lhs == 2 then
+            varargout = list(y, hd);
+        elseif lhs == 1 then
+            varargout = list(y);
+        else
+            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 1, 2));
         end
 
-    elseif (varargin(1)=="root"), // Case root (dasrt)
-        [lhs,rhs]=argn();
-
-        if (rhs==7),// Call without optional arguments
-            [typ,x0,t0,t,res,ng,surface]=varargin(:)
-            if (lhs==2),
-                [y,nn]=dasrt(x0,t0,t,res,ng,surface,%DAEOPTIONS);
-            elseif (lhs==3),
-                [y,nn,hd]=dasrt(x0,t0,t,res,ng,surface,%DAEOPTIONS);
-            else,
-                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+    elseif varargin(1) == "root" then // Case root (dasrt)
+        [lhs, rhs] = argn();
+
+        if rhs == 7 then // Call without optional arguments
+            [typ, x0, t0, t, res, ng, surface] = varargin(:)
+            if lhs == 2 then
+                [y, nn] = dasrt(x0, t0, t, res, ng, surface, %DAEOPTIONS);
+            elseif lhs == 3 then
+                [y, nn, hd] = dasrt(x0, t0, t, res, ng, surface, %DAEOPTIONS);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
             end
 
-        elseif (rhs==11),// Call with all the optional arguments
-            [typ,x0,t0,t,rtol,atol,res,jac,ng,surface,hd]=varargin(:)
-            if (lhs==2),
-                [y,nn]=dasrt(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS,hd);
-            elseif (lhs==3),
-                [y,nn,hd]=dasrt(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS,hd);
-            else,
-                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+        elseif rhs == 11 then // Call with all the optional arguments
+            [typ, x0, t0, t, rtol, atol, res, jac, ng, surface, hd] = varargin(:)
+            if lhs == 2 then
+                [y, nn] = dasrt(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS,hd);
+            elseif lhs == 3 then
+                [y, nn, hd] = dasrt(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, hd);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
             end
 
-        elseif (rhs==8),
+        elseif rhs == 8 then
 
-            if (type(varargin(5))==1),
-                [typ,x0,t0,t,atol,res,ng,surface]=varargin(:)
-                if (lhs==2),
-                    [y,nn]=dasrt(x0,t0,t,atol,res,ng,surface,%DAEOPTIONS);
-                elseif (lhs==3),
-                    [y,nn,hd]=dasrt(x0,t0,t,atol,res,ng,surface,%DAEOPTIONS);
-                else,
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+            if type(varargin(5)) == 1 then
+                [typ, x0, t0, t, atol, res, ng, surface] = varargin(:)
+                if lhs == 2 then
+                    [y, nn] = dasrt(x0, t0, t, atol, res, ng, surface,%DAEOPTIONS);
+                elseif lhs == 3 then
+                    [y, nn, hd] = dasrt(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                 end
 
-            else,
+            else
 
-                if (type(varargin(8))==1),
-                    [typ,x0,t0,t,res,ng,surface,hd]=varargin(:)
-                    if (lhs==2),
-                        [y,nn]=dasrt(x0,t0,t,res,ng,surface,%DAEOPTIONS,hd);
-                    elseif (lhs==3),
-                        [y,nn,hd]=dasrt(x0,t0,t,res,ng,surface,%DAEOPTIONS,hd);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                if type(varargin(8)) == 1 then
+                    [typ, x0, t0, t, res, ng, surface, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = dasrt(x0, t0, t, res, ng, surface, %DAEOPTIONS, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = dasrt(x0, t0, t, res, ng, surface, %DAEOPTIONS, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                     end
 
-                else,
-                    [typ,x0,t0,t,res,jac,ng,surface]=varargin(:)
-                    if (lhs==2),
-                        [y,nn]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS);
-                    elseif (lhs==3),
-                        [y,nn,hd]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS);
-                    else,
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                else
+                    [typ, x0, t0, t, res, jac, ng, surface] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                     end
 
                 end
 
             end
 
-        elseif (rhs==9),
+        elseif rhs == 9 then
 
-            if (type(varargin(5))==1),
+            if type(varargin(5)) == 1 then
 
-                if (type(varargin(6))==1),
-                    [typ,x0,t0,t,rtol,atol,res,ng,surface]=varargin(:)
-                    if (lhs==2),
-                        [y,nn]=dasrt(x0,t0,t,atol,rtol,res,ng,surface,%DAEOPTIONS);
-                    elseif (lhs==3),
-                        [y,nn,hd]=dasrt(x0,t0,t,atol,rtol,res,ng,surface,%DAEOPTIONS);
+                if type(varargin(6)) == 1 then
+                    [typ, x0, t0, t, rtol, atol, res, ng, surface] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = dasrt(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = dasrt(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS);
                     else
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                     end
 
-                else,
+                else
 
-                    if (type(varargin(9))==1),
-                        [typ,x0,t0,t,atol,res,ng,surface,hd]=varargin(:)
-                        if (lhs==2),
-                            [y,nn]=dasrt(x0,t0,t,atol,res,ng,surface,%DAEOPTIONS,hd);
-                        elseif (lhs==3),
-                            [y,nn,hd]=dasrt(x0,t0,t,atol,res,ng,surface,%DAEOPTIONS,hd);
+                    if type(varargin(9)) == 1 then
+                        [typ, x0, t0, t, atol, res, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = dasrt(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = dasrt(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, hd);
                         else
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                         end
 
-                    else,
-                        [typ,x0,t0,t,atol,res,jac,ng,surface]=varargin(:)
-                        if (lhs==2),
-                            [y,nn]=dasrt(x0,t0,t,atol,res,jac,ng,surface,%DAEOPTIONS);
-                        elseif (lhs==3),
-                            [y,nn,hd]=dasrt(x0,t0,t,atol,res,jac,ng,surface,%DAEOPTIONS);
+                    else
+                        [typ, x0, t0, t, atol, res, jac, ng, surface] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = dasrt(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = dasrt(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS);
                         else
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                         end
 
                     end
                 end
 
-            else,
-                [typ,x0,t0,t,res,jac,ng,surface,hd]=varargin(:)
-                if (lhs==2),
-                    [y,nn]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS,hd);
-                elseif (lhs==3),
-                    [y,nn,hd]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS,hd);
+            else
+                [typ, x0, t0, t, res, jac, ng, surface, hd] = varargin(:)
+                if lhs == 2 then
+                    [y, nn] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
+                elseif lhs == 3 then
+                    [y, nn, hd] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
                 else
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                 end
             end
 
-        elseif (rhs==10),
+        elseif rhs == 10 then
 
-            if (type(varargin(5))==1),
+            if type(varargin(5)) == 1 then
 
-                if (type(varargin(6))==1),
+                if type(varargin(6)) == 1 then
 
-                    if (type(varargin(10))==1),
-                        [typ,x0,t0,t,rtol,atol,res,ng,surface,hd]=varargin(:)
-                        if (lhs==2),
-                            [y,nn]=dasrt(x0,t0,t,atol,rtol,res,ng,surface,%DAEOPTIONS,hd);
-                        elseif (lhs==3),
-                            [y,nn,hd]=dasrt(x0,t0,t,atol,rtol,res,ng,surface,%DAEOPTIONS,hd);
+                    if type(varargin(10)) == 1 then
+                        [typ, x0, t0, t, rtol, atol, res, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = dasrt(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = dasrt(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, hd);
                         else
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                         end
 
-                    else,
-                        [typ,x0,t0,t,rtol,atol,res,jac,ng,surface]=varargin(:)
-                        if (lhs==2),
-                            [y,nn]=dasrt(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS);
-                        elseif (lhs==3),
-                            [y,nn,hd]=dasrt(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS);
+                    else
+                        [typ, x0, t0, t, rtol, atol, res, jac, ng, surface] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = dasrt(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = dasrt(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS);
                         else
-                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                         end
                     end
 
-                else,
-                    [typ,x0,t0,t,atol,res,jac,ng,surface,hd]=varargin(:)
-                    if (lhs==2),
-                        [y,nn]=dasrt(x0,t0,t,atol,res,jac,ng,surface,%DAEOPTIONS,hd);
-                    elseif (lhs==3),
-                        [y,nn,hd]=dasrt(x0,t0,t,atol,res,jac,ng,surface,%DAEOPTIONS,hd);
+                else
+                    [typ, x0, t0, t, atol, res, jac, ng, surface, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = dasrt(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = dasrt(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, hd);
                     else
-                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                     end
 
                 end
 
-            else,
-                [typ,x0,t0,t,res,jac,ng,surface,hd]=varargin(:)
-                if (lhs==2),
-                    [y,nn]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS,hd);
-                elseif (lhs==3),
-                    [y,nn,hd]=dasrt(x0,t0,t,res,jac,ng,surface,%DAEOPTIONS,hd);
+            else
+                [typ, x0, t0, t, res, jac, ng, surface, hd] = varargin(:)
+                if lhs == 2 then
+                    [y, nn] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
+                elseif lhs == 3 then
+                    [y, nn, hd] = dasrt(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
                 else
-                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
                 end
             end
 
-        else,
-            error(sprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dae",7,11));
+        else
+            error(sprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dae", 7, 11));
+        end
+        if %DAEOPTIONS(2) == 0 then
+            [r, c] = size(y);
+            y = y([2:1:r], :);
         end
+        if lhs == 2 then
+            varargout = list(y, nn);
+        elseif lhs == 3 then
+            varargout = list(y, nn, hd);
+        else
+            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+        end
+
+    elseif varargin(1) == "root2" then // Case root2 (daskr)
+        [lhs, rhs] = argn();
+
+        if rhs == 7 then // Call without optional arguments
+            [typ, x0, t0, t, res, ng, surface] = varargin(:)
+            if lhs == 2 then
+                [y, nn] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS);
+            elseif lhs == 3 then
+                [y, nn, hd] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+            end
+
+        elseif rhs == 13 then // Call with all the optional arguments
+            [typ, x0, t0, t, rtol, atol, res, jac, ng, surface, psol, pjac, hd] = varargin(:)
+            if lhs == 2 then
+                [y, nn] = daskr(x0,t0,t,atol,rtol,res,jac,ng,surface,%DAEOPTIONS, psol, pjac, hd);
+            elseif lhs == 3 then
+                [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+            end
+
+        elseif rhs == 8 then
+
+            if type(varargin(5)) == 1 then
+                [typ, x0, t0, t, atol, res, ng, surface] = varargin(:)
+                if lhs == 2 then
+                    [y, nn] = daskr(x0, t0, t, atol, res, ng, surface,%DAEOPTIONS);
+                elseif lhs == 3 then
+                    [y, nn, hd] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                end
+
+            else
+
+                if type(varargin(8)) == 1 then
+                    [typ, x0, t0, t, res, ng, surface, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+
+                else
+                    [typ, x0, t0, t, res, jac, ng, surface] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+
+                end
+
+            end
+
+        elseif rhs == 9 then
+
+            if type(varargin(5)) == 1 then
+
+                if type(varargin(6)) == 1 then
+                    [typ, x0, t0, t, rtol, atol, res, ng, surface] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+
+                else
+
+                    if type(varargin(9)) == 1 then
+                        [typ, x0, t0, t, atol, res, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+
+                    else
+                        [typ, x0, t0, t, atol, res, jac, ng, surface] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+
+                    end
+                end
+
+            else
+                if type(varargin(9)) == 1 then
+                    [typ, x0, t0, t, res, jac, ng, surface, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+                else
+                    [typ, x0, t0, t, res, ng, surface, psol, pjac] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+                end
+
+            end
+
+        elseif rhs == 10 then
 
-        if (%DAEOPTIONS(2)==0),
-            [r,c]=size(y);
-            y=y([2:1:r],:);
+            if type(varargin(5)) == 1 then
+
+                if type(varargin(6)) == 1 then
+
+                    if type(varargin(10)) == 1 then
+                        [typ, x0, t0, t, rtol, atol, res, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+
+                    else
+                        [typ, x0, t0, t, rtol, atol, res, jac, ng, surface] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    end
+
+                else
+                    if type(varargin(10)) == 1 then
+                        [typ, x0, t0, t, atol, res, jac, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    else
+                        [typ, x0, t0, t, atol, res, ng, surface, psol, pjac] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    end
+
+                end
+
+            else
+                if type(varargin(10)) == 1 then
+                    [typ, x0, t0, t, res, ng, surface, psol, pjac, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+                else
+                    [typ, x0, t0, t, res, jac, ng, surface, psol, pjac] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+                end
+            end
+
+        elseif rhs == 11 then
+
+            if type(varargin(5)) == 1 then
+
+                if type(varargin(6)) == 1 then
+
+                    if type(varargin(11)) == 1 then
+                        [typ, x0, t0, t, rtol, atol, res, jac, ng, surface, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+
+                    else
+                        [typ, x0, t0, t, rtol, atol, res, ng, surface, psol, pjac] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, psol, pjac);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    end
+
+                else
+                    if type(varargin(11)) == 1 then
+                        [typ, x0, t0, t, atol, res, ng, surface, psol, pjac, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    else
+                        [typ, x0, t0, t, atol, res, jac, ng, surface, psol, pjac] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    end
+                end
+
+            else
+                [typ, x0, t0, t, res, jac, ng, surface, psol, pjac, hd] = varargin(:)
+                if lhs == 2 then
+                    [y, nn] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                elseif lhs == 3 then
+                    [y, nn, hd] = daskr(x0, t0, t, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                else
+                    error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                end
+            end
+
+        elseif rhs == 12 then
+
+            if type(varargin(5)) == 1 then
+
+                if type(varargin(6)) == 1 then
+
+                    if type(varargin(12)) == 1 then
+                        [typ, x0, t0, t, rtol, atol, res, ng, surface, psol, pjac, hd] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+
+                    else
+                        [typ, x0, t0, t, rtol, atol, res, jac, ng, surface, psol, pjac] = varargin(:)
+                        if lhs == 2 then
+                            [y, nn] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                        elseif lhs == 3 then
+                            [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac);
+                        else
+                            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                        end
+                    end
+
+                else
+                    [typ, x0, t0, t, atol, res, jac, ng, surface, psol, pjac, hd] = varargin(:)
+                    if lhs == 2 then
+                        [y, nn] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                    elseif lhs == 3 then
+                        [y, nn, hd] = daskr(x0, t0, t, atol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+                    else
+                        error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+                    end
+                end
+            end
+
+        elseif rhs == 13 then
+
+            [typ, x0, t0, t, rtol, atol, res, jac, ng, surface, psol, pjac, hd] = varargin(:)
+            if lhs == 2 then
+                [y, nn] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+            elseif lhs == 3 then
+                [y, nn, hd] = daskr(x0, t0, t, atol, rtol, res, jac, ng, surface, %DAEOPTIONS, psol, pjac, hd);
+            else
+                error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
+            end
+
+        else
+            error(sprintf(gettext("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dae", 7, 13));
+        end
+        if %DAEOPTIONS(2) == 0 then
+            [r, c] = size(y);
+            y = y([2:1:r], :);
         end
-        if (lhs==2),
-            varargout=list(y,nn);
-        elseif (lhs==3),
-            varargout=list(y,nn,hd);
-        else,
-            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae",2,3));
+        if lhs == 2 then
+            varargout = list(y, nn);
+        elseif lhs == 3 then
+            varargout = list(y, nn, hd);
+        else
+            error(sprintf(gettext("%s: Wrong number of output argument(s): %d or %d expected.\n"), "dae", 2, 3));
         end
     else
-        error(sprintf(gettext("%s: Invalid option %s: real matrix expected.\n"),"dae","root"));
+        error(sprintf(gettext("%s: Invalid option %s: real matrix expected.\n"), "dae", "root"));
     end
 
 endfunction
index 40a4844..9cf3c76 100644 (file)
@@ -1,4 +1,5 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - Scilab Enterprises - Paul Bignier : added daskr options
 // Copyright (C) 2008 - INRIA - Sabine GAUZERE
 // ...
 //
@@ -20,7 +21,7 @@ function [%DAEOPTIONS] = daeoptions(%DAEOPTIONS)
     //options = [1,0,0,%inf,0,2,500,12,5,0,-1,-1];
     options = list ([],0,[],[],[],0,0);
     //default = [string(options(1:10)),sci2exp(options(11:12))]
-    default = ["[]","0","[]","[]","[]","0","0"];
+    default = ["0","[]","0","[]","[]","[]","0","0","0","[]","0","0","[]","[]","1"];
     //default(find(default=='Inf'))="%inf"
 
     if argn(2)>0 then
@@ -28,12 +29,17 @@ function [%DAEOPTIONS] = daeoptions(%DAEOPTIONS)
     end
 
     //lab_=[string(options(1:10)),sci2exp(options(11:12))]
-    lab_= ["[]","0","[]","[]","[]","0","0"];
+    lab_= ["0","[]","0","[]","[]","[]","0","0","0","[]","0","0","[]","[]","1"];
     //lab_(find(lab_=="[]"))="[]"
 
 
     chapeau=["Defining %DAEOPTIONS variable";
     "*****************************";
+    "Meaning of solver:";
+    "solver is a real scalar which selects the solver";
+    "0 : use dassl/dasrt";
+    "1 : use daskr";
+    " ";
     "Meaning of tstop:";
     "tstop is a real scalar which gives the maximum time";
     "for which g is allowed to be evaluated";
@@ -55,28 +61,89 @@ function [%DAEOPTIONS] = daeoptions(%DAEOPTIONS)
     " ";
     "Meaning of isest:";
     "0 : if g(t0,y0,ydot0)=0";
-    "1 : if ydot0 is just an estimation";
+    "1 : if ydot0 is just an estimation and you want to use dassl/dasrt";
+    "[+-1,...,+-1]: if ydot0 is just an estimation and you want to use daskr, with:";
+    "1 if y(i) is a differential variable and";
+    "-1 if y(i) is an algebraic variable";
+    "(if its derivatives do not appear explicitly in the function g(t, y, ydot)).";
     " "
+    "The following values are only used by daskr (solver = 1)";
+    " ";
+    "Meaning of method:";
+    "0 : use GMRes Krylov method and provide a psol routine in dae";
+    "1 : use direct method";
+    " ";
+    "Meaning of Kry_params:";
+    "Treat as dummy argument if you have set method=0. Otherwise, set :";
+    "[] : default parameters";
+    "[maxl kmp nrmax epli] : where";
+    "- maxl  = maximum number of iterations in the GMRes algorithm (default min(5, neq)),";
+    "- kmp   = number of vectors on which orthogonalization is done in the GMRes algorithm (default maxl),";
+    "- nrmax = maximum number of restarts of the GMRes algorithm per nonlinear iteration (default 5),";
+    "- epli  = convergence test constant in GMRes algorithm (default 0.05).";
+    " ";
+    "Meaning of init:";
+    "Treat as dummy argument if you have set isest=0. Otherwise, set :";
+    "0 : stop after initial values computation";
+    "1 : proceed to integration";
+    " ";
+    "Meaning of precond:";
+    "Treat as dummy argument if you have set method=0. Otherwise, set :";
+    "0 : specify a specific pjac routine in dae()";
+    "1 : use the default";
+    " ";
+    "Meaning of control:";
+    "[] : if you wish to control errors locally on all the variables then set to [].";
+    "[+-1,...,+-1] :  1 if y(i) is a differential variable and -1 if y(i) is an algebraic variable";
+    "(if its derivatives do not appear explicitly in the function g(t, y, ydot)).";
+    " ";
+    "Meaning of heuristic:";
+    "Treat as dummy argument if you have set isest=0. Otherwise, set :";
+    "[] : default parameters";
+    "[mxnit mxnj mxnh lsoff stptol epinit] : where";
+    "mxnj   = maximum number of Jacobian or preconditioner evaluations (default 6 if info(8)=0, 2 otherwise),";
+    "lsoff  = flag to turn off the linesearch algorithm (lsoff = 0 means linesearch is on, lsoff = 1 means it is turned off) (default 0)";
+    "mxnj   = maximum number of Jacobian or preconditioner evaluations (default 6 if info(8)=0, 2 otherwise),";
+    "stptol = minimum scaled step in linesearch algorithm (default (unit roundoff)^(2/3)),";
+    "epinit = swing factor in the Newton iteration convergence test (default 0.01).";
+    " ";
+    "Meaning of verbosity:";
+    "0 : standard printing";
+    "1 : minimal printing";
+    "2 : full printing";
+    " ";
     "Default values are given in square brackets"
     "If the function is called without argument, default values"+...
     " are used"
     ]
 
-    dims = list("vec",-1,"vec",1,"vec",-1,"vec",-1,"vec",-1,..
-    "vec",1,"vec",1);
+    dims = list("vec",1,"vec",-1,"vec",1,"vec",-1,"vec",-1,"vec",-1,..
+    "vec",1,"vec",-1,"vec",1,"vec",-1,"vec",1,"vec",1,"vec",-1,"vec",-1,"vec",1);
 
 
-    labels = ["tstop (maximum time) ","imode (assumes imode = 0 or 1)",...
+    labels = ["solver (assumes solver = 0 or 1) ",...
+    "tstop (maximum time) ","imode (assumes imode = 0 or 1)",...
     "band ()",...
     "maxstep (max step size)","stepin (initial step size)",...
     "nonneg (assumes nonneg = 0 or 1)",...
-    "isest (assumes isest = 0 or 1)"] +"    ["+default+"]";
+    "isest (assumes isest = 0, 1 or [+-1,...])",...
+    "method (assumes method = 0 or 1)","Kry_params",...
+    "init (assumes init = 0 or 1)","precond (assumes precond = 0 or 1)",...
+    "control","heuristic","verbosity (assumes verbosity = 0 or 1)"] +"    ["+default+"]";
+
 
-    [ok,tstop,imode,band,maxstep,stepin,nonneg,isest] = getvalue(chapeau,labels,dims,lab_);
+
+    [solver,ok,tstop,imode,band,maxstep,stepin,nonneg,isest,...
+    method,Kry_params,init,precond,control,heuristic,verbosity] = getvalue(chapeau,labels,dims,lab_);
     //ml = mlmu(1);
     //mu = mlmu(2);
 
-    DAEOPTIONS = list(tstop,imode,band,maxstep,stepin,nonneg,isest);
+    if solver == 0 then
+        DAEOPTIONS = list(tstop,imode,band,maxstep,stepin,nonneg,isest);
+    else
+        DAEOPTIONS = list(tstop,imode,band,maxstep,stepin,nonneg,isest,...
+        method,Kry_params,init,precond,control,heuristic,verbosity);
+    end
 
     if DAEOPTIONS<>list() then
         %DAEOPTIONS=DAEOPTIONS
index 0cc7553..15ddc1b 100644 (file)
@@ -282,7 +282,8 @@ c     --   subvariable consistent(info) --
       il10e7=iadr(l10+istk(il10+1+7)-1)
       m10e7 =istk(il10e7+2)*istk(il10e7+2)
       l10e7 = sadr(il10e7+4)
-      if(m10e7.eq.0) then
+      if(m10e7.eq.0.or.(m10e7.eq.1.and.stk(l10e7).eq.0)) then
+c        info is then [] or [0]
          info(11)=0
       else
 c        info then looks like list(..., [+-1 +-1 ... +-1 +-1],...)
@@ -860,4 +861,3 @@ c     Remise en place de la pile
  1150 call unsfdcopy(lw-lw0,stk(lw0),1,stk(l0),1)
       return
       end
-
index 535494b..a593f9c 100644 (file)
@@ -1,5 +1,6 @@
 // ===========================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
 // Copyright (C) 2008 - INRIA - Sabinere Gauzere
 //
 //  This file is distributed under the same license as the Scilab package.
@@ -92,7 +93,7 @@ y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
 t0=0;t=[0.01,0.1,1,10,100];
 rtol=0;atol=1.d-6;
 y=dae([y0,y0d],t0,t,rtol,atol,'dres2');
-//                 Root finder
+//                 Root finder (dasrt)
 //
 //-----------------------------------------------------------------------
 // First problem.
@@ -300,3 +301,162 @@ y0=ones(n,1);yd0=0*y0;
 yb=dae([y0,yd0],0,0:0.1:10,'myres','myjac');
 a = norm(y-yb);
 if (a > %eps * 1e5) then bugmes();quit;end
+//                 Root finder (daskr)
+//
+y0=1;t=2:6;t0=1;y0d=3;
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1','psol1','pjac1');
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1','psol1','pjac1');
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+%DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1');
+assert_checkalmostequal(nn(1),2.53,0.001);
+// 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')
+Warning : redefining function: res1                    . Use funcprot(0) to avoid this message
+
+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;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.53,0.001);
+// Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
+// pjac uses the macro res1 defined above.
+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
+            ires = -1;
+            return;
+        end
+        delinv = 1/del;
+        for j=1:neq
+            wp(nrow+j) = delinv*(e(j)-savr(j));
+            if isnan(wp(nrow+j)) then
+                ires = -1;
+                return;
+            end
+            iwp(nrow+j, 1) = i;
+            iwp(nrow+j, 2) = j;
+        end
+        nrow = nrow + neq;
+        y(i) = ysave;
+        ydot(i) = ypsave;
+    end
+endfunction
+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
+y0=1;t=2:6;t0=1;y0d=3;
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,'gr1',psol,pjac);
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,'gr1',psol,pjac);
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,'gr1',psol,pjac);
+assert_checkalmostequal(nn(1),2.53,0.001);
+//C
+//C-----------------------------------------------------------------------
+//C Second problem (Van Der Pol oscillator).
+//C The initial value problem is..
+//C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
+//C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
+//C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
+//C The root function is  G = Y1.
+//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-----------------------------------------------------------------------
+%DAEOPTIONS=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]=dae("root2",[y0,y0d],t0,t,rtol,atol,'res2','jac2',ng,'gr2');
+assert_checkalmostequal(nn(1),81.163512,0.001);
+deff('[delta,ires]=res2(t,y,ydot)',...
+'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,'jac2',ng,'gr2');
+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]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,'gr2');
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2);
+// Same problem, with psol and pjac example routines
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,'gr2','psol1','pjac1');
+assert_checkalmostequal(nn(1),81.163512,0.009);
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,'psol1','pjac1');
+assert_checkalmostequal(nn(1),81.163512,0.009);
+// 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]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,'gr2',psol,pjac);
+assert_checkalmostequal(nn(1),81.163512,0.003);
+deff('s=gr2(t,y,yd)','s=y(1)')
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,psol,pjac);
+assert_checkalmostequal(nn(1),81.163512,0.003);
+%DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+//           Hot Restart
+[yy,nn,hotd]=dae("root2",[y0,y0d],t0,t,rtol,atol,'res2','jac2',ng,'gr2');
+t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
+[yy,nn,hotd]=dae("root2",[y01,y0d1],t01,t,rtol,atol,'res2','jac2',ng,'gr2',hotd);
+assert_checkalmostequal(nn(1),162.57763,0.004);
index 0095c1e..a89d494 100644 (file)
@@ -1,5 +1,6 @@
 // ===========================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier : added "root2" (daskr)
 // Copyright (C) 2008 - INRIA - Sabinere Gauzere
 //
 //  This file is distributed under the same license as the Scilab package.
@@ -20,34 +21,34 @@ ilib_verbose(0);
 t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
 info=list([],0,[],[],[],0,0);
 //    Calling Scilab functions
-deff('[r,ires]=dres1(t,y,ydot)','r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0')
-deff('pd=djac1(t,y,ydot,cj)','pd=[cj+10,0;1,1]')
+deff("[r,ires]=dres1(t,y,ydot)","r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0")
+deff("pd=djac1(t,y,ydot,cj)","pd=[cj+10,0;1,1]")
 //   scilab function, without jacobian
 yy0=dae([y0,y0d],t0,t,dres1);
 //   scilab functions, with jacobian
 yy1=dae([y0,y0d],t0,t,dres1,djac1);
 // fortran routine, without jacobian
-yy2=dae([y0,y0d],t0,t,'dres1');   //=yy0
+yy2=dae([y0,y0d],t0,t,"dres1");   //=yy0
 if norm(yy2-yy0,1)>1E-5 then pause,end
 // fortran routines, with jacobian
-yy3=dae([y0,y0d],t0,t,'dres1','djac1');  //=yy1
+yy3=dae([y0,y0d],t0,t,"dres1","djac1");  //=yy1
 if norm(yy3-yy1,1)>1E-5 then pause,end
-yy3bis=dae([y0,y0d],t0,t,'dres1',djac1);
+yy3bis=dae([y0,y0d],t0,t,"dres1",djac1);
 // call fortran dres1 and scilab's djac1
-yy3ter=dae([y0,y0d],t0,t,dres1,'djac1');
+yy3ter=dae([y0,y0d],t0,t,dres1,"djac1");
 //
 // with specific atol and rtol parameters
 atol=1.d-6;rtol=0;
 yy4=dae([y0,y0d],t0,t,rtol,atol,dres1);
-yy5=dae([y0,y0d],t0,t,rtol,atol,'dres1'); //=yy4
+yy5=dae([y0,y0d],t0,t,rtol,atol,"dres1"); //=yy4
 if norm(yy5-yy4,1)>1E-9 then pause,end
 yy6=dae([y0,y0d],t0,t,rtol,atol,dres1,djac1);
-yy7=dae([y0,y0d],t0,t,rtol,atol,'dres1','djac1'); //==yy6
+yy7=dae([y0,y0d],t0,t,rtol,atol,"dres1","djac1"); //==yy6
 if norm(yy7-yy6,1)>1E-12 then pause,end
 //
 //   Testing E xdot - A x=0
 //   x(0)=x0;   xdot(0)=xd0
-rand('seed',0);
+rand("seed",0);
 nx=5;
 E=rand(nx,1)*rand(1,nx);
 A=rand(nx,nx);
@@ -57,7 +58,7 @@ pp=Si*E;
 [q,m]=fullrf(pp);
 x0=q(:,1);
 x0d=pinv(E)*A*x0;
-deff('[r,ires]=g(t,x,xdot)','r=E*xdot-A*x;ires=0');
+deff("[r,ires]=g(t,x,xdot)","r=E*xdot-A*x;ires=0");
 t=[1,2,3];t0=0;
 %DAEOPTIONS=list([],0,[],[],[],0,0);
 x=dae([x0,x0d],t0,t,g);
@@ -70,10 +71,10 @@ t=1.5409711;
 ww=dae([x0,x0d],t0,t,g);
 ww=[t;ww];
 if abs(ww(5)-1)>0.001 then pause,end
-deff('[rt]=surface(t,y,yd)','rt=y(4)-1');
+deff("[rt]=surface(t,y,yd)","rt=y(4)-1");
 nbsurf=1;
 [yyy,nnn]=dae("root",[x0,x0d],t0,t,g,nbsurf,surface);
-deff('pd=j(t,y,ydot,cj)','pd=cj*E-A');
+deff("pd=j(t,y,ydot,cj)","pd=cj*E-A");
 x=dae([x0,x0d],t0,t,g,j);
 x=[t;x];
 x2=x(2:nx+1,1);
@@ -94,9 +95,9 @@ delta=0*y0;
 y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
 t0=0;t=[0.01,0.1,1,10,100];
 rtol=0;atol=1.d-6;
-y=dae([y0,y0d],t0,t,rtol,atol,'dres2');
+y=dae([y0,y0d],t0,t,rtol,atol,"dres2");
 
-//                 Root finder
+//                 Root finder (dasrt)
 //
 //-----------------------------------------------------------------------
 // First problem.
@@ -110,16 +111,16 @@ y=dae([y0,y0d],t0,t,rtol,atol,'dres2');
 y0=1;t=2:6;t0=1;y0d=3;
 %DAEOPTIONS=list([],0,[],[],[],0,0);
 atol=1.d-6;rtol=0;ng=2;
-[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1');
+[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 if abs(nn(1)-2.47)>0.001 then pause,end
 y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
-[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1');
+[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 if abs(nn(1)-2.5)>0.001 then pause,end
 y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
-[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,'res1',ng,'gr1');
+[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
 if abs(nn(1)-2.53)>0.001 then pause,end
-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]')
+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;
 %DAEOPTIONS=list([],0,[],[],[],0,0);
 atol=1.d-6;rtol=0;ng=2;
@@ -145,83 +146,81 @@ if abs(nn(1)-2.53)>0.001 then pause,end
 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;
 %DAEOPTIONS=list([],0,[],[],[],0,0);
-[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,'res2','jac2',ng,'gr2');
+[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 if abs(nn(1)-81.163512)>0.001 then pause,end
-deff('[delta,ires]=res2(t,y,ydot)',...
-'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
-[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,'jac2',ng,'gr2');
-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]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,'gr2');
-deff('s=gr2(t,y,yd)','s=y(1)')
+deff("[delta,ires]=res2(t,y,ydot)",...
+"ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")
+[yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,"jac2",ng,"gr2");
+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]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2");
+deff("s=gr2(t,y,yd)","s=y(1)")
 [yy,nn]=dae("root",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2);
 //           Hot Restart
-[yy,nn,hotd]=dae("root",[y0,y0d],t0,t,rtol,atol,'res2','jac2',ng,'gr2');
+[yy,nn,hotd]=dae("root",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(1:2,qq);y0d1=yy(2:3,qq);
 %DAEOPTIONS=list([],0,[],[],[],0,0);
-[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,rtol,atol,'res2','jac2',ng,'gr2',hotd);
+[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,rtol,atol,"res2","jac2",ng,"gr2",hotd);
 if abs(nn(1)-162.57763)>0.004 then pause,end
 
 //same with C code
-
 ilib_verbose(0);
 
 cd TMPDIR;
 mkdir("dae_test1");
 cd("dae_test1");
 
-code=['#include <math.h>'
-      'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
-      '{res[0] = yd[0] - y[1];'
-      ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
-      ' '
-      'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
-      '{pd[0]=*cj - 0.0;'
-      ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
-      ' pd[2]=    - 1.0;'
-      ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
-      ' '
-      'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
-      '{ groot[0] = y[0];}'];
-mputl(code,'t22.c') ;
-ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c','','c');
-exec('loader.sce');
+code=["#include <math.h>"
+"void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)"
+"{res[0] = yd[0] - y[1];"
+" res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}"
+" "
+"void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)"
+"{pd[0]=*cj - 0.0;"
+" pd[1]=    - (-200.0*y[0]*y[1] - 1.0);"
+" pd[2]=    - 1.0;"
+" pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}"
+" "
+"void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)"
+"{ groot[0] = y[0];}"];
+mputl(code,"t22.c") ;
+ilib_for_link(["res22" "jac22" "gr22"],"t22.c","","c");
+exec("loader.sce");
 
 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;
 %DAEOPTIONS =list([],0,[],[],[],0,0);
- //hot restart
+//hot restart
 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
-[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',hotd);
+[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",hotd);
 
 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;
 %DAEOPTIONS =list([],0,[],[],[],0,0);
- [yy,nn]=dae("root",[y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22');
+[yy,nn]=dae("root",[y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22");
 //hot restart
-[yy,nn,hotd]=dae("root",[y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22');
+[yy,nn,hotd]=dae("root",[y0,y0d],t0,t,atol,rtol,"res22","jac22",ng,"gr22");
 t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
-[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',hotd);
+[yy,nn,hotd]=dae("root",[y01,y0d1],t01,t,atol,rtol,"res22","jac22",ng,"gr22",hotd);
 
 //banded systems
-
 A=[-17,6,3,0,0,0,0,0,0,0;
-   8,-12,9,4,0,0,0,0,0,0;
-   0,7,-17,3,8,0,0,0,0,0;
-   0,0,3,-13,2,2,0,0,0,0;
-   0,0,0,4,-18,6,4,0,0,0;
-   0,0,0,0,7,-13,7,0,0,0;
-   0,0,0,0,0,7,-16,7,9,0;
-   0,0,0,0,0,0,5,-17,8,1;
-   0,0,0,0,0,0,0,4,-14,8;
-   0,0,0,0,0,0,0,0,10,-10];
+8,-12,9,4,0,0,0,0,0,0;
+0,7,-17,3,8,0,0,0,0,0;
+0,0,3,-13,2,2,0,0,0,0;
+0,0,0,4,-18,6,4,0,0,0;
+0,0,0,0,7,-13,7,0,0,0;
+0,0,0,0,0,7,-16,7,9,0;
+0,0,0,0,0,0,5,-17,8,1;
+0,0,0,0,0,0,0,4,-14,8;
+0,0,0,0,0,0,0,0,10,-10];
 n=size(A,1);
 //Full jacobian case for reference
 function [r,ires]=res(t,y,yd)
-  r=yd-A*y; ires=0
+    r=yd-A*y; ires=0
 endfunction
 function pd=jac(x,y,yd,cj)
-  pd=A+cj*eye()
+    pd=A+cj*eye()
 endfunction
 
 y0=ones(n,1);yd0=0*y0;
@@ -240,78 +239,240 @@ norm(y1-yb1);
 cd TMPDIR;
 mkdir("dae_test2");
 cd("dae_test2");
-ccode=['#include <math.h>'
-       'void myres(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
-       '{'
-       '  *ires =0;'
-       '  res[0]=yd[0]+17.0*y[0]- 6.0*y[1]- 3.0*y[2];'
-       '  res[1]=yd[1]-8.0*y[0]+12.0*y[1]- 9.0*y[2]- 4.0*y[3];'
-       '  res[2]=yd[2]          -7.0*y[1]+17.0*y[2]- 3.0*y[3]- 8.0*y[4];'
-       '  res[3]=yd[3]                    -3.0*y[2]+13.0*y[3]- 2.0*y[4]- 2.0*y[5];'
-       '  res[4]=yd[4]                              -4.0*y[3]+18.0*y[4]- 6.0*y[5]- 4.0*y[6];'
-       '  res[5]=yd[5]                                        -7.0*y[4]+13.0*y[5]- 7.0*y[6]- 0.0*y[7];'
-       '  res[6]=yd[6]                                                  -7.0*y[5]+16.0*y[6]- 7.0*y[7]- 9.0*y[8];'
-       '  res[7]=yd[7]                                                            -5.0*y[6]+17.0*y[7]- 8.0*y[8]- 1.0*y[9];'
-       '  res[8]=yd[8]                                                                      -4.0*y[7]+14.0*y[8]- 8.0*y[9];'
-       '  res[9]=yd[9]                                                                               -10.0*y[8]+10.0*y[9];'
-       '}'
-       'void myjac(double *t,double *y,double *yd,double *res,double *cj,double *rpar,int *ipar)'
-       '{'
-       '  res[0]=0.0;'
-       '  res[1]=0.0;'
-       '  res[2]=0.0;'
-       '  res[3]=-17.0+*cj;'
-       '  res[4]=8.0;'
-       '  res[5]=0.0;'
-       '  res[6]=0.0;'
-       '  res[7]=6.0;'
-       '  res[8]=-12.0+*cj;'
-       '  res[9]=7.0;'
-       '  res[10]=0.0;'
-       '  res[11]=3.0;'
-       '  res[12]=9.0;'
-       '  res[13]=-17.0+*cj;'
-       '  res[14]=3.0;'
-       '  res[15]=0.0;'
-       '  res[16]=4.0;'
-       '  res[17]=3.0;'
-       '  res[18]=-13.0+*cj;'
-       '  res[19]=4.0;'
-       '  res[20]=0.0;'
-       '  res[21]=8.0;'
-       '  res[22]=2.0;'
-       '  res[23]=-18.0+*cj;'
-       '  res[24]=7.0;'
-       '  res[25]=0.0;'
-       '  res[26]=2.0;'
-       '  res[27]=6.0;'
-       '  res[28]=-13.0+*cj;'
-       '  res[29]=7.0;'
-       '  res[30]=0.0;'
-       '  res[31]=4.0;'
-       '  res[32]=7.0;'
-       '  res[33]=-16.0+*cj;'
-       '  res[34]=5.0;'
-       '  res[35]=0.0;'
-       '  res[36]=0.0;'
-       '  res[37]=7.0;'
-       '  res[38]=-17.0+*cj;'
-       '  res[39]=4.0;'
-       '  res[40]=0.0;'
-       '  res[41]=9.0;'
-       '  res[42]=8.0;'
-       '  res[43]=-14.0+*cj;'
-       '  res[44]=10.0;'
-       '  res[45]=0.0;'
-       '  res[46]=1.0;'
-       '  res[47]=8.0;'
-       '  res[48]=-10.0+*cj;'
-       '  res[49]=0.0;'
-       '}'];
-mputl(ccode,'band.c'); //create the C file of myjac
-ilib_for_link(['myres','myjac'],'band.c','','c');//compile
-exec('loader.sce');
+ccode=["#include <math.h>"
+"void myres(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)"
+"{"
+"  *ires =0;"
+"  res[0]=yd[0]+17.0*y[0]- 6.0*y[1]- 3.0*y[2];"
+"  res[1]=yd[1]-8.0*y[0]+12.0*y[1]- 9.0*y[2]- 4.0*y[3];"
+"  res[2]=yd[2]          -7.0*y[1]+17.0*y[2]- 3.0*y[3]- 8.0*y[4];"
+"  res[3]=yd[3]                    -3.0*y[2]+13.0*y[3]- 2.0*y[4]- 2.0*y[5];"
+"  res[4]=yd[4]                              -4.0*y[3]+18.0*y[4]- 6.0*y[5]- 4.0*y[6];"
+"  res[5]=yd[5]                                        -7.0*y[4]+13.0*y[5]- 7.0*y[6]- 0.0*y[7];"
+"  res[6]=yd[6]                                                  -7.0*y[5]+16.0*y[6]- 7.0*y[7]- 9.0*y[8];"
+"  res[7]=yd[7]                                                            -5.0*y[6]+17.0*y[7]- 8.0*y[8]- 1.0*y[9];"
+"  res[8]=yd[8]                                                                      -4.0*y[7]+14.0*y[8]- 8.0*y[9];"
+"  res[9]=yd[9]                                                                               -10.0*y[8]+10.0*y[9];"
+"}"
+"void myjac(double *t,double *y,double *yd,double *res,double *cj,double *rpar,int *ipar)"
+"{"
+"  res[0]=0.0;"
+"  res[1]=0.0;"
+"  res[2]=0.0;"
+"  res[3]=-17.0+*cj;"
+"  res[4]=8.0;"
+"  res[5]=0.0;"
+"  res[6]=0.0;"
+"  res[7]=6.0;"
+"  res[8]=-12.0+*cj;"
+"  res[9]=7.0;"
+"  res[10]=0.0;"
+"  res[11]=3.0;"
+"  res[12]=9.0;"
+"  res[13]=-17.0+*cj;"
+"  res[14]=3.0;"
+"  res[15]=0.0;"
+"  res[16]=4.0;"
+"  res[17]=3.0;"
+"  res[18]=-13.0+*cj;"
+"  res[19]=4.0;"
+"  res[20]=0.0;"
+"  res[21]=8.0;"
+"  res[22]=2.0;"
+"  res[23]=-18.0+*cj;"
+"  res[24]=7.0;"
+"  res[25]=0.0;"
+"  res[26]=2.0;"
+"  res[27]=6.0;"
+"  res[28]=-13.0+*cj;"
+"  res[29]=7.0;"
+"  res[30]=0.0;"
+"  res[31]=4.0;"
+"  res[32]=7.0;"
+"  res[33]=-16.0+*cj;"
+"  res[34]=5.0;"
+"  res[35]=0.0;"
+"  res[36]=0.0;"
+"  res[37]=7.0;"
+"  res[38]=-17.0+*cj;"
+"  res[39]=4.0;"
+"  res[40]=0.0;"
+"  res[41]=9.0;"
+"  res[42]=8.0;"
+"  res[43]=-14.0+*cj;"
+"  res[44]=10.0;"
+"  res[45]=0.0;"
+"  res[46]=1.0;"
+"  res[47]=8.0;"
+"  res[48]=-10.0+*cj;"
+"  res[49]=0.0;"
+"}"];
+mputl(ccode,"band.c"); //create the C file of myjac
+ilib_for_link(["myres","myjac"],"band.c","","c");//compile
+exec("loader.sce");
 y0=ones(n,1);yd0=0*y0;
-yb=dae([y0,yd0],0,0:0.1:10,'myres','myjac');
+yb=dae([y0,yd0],0,0:0.1:10,"myres","myjac");
 a = norm(y-yb);
 if (a > %eps * 1e5) then pause,end
+
+//                 Root finder (daskr)
+//
+y0=1;t=2:6;t0=1;y0d=3;
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1","psol1","pjac1");
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1","psol1","pjac1");
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+%DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res1",ng,"gr1");
+assert_checkalmostequal(nn(1),2.53,0.001);
+
+// 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;
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
+assert_checkalmostequal(nn(1),2.53,0.001);
+
+// Same problem, but using macros for the preconditioner evaluation and application functions 'pjac' and 'psol'
+// pjac uses the macro res1 defined above.
+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
+            ires = -1;
+            return;
+        end
+        delinv = 1/del;
+        for j=1:neq
+            wp(nrow+j) = delinv*(e(j)-savr(j));
+            if isnan(wp(nrow+j)) then
+                ires = -1;
+                return;
+            end
+            iwp(nrow+j, 1) = i;
+            iwp(nrow+j, 2) = j;
+        end
+        nrow = nrow + neq;
+        y(i) = ysave;
+        ydot(i) = ypsave;
+    end
+endfunction
+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
+y0=1;t=2:6;t0=1;y0d=3;
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+atol=1.d-6;rtol=0;ng=2;
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
+assert_checkalmostequal(nn(1),2.47,0.001);
+y0=yy(1,2);y0d=yy(2,2);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
+assert_checkalmostequal(nn(1),2.5,0.001);
+y0=yy(1,1);y0d=yy(2,1);t0=nn(1);t=[3,4,5,6];
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res1,ng,"gr1",psol,pjac);
+assert_checkalmostequal(nn(1),2.53,0.001);
+//C
+//C-----------------------------------------------------------------------
+//C Second problem (Van Der Pol oscillator).
+//C The initial value problem is..
+//C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
+//C   Y1(0) = 2,  Y2(0) = 0,  0 .LE. T .LE. 200
+//C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
+//C The root function is  G = Y1.
+//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-----------------------------------------------------------------------
+%DAEOPTIONS=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]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
+assert_checkalmostequal(nn(1),81.163512,0.001);
+
+deff("[delta,ires]=res2(t,y,ydot)",...
+"ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,"jac2",ng,"gr2");
+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]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2");
+deff("s=gr2(t,y,yd)","s=y(1)")
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2);
+
+// Same problem, with psol and pjac example routines
+%DAEOPTIONS=list([],0,[],[],[],0,[],1,[],0,1,[],[],1);
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2","psol1","pjac1");
+assert_checkalmostequal(nn(1),81.163512,0.009);
+deff("s=gr2(t,y,yd)","s=y(1)")
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,"psol1","pjac1");
+assert_checkalmostequal(nn(1),81.163512,0.009);
+
+// 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]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,"gr2",psol,pjac);
+assert_checkalmostequal(nn(1),81.163512,0.003);
+deff("s=gr2(t,y,yd)","s=y(1)")
+[yy,nn]=dae("root2",[y0,y0d],t0,t,rtol,atol,res2,jac2,ng,gr2,psol,pjac);
+assert_checkalmostequal(nn(1),81.163512,0.003);
+%DAEOPTIONS=list([],0,[],[],[],0,[],0,[],0,0,[],[],1);
+
+//           Hot Restart
+[yy,nn,hotd]=dae("root2",[y0,y0d],t0,t,rtol,atol,"res2","jac2",ng,"gr2");
+t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
+[yy,nn,hotd]=dae("root2",[y01,y0d1],t01,t,rtol,atol,"res2","jac2",ng,"gr2",hotd);
+assert_checkalmostequal(nn(1),162.57763,0.004);
index 4e9e1b4..c68659c 100644 (file)
@@ -709,6 +709,7 @@ csim_3.png=c1b39c2350976c73e2b641855d55ddba
 csim_4.png=2795df06343733566b9081791d821f28
 csim_5.png=ac09bcd32af2fb43ae1028732ebbc1b7
 dae_1.png=f2902f78357a77a69afa81370a7fae7b
+dae_2.png=27ba222bd0de354067def52e330cfbad
 damp_en_US_1.png=e81873a345f84e42ce2c4f4779d3421d
 damp_fr_FR_1.png=e81873a345f84e42ce2c4f4779d3421d
 damp_ja_JP_1.png=e81873a345f84e42ce2c4f4779d3421d
diff --git a/scilab/modules/helptools/images/dae_2.png b/scilab/modules/helptools/images/dae_2.png
new file mode 100644 (file)
index 0000000..4ce2eea
Binary files /dev/null and b/scilab/modules/helptools/images/dae_2.png differ