Differential_equations: improvement for daskr 41/12041/3
Paul BIGNIER [Tue, 16 Jul 2013 12:58:18 +0000 (14:58 +0200)]
Similarly to scicos.c's 'jacpsol' and 'psol', doing preconditioner factorization in the evaluation function instead of the solving one.

Change-Id: I332594a2892348ad4f7f102852edb0697429cb9a

scilab/modules/differential_equations/help/en_US/daskr.xml
scilab/modules/differential_equations/help/fr_FR/daskr.xml [new file with mode: 0644]
scilab/modules/differential_equations/sci_gateway/c/Ex-daskr.c
scilab/modules/differential_equations/tests/unit_tests/daskr.dia.ref
scilab/modules/differential_equations/tests/unit_tests/daskr.tst

index 34c7f7c..e1f710a 100644 (file)
@@ -1,7 +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
+ * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
  *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
                 <listitem>
                     <para>
                         is either <literal>y0</literal> (<literal>ydot0</literal> is
-                        estimated by <literal>dassl</literal> with zero as first estimate)
+                        estimated by <literal>daskr</literal> with zero as first estimate)
                         or the matrix <literal>[y0 ydot0]</literal>.
-                        <literal>g(t,y0,ydot0)</literal> must be equal to zero. If you only
-                        know an estimate of <literal>ydot0</literal> set
+                        <literal>g(t, y0, ydot0)</literal> must be equal to zero.
+                        If you only know an estimate of <literal>ydot0</literal>, set
                         <literal>info(7)=1</literal>.
                     </para>
                     <variablelist>
@@ -44,8 +44,7 @@
                             <term>ydot0</term>
                             <listitem>
                                 <para>a real column vector of the time derivative of
-                                    <literal>y</literal> at <literal>t0</literal> (may be an
-                                    estimate).
+                                    <literal>y</literal> at <literal>t0</literal> (may be an estimate).
                                 </para>
                             </listitem>
                         </varlistentry>
                 <term>t</term>
                 <listitem>
                     <para>a real scalar or vector. Gives instants for which you want the
-                        solution. Note that you can get solution at each dassl's step point
+                        solution. Note that you can get solution at each daskr's step point
                         by setting <literal>info(2)=1</literal>.
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>nn</term>
-                <listitem>
-                    <para>
-                        a vector with two entries <literal>[times num]</literal>
-                        <literal>times</literal> is the value of the time at which the
-                        surface is crossed, <literal>num</literal> is the number of the
-                        crossed surface.
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
                 <term>atol, rtol</term>
                 <listitem>
                     <para>real scalars or column vectors of same size as
-                        <literal>y</literal>. <literal>atol, rtol</literal> give respectively
-                        absolute and relative error tolerances of solution. If vectors the
-                        tolerances are specified for each component of
+                        <literal>y</literal>. <literal>atol</literal> and <literal>rtol</literal> give respectively
+                        absolute and relative error tolerances of solution.
+                        If vectors, then the tolerances are specified for each component of
                         <literal>y</literal>.
                     </para>
                 </listitem>
                 <term>res</term>
                 <listitem>
                     <para>
-                        <link linkend="external">external</link> (function or list or string). Computes the value of
-                        <literal>g(t,y,ydot)</literal>. It may be :
+                        <link linkend="external">external</link> (function, list or string). Computes the value of
+                        <literal>g(t, y, ydot)</literal>. It may be :
                     </para>
                     <itemizedlist>
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>[r,ires]=res(t,y,ydot)</literal> and
-                                <literal>res</literal> must return the residual
-                                <literal>r=g(t,y,ydot)</literal> and error flag
+                                <literal>[r, ires] = res(t, y, ydot)</literal> and must return the residual
+                                <literal>r = g(t, y, ydot)</literal> and error flag
                                 <literal>ires</literal>. <literal>ires = 0</literal> if
                                 <literal>res</literal> succeeds to compute <literal>r</literal>,
-                                <literal>=-1</literal> if residual is locally not defined for
-                                <literal>(t,y,ydot)</literal>, <literal>=-2</literal> if
+                                <literal>ires = -1</literal> if residual is locally not defined for
+                                <literal>(t, y, ydot)</literal>, <literal>ires = -2</literal> if
                                 parameters are out of admissible range.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A list.</para>
-                            <para>This form allows to pass parameters other than t,y,ydot to
-                                the function. It must be as follows:
+                            <para>This form allows to pass parameters other than t, y, ydot to the function.
+                                It must be as follows:
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-list(res,x1,x2,...)
+list(res, x1, x2, ...)
  ]]></programlisting>
                             <para>where the calling sequence of the function
                                 <literal>res</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-r=res(t,y,ydot,x1,x2,...)
+r = res(t, y, ydot, x1, x2, ...)
  ]]></programlisting>
                             <para>
                                 <literal>res</literal> still returns
-                                <literal>r=g(t,y,ydot)</literal> as a function of
-                                <literal>(t,y,ydot,x1,x2,...)</literal>.
+                                <literal>r = g(t, y, ydot)</literal> as a function of
+                                <literal>(t, y, ydot, x1, x2, ...)</literal>.
                             </para>
                             <para>Warning: this form must not be used if there is no extra
-                                argument to pass to <literal>the function.</literal>
+                                argument to pass to the function.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A string.</para>
-                            <para>It must refer to the name of a C or Fortran subroutine
+                            <para>It must refer to the name of a C function or a Fortran subroutine
                                 linked with Scilab.
                             </para>
-                            <para>In C The calling sequence must be:</para>
+                            <para>In C, the calling sequence must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 void res(double *t, double y[], double yd[], double r[],
          int *ires, double rpar[], int ipar[])
  ]]></programlisting>
-                            <para>In Fortran it must be:</para>
+                            <para>In Fortran, it must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-subroutine res(t,y,yd,r,ires,rpar,ipar)
-double precision t, y(*),yd(*),r(*),rpar(*)
-integer ires,ipar(*)
+subroutine res(t, y, yd, r, ires, rpar, ipar)
+double precision t, y(*), yd(*),r(*),rpar(*)
+integer ires, ipar(*)
  ]]></programlisting>
                             <para>
-                                The <literal>rpar</literal> and <literal>ipar</literal> arrays must be present but cannot be
-                                used.
+                                The <literal>rpar</literal> and <literal>ipar</literal>
+                                arrays must be present but cannot be used.
                             </para>
                         </listitem>
                     </itemizedlist>
@@ -161,17 +148,16 @@ integer ires,ipar(*)
                 <term>jac</term>
                 <listitem>
                     <para>
-                        <link linkend="external">external</link> (function or list or string). Computes the value of
-                        <literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
+                        <link linkend="external">external</link> (function, list or string).
+                        Computes the value of <literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
                         <literal>cj</literal>.
                     </para>
                     <itemizedlist>
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>r=jac(t,y,ydot,cj)</literal> and the
-                                <literal>jac</literal> function must return
-                                <literal>r=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot</literal> where
+                                <literal>r = jac(t, y, ydot, cj)</literal> and must return
+                                <literal>r = dg(t, y, ydot)/dy + cj*dg(t, y, ydot)/dydot</literal> where
                                 <literal>cj</literal> is a real scalar.
                             </para>
                         </listitem>
@@ -179,34 +165,33 @@ integer ires,ipar(*)
                             <para>A list.</para>
                             <para>It must be as follows</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-list(jac,x1,x2,...)
+list(jac, x1, x2, ...)
  ]]></programlisting>
                             <para>where the calling sequence of the function
                                 <literal>jac</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-r=jac(t,y,ydot,cj,x1,x2,...)
+r=jac(t, y, ydot, cj, x1, x2,...)
  ]]></programlisting>
                             <para>
                                 <literal>jac</literal> still returns
                                 <literal>dg/dy + cj*dg/dydot</literal> as a function of
-                                <literal>(t,y,ydot,cj,x1,x2,...)</literal>.
+                                <literal>(t, y, ydot, cj, x1, x2, ...)</literal>.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A character string.</para>
-                            <para>It must refer to the name of a Fortran subroutine linked
-                                with scilab
+                            <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
                             </para>
-                            <para>In C The calling sequence must be:</para>
+                            <para>In C, the calling sequence must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 void jac(double *t, double y[], double yd[], double pd[],
          double *cj, double rpar[], int ipar[])
  ]]></programlisting>
-                            <para>In Fortran it must be:</para>
+                            <para>In Fortran, it must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-subroutine jac(t,y,yd,pd,cj,rpar,ipar)
-double precision t, y(*),yd(*),pd(*),cj,rpar(*)
+subroutine jac(t, y, yd, pd, cj, rpar, ipar)
+double precision t, y(*), yd(*), pd(*), cj, rpar(*)
 integer ipar(*)
  ]]></programlisting>
                         </listitem>
@@ -217,8 +202,8 @@ integer ipar(*)
                 <term>surf</term>
                 <listitem>
                     <para>
-                        <link linkend="external">external</link> (function or list or string). Computes the value of
-                        the column vector <literal>surf(t,y)</literal> with
+                        <link linkend="external">external</link> (function, list or string).
+                        Computes the value of the column vector <literal>surf(t, y)</literal> with
                         <literal>ng</literal> components. Each component defines a surface.
                         It may be defined by:
                     </para>
@@ -226,36 +211,35 @@ integer ipar(*)
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>surf(t,y)</literal>
+                                <literal>surf(t, y)</literal>
                             </para>
                         </listitem>
                         <listitem>
                             <para>A list.</para>
                             <para>It must be as follows</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-list(surf,x1,x2,...)
+list(surf, x1, x2, ...)
  ]]></programlisting>
                             <para>where the calling sequence of the function
                                 <literal>surf</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-r=surf(t,y,x1,x2,...)
+r = surf(t, y, x1, x2, ...)
  ]]></programlisting>
                         </listitem>
                         <listitem>
                             <para>A character string.</para>
-                            <para>It must refer to the name of a Fortran subroutine linked
-                                with scilab.
+                            <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab.
                             </para>
-                            <para>In C the calling sequence must be:</para>
+                            <para>In C, the calling sequence must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 void surf(int *ny, double *t, double y[], int *ng, double gout[])
  ]]></programlisting>
-                            <para>In Fortran it must be:</para>
+                            <para>In Fortran, it must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-subroutine surf(ny,t,y,ng,gout)
-double precision t, y(*),gout(*)
-integer ny,ng
+subroutine surf(ny, t, y, ng, gout)
+double precision t, y(*), gout(*)
+integer ny, ng
  ]]></programlisting>
                         </listitem>
                     </itemizedlist>
@@ -265,17 +249,16 @@ integer ny,ng
                 <term>info</term>
                 <listitem>
                     <para>
-                        list which contains <literal>14</literal> elements. Default
-                        value is <literal>list([],0,[],[],[],0,[],0,[],0,[],[],[],1)</literal>.
+                        list which contains <literal>14</literal> elements. Default value is
+                        <literal>list([], 0, [], [], [], 0, [], 0, [], 0, [], [], [], 1)</literal>.
                     </para>
                     <variablelist>
                         <varlistentry>
                             <term>info(1)</term>
                             <listitem>
                                 <para>real scalar which gives the maximum time for which
-                                    <literal>g</literal> is allowed to be evaluated or an empty
-                                    matrix <literal>[]</literal> if no limits imposed for
-                                    time.
+                                    <literal>g</literal> is allowed to be evaluated or an empty matrix
+                                    <literal>[]</literal> if no limits imposed for time.
                                 </para>
                             </listitem>
                         </varlistentry>
@@ -283,10 +266,10 @@ integer ny,ng
                             <term>info(2)</term>
                             <listitem>
                                 <para>
-                                    flag which indicates if <literal>dassl</literal> returns
-                                    its intermediate computed values (<literal>flag=1</literal>)
+                                    flag which indicates if <literal>daskr</literal> returns
+                                    its intermediate computed values (<literal>= 1</literal>)
                                     or only the user specified time point values
-                                    (<literal>flag=0</literal>).
+                                    (<literal>= 0</literal>).
                                 </para>
                             </listitem>
                         </varlistentry>
@@ -296,10 +279,8 @@ integer ny,ng
                                 <para>
                                     <literal>2</literal> components vector which give the
                                     definition <literal>[ml,mu]</literal> of band matrix computed
-                                    by <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) =
-                                        "dg(i)/dy(j)+cj*dg(i)/dydot(j)"
-                                    </literal>
-                                    . If <literal>jac</literal> returns a full matrix set
+                                    by <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"</literal>.
+                                    If <literal>jac</literal> returns a full matrix set
                                     <literal>info(3)=[]</literal>.
                                     Treat as dummmy if
                                     <literal>info(8)=1</literal>.
@@ -337,12 +318,12 @@ integer ny,ng
                             <listitem>
                                 <para>
                                     if ydot0 is set so that
-                                    <literal>g(t0,y0,ydot0)=0</literal> then set
+                                    <literal>g(t0, y0, ydot0) = 0</literal> then set
                                     <literal>info(7)=[]</literal>. Otherwise, set
                                     <literal>info(7)=[+-1, ..., +-1]</literal>, with
-                                    <literal>info(7)(i)=1</literal> if y(i) is a differential variable and
-                                    <literal>info(7)(i)=-1</literal> if y(i) is an algebraic variable
-                                    (if its derivatives do not appear explicitly in the function g(t,y,ydot)).
+                                    <literal>info(7)(i) = 1</literal> if y(i) is a differential variable and
+                                    <literal>info(7)(i) = -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>
@@ -350,8 +331,8 @@ integer ny,ng
                             <term>info(8)</term>
                             <listitem>
                                 <para>
-                                    direct / Krylov. Set <literal>info(8)=1</literal> and provide a routine in <literal>psol</literal>
-                                    if you want the solver to use Krylov iterations, else (dassl's direct method) set
+                                    direct / Krylov method. Set <literal>info(8)=1</literal> and provide a routine in <literal>psol</literal>
+                                    if you want the solver to use Krylov iterations, else (daskr's direct method) set
                                     <literal>info(8)=0</literal>.
                                 </para>
                             </listitem>
@@ -362,22 +343,23 @@ integer ny,ng
                                 <para>
                                     Krylov parameters. Treat as dummy argument if you have set
                                     <literal>info(8)=0</literal>. Otherwise, set
+                                    <literal>info(9)=[]</literal> or
                                     <literal>info(9)=[maxl kmp nrmax epli]</literal>, where:
                                 </para>
                                 <para>
-                                    - maxl = maximum number of iterations in the SPIGMR algorithm (default
-                                    <literal>min(5,neq)</literal>),
+                                    - maxl = maximum number of iterations in the GMRes algorithm (default
+                                    <literal>min(5, neq)</literal>),
                                 </para>
                                 <para>
-                                    - kmp = number of vectors on which orthogonalization is done in the SPIGMR algorithm
+                                    - kmp = number of vectors on which orthogonalization is done in the GMRes algorithm
                                     (default maxl),
                                 </para>
                                 <para>
-                                    - nrmax = maximum number of restarts of the SPIGMR algorithm per nonlinear iteration
+                                    - nrmax = maximum number of restarts of the GMRes algorithm per nonlinear iteration
                                     (default <literal>5</literal>),
                                 </para>
                                 <para>
-                                    - epli = convergence test constant in SPIGMR algorithm (default <literal>0.05</literal>).
+                                    - epli = convergence test constant in GMRes algorithm (default <literal>0.05</literal>).
                                 </para>
                             </listitem>
                         </varlistentry>
@@ -412,9 +394,9 @@ integer ny,ng
                                     if you wish to control errors locally on all the variables then set
                                     <literal>info(12)=[]</literal>. Otherwise, set
                                     <literal>info(12)=[+-1, ..., +-1]</literal>, with
-                                    <literal>info(12)(i)=1</literal> if y(i) is a differential variable and
-                                    <literal>info(12)(i)=-1</literal> if y(i) is an algebraic variable
-                                    (if its derivatives do not appear explicitly in the function g(t,y,ydot)).
+                                    <literal>info(12)(i) = 1</literal> if y(i) is a differential variable and
+                                    <literal>info(12)(i) = -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>
@@ -424,6 +406,7 @@ integer ny,ng
                                 <para>
                                     heuristic parameters. Treat as dummy argument if
                                     <literal>info(7)=[]</literal>. Otherwise, set
+                                    <literal>info(13)=[]</literal> or
                                     <literal>info(13)=[mxnit mxnj mxnh lsoff stptol epinit]</literal>, where:
                                 </para>
                                 <para>
@@ -435,7 +418,7 @@ integer ny,ng
                                     <literal>6</literal> if <literal>info(8)=0</literal>, <literal>2</literal> otherwise),
                                 </para>
                                 <para>
-                                    - mxnh = maximum number of values of the artificial stepsize parameter H to be tried if info(7) &#8800; [] (default
+                                    - mxnh = maximum number of values of the artificial stepsize parameter <literal>h</literal> to be tried if info(7) &#8800; [] (default
                                     <literal>5</literal>),
                                 </para>
                                 <para>
@@ -467,30 +450,29 @@ integer ny,ng
                 <term>psol</term>
                 <listitem>
                     <para>
-                        <link linkend="external">external</link> (function or list or string). Solves a linear system
-                        <literal>P*x = b</literal>, with P being the preconditioner that routine <literal>pjac</literal>
-                        computed beforehand and stocked in wp and iwp.
+                        <link linkend="external">external</link> (function, list or string). Solves a linear system
+                        <literal>P*x = b</literal>, with P being the factored preconditioner that routine <literal>pjac</literal>
+                        computed beforehand and stored in <literal>wp</literal> and <literal>iwp</literal>.
                     </para>
                     <itemizedlist>
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
-                                <literal>[r, ier] = psol(wp, iwp, b)</literal> and the
-                                <literal>psol</literal> function must return the solution of that system in
+                                <literal>[r, ier] = psol(wp, iwp, b)</literal> and must return the solution of the system in
                                 <literal>r</literal> and an error flag <literal>ier</literal>.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A list.</para>
-                            <para>It must be as follows</para>
+                            <para>It must be as follows:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-list(psol,x1,x2,...)
+list(psol, x1, x2, ...)
  ]]></programlisting>
-                            <para>where the calling sequence of the function
-                                <literal>psol</literal> is now
+                            <para>
+                                where the calling sequence of <literal>psol</literal> is now
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-psol(wp,iwp,b,x1,x2,...)
+psol(wp, iwp, b, x1, x2, ...)
  ]]></programlisting>
                             <para>
                                 <literal>psol</literal> still returns the solution in <literal>r</literal>.
@@ -498,23 +480,23 @@ psol(wp,iwp,b,x1,x2,...)
                         </listitem>
                         <listitem>
                             <para>A character string.</para>
-                            <para>It must refer to the name of a Fortran subroutine linked
-                                with scilab
+                            <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
                             </para>
-                            <para>In C The calling sequence must be:</para>
+                            <para>In C, the calling sequence must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 void psol (int*neq, double*t, double*y, double*ydot, double*savr,
   double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
  ]]></programlisting>
-                            where the arrays <literal>wp</literal> and <literal>iwp</literal> contain matrix elements of preconditioner
-                            <literal>P</literal> in compressed row format.
-                            <para>In Fortran it must be:</para>
+                            where the arrays <literal>wp</literal> and <literal>iwp</literal> contain matrix elements of LU-factored preconditioner
+                            <literal>P</literal>, <literal>wp</literal> being the values and
+                            <literal>iwp</literal> the pivots used in the factorization.
+                            <para>In Fortran, it must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
   wp, iwp, b, eplin, ier, rpar, ipar)
-double precision t,y(*),ydot(*),savr(*),wk(*),cj,wght(*),wp(*),
-  b(*),eplin,rpar(*)
-integer neq,iwp(*),ier,ipar(*)
+double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
+  b(*), eplin, rpar(*)
+integer neq, iwp(*), ier, ipar(*)
  ]]></programlisting>
                         </listitem>
                     </itemizedlist>
@@ -524,47 +506,46 @@ integer neq,iwp(*),ier,ipar(*)
                 <term>pjac</term>
                 <listitem>
                     <para>
-                        <link linkend="external">external</link> (function or list or string). Computes the value of
+                        <link linkend="external">external</link> (function, list or string). Computes the value of
                         <literal>dg/dy + cj*dg/dydot</literal> for a given value of parameter
-                        <literal>cj</literal> and LU-factorizes it in double and int arrays.
+                        <literal>cj</literal> and LU-factorizes it in two arrays, real and integer.
                     </para>
                     <itemizedlist>
                         <listitem>
                             <para>A Scilab function.</para>
                             <para>Its calling sequence must be
                                 <literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> and in return,
-                                the arrays wp and iwp must contain all preconditioner information in compressed sparse row format.
+                                the arrays <literal>wp</literal> and <literal>iwp</literal> must contain all factored preconditioner information.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A list.</para>
                             <para>It must be as follows</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-list(pjac,x1,x2,...)
+list(pjac, x1, x2, ...)
  ]]></programlisting>
-                            <para>where the calling sequence of the function
-                                <literal>pjac</literal> is now
+                            <para>
+                                where the calling sequence of <literal>pjac</literal> is
                             </para>
                             <programlisting role="no-scilab-exec"><![CDATA[
-pjac(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)
+pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
  ]]></programlisting>
                             <para>
                                 <literal>pjac</literal> still returns factorized
                                 <literal>dg/dy + cj*dg/dydot</literal> as a function of
-                                <literal>(neq,t,y,ydot,h,cj,rewt,savr,x1,x2,...)</literal>.
+                                <literal>(neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...)</literal>.
                             </para>
                         </listitem>
                         <listitem>
                             <para>A character string.</para>
-                            <para>It must refer to the name of a Fortran subroutine linked
-                                with scilab
+                            <para>It must refer to the name of a C function or a Fortran subroutine linked with Scilab
                             </para>
-                            <para>In C The calling sequence must be:</para>
+                            <para>In C, the calling sequence must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
 double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
  ]]></programlisting>
-                            <para>In Fortran it must be:</para>
+                            <para>In Fortran, it must be:</para>
                             <programlisting role="no-scilab-exec"><![CDATA[
 subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
   wk, h, cj, wp, iwp, ier, rpar, ipar)
@@ -580,7 +561,7 @@ integer ires, neq, iwp(*), ier, ipar(*)
                 <term>hd</term>
                 <listitem>
                     <para>
-                        real vector which allows to store the <literal>dassl</literal>
+                        real vector which allows to store the <literal>daskr</literal>
                         context and to resume integration.
                     </para>
                 </listitem>
@@ -589,8 +570,20 @@ integer ires, neq, iwp(*), ier, ipar(*)
                 <term>r</term>
                 <listitem>
                     <para>
-                        real matrix . Each column is the vector <literal>[t;x(t);xdot(t)]</literal> where
-                        <literal>t</literal> is time index for which the solution had been computed.
+                        real matrix. Each column is the vector <literal>[t; x(t); xdot(t)]</literal> where
+                        <literal>t</literal> is the time index for which the solution has been computed,
+                        <literal>x(t)</literal> is the value of the computed solution,
+                        <literal>xdot(t)</literal> is the derivative of the computed solution.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>nn</term>
+                <listitem>
+                    <para>
+                        a vector with two entries <literal>[times num]</literal>,
+                        <literal>times</literal> is the value of the time at which the surface is crossed,
+                        <literal>num</literal> is the number of the crossed surface.
                     </para>
                 </listitem>
             </varlistentry>
@@ -598,10 +591,10 @@ integer ires, neq, iwp(*), ier, ipar(*)
     </refsection>
     <refsection>
         <title>Description</title>
-        <para>Solution of the implicit differential equation.</para>
+        <para>Solution of the implicit differential equation:</para>
         <programlisting role="no-scilab-exec"><![CDATA[
-g(t,y,ydot) = 0
-y(t0) = y0  and   ydot(t0) = ydot0
+g(t, y, ydot) = 0
+y(t0) = y0  and  ydot(t0) = ydot0
  ]]></programlisting>
         <para>Returns the surface crossing instants and the number of the surface
             reached in <literal>nn</literal>.
@@ -611,16 +604,16 @@ y(t0) = y0  and   ydot(t0) = ydot0
     <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[
-// dy/dt = ((2*log(y)+8)/t -5)*y,  y(1) = 1,  1 <= t <=6
-// g1 = ((2*log(y)+8)/t - 5)*y
+// dy/dt = ((2*log(y)+8)/t-5)*y,  y(1) = 1,  1 <= t <= 6
+// g1 = ((2*log(y)+8)/t-5)*y
 // g2 = log(y) - 2.2491
 y0 = 1; t = 2:6; t0 = 1; y0d = 3;
 atol = 1.d-6; rtol = 0; ng = 2;
 
-deff('[delta,ires] = res1(t,y,ydot)','ires = 0; delta = ydot-((2*log(y)+8)/t-5)*y')
-deff('[rts] = gr1(t,y)','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)', 'rts = [((2*log(y)+8)/t-5)*y; log(y)-2.2491]')
 
-[yy,nn] = daskr([y0,y0d],t0,t,atol,rtol,res1,ng,gr1);
+[yy, nn] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
 nn
 
 // Should return nn = [2.4698972 2]
diff --git a/scilab/modules/differential_equations/help/fr_FR/daskr.xml b/scilab/modules/differential_equations/help/fr_FR/daskr.xml
new file mode 100644 (file)
index 0000000..f324865
--- /dev/null
@@ -0,0 +1,660 @@
+<?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
+ *
+ * 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
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="daskr" xml:lang="fr">
+    <refnamediv>
+        <refname>daskr</refname>
+        <refpurpose>solveur de DAE avec traversées de zéros</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Séquence d'appel</title>
+        <synopsis>[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])</synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Paramètres</title>
+        <variablelist>
+            <varlistentry>
+                <term>x0</term>
+                <listitem>
+                    <para>
+                        représente soit <literal>y0</literal> (<literal>ydot0</literal> sera
+                        estimé par <literal>daskr</literal> avec zéro comme première estimation),
+                        soir la matrice <literal>[y0 ydot0]</literal>.
+                        <literal>g(t, y0, ydot0)</literal> doit être égal à zéro. Si vous ne connaissez
+                        qu'une estimation de <literal>ydot0</literal>, assignez
+                        <literal>info(7)=1</literal>.
+                    </para>
+                    <variablelist>
+                        <varlistentry>
+                            <term>y0</term>
+                            <listitem>
+                                <para>vecteur colonne réel des conditions initiales.</para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>ydot0</term>
+                            <listitem>
+                                <para>vecteur colonne réel de la dérivée en temps de
+                                    <literal>y</literal> à <literal>t0</literal> (peut être une estimation.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                    </variablelist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>t0</term>
+                <listitem>
+                    <para>réel, temps initial.</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>t</term>
+                <listitem>
+                    <para>réel, scalaire ou vecteur. Temps auxquels la solution est désirée.
+                        Notez que vous pouvez obtenir la solution à chaque étape de daskr en fixant
+                        <literal>info(2)=1</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>atol, rtol</term>
+                <listitem>
+                    <para>réels scalaires ou vecteurs colonnes de même taille que
+                        <literal>y</literal>. <literal>atol</literal> et <literal>rtol</literal> représentent
+                        les tolérances d'erreur absolue et relative de la solution.
+                        Si ce sont des vecteurs, alors les tolérances sont spécifiées pour chaque composante de
+                        <literal>y</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>res</term>
+                <listitem>
+                    <para>
+                        <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+                        Calcule la valeur de
+                        <literal>g(t, y, ydot)</literal>. Elle peut être :
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>Une fonction Scilab.</para>
+                            <para>Sa séquence d'appel doit être
+                                <literal>[r, ires] = res(t, y, ydot)</literal> et doit retourner le résidu
+                                <literal>r = g(t, y, ydot)</literal> et un drapeau d'erreur
+                                <literal>ires</literal>. <literal>ires = 0</literal> si
+                                <literal>res</literal> a correctement calculé <literal>r</literal>,
+                                <literal>ires = -1</literal> si le résidu est localement non défini pour
+                                <literal>(t, y, ydot)</literal>, <literal>ires = -2</literal> si
+                                des paramètres sont hors du champ admissible.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une liste.</para>
+                            <para>Cette forme permet de passer des paramètres autres que t, y, ydot à la fonction.
+                                Elle doit se présenter comme suit :
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+list(res, x1, x2, ...)
+ ]]></programlisting>
+                            <para>où la séquence d'appel de la fonction
+                                <literal>res</literal> est maintenant
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+r = res(t, y, ydot, x1, x2, ...)
+ ]]></programlisting>
+                            <para>
+                                <literal>res</literal> retourne toujours
+                                <literal>r = g(t, y, ydot)</literal> comme fonction de
+                                <literal>(t, y, ydot, x1, x2, ...)</literal>.
+                            </para>
+                            <para>Attention : cette forme ne doit pas être utilisée s'il n'y pas
+                                d'argument additionnel à passer à la fonction.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une chaîne de caractères.</para>
+                            <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran
+                                reliée à Scilab.
+                            </para>
+                            <para>En C, la séquence d'appel doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+void res(double *t, double y[], double yd[], double r[],
+         int *ires, double rpar[], int ipar[])
+ ]]></programlisting>
+                            <para>En Fortran, elle doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+subroutine res(t, y, yd, r, ires, rpar, ipar)
+double precision t, y(*), yd(*),r(*),rpar(*)
+integer ires, ipar(*)
+ ]]></programlisting>
+                            <para>
+                                Les tableaux <literal>rpar</literal> et <literal>ipar</literal>
+                                doivent être présents mais ne peuvent pas être utilisés.
+                            </para>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>jac</term>
+                <listitem>
+                    <para>
+                        <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+                        Calcule la valeur de <literal>dg/dy + cj*dg/dydot</literal> pour une valeur donnée du paramètre
+                        <literal>cj</literal>.
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>Une fonction Scilab.</para>
+                            <para>Sa séquence d'appel doit être
+                                <literal>r = jac(t, y, ydot, cj)</literal> et doit retourner
+                                <literal>r = dg(t, y, ydot)/dy + cj*dg(t, y, ydot)/dydot</literal> où
+                                <literal>cj</literal> est un scalaire réel.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une liste.</para>
+                            <para>Elle doit se présenter comme suit :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+list(jac, x1, x2, ...)
+ ]]></programlisting>
+                            <para>où la séquence d'appel de la fonction
+                                <literal>jac</literal> est désormais
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+r=jac(t, y, ydot, cj, x1, x2,...)
+ ]]></programlisting>
+                            <para>
+                                <literal>jac</literal> retourne toujours
+                                <literal>dg/dy + cj*dg/dydot</literal> comme fonction de
+                                <literal>(t, y, ydot, cj, x1, x2,...)</literal>.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une chaîne de caractères.</para>
+                            <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.
+                            </para>
+                            <para>En C, la séquence d'appel doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+void jac(double *t, double y[], double yd[], double pd[],
+         double *cj, double rpar[], int ipar[])
+ ]]></programlisting>
+                            <para>En Fortran, elle doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+subroutine jac(t, y, yd, pd, cj, rpar, ipar)
+double precision t, y(*), yd(*), pd(*), cj, rpar(*)
+integer ipar(*)
+ ]]></programlisting>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>surf</term>
+                <listitem>
+                    <para>
+                        <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+                        Calcule la valeur du vecteur colonne <literal>surf(t, y)</literal> à
+                        <literal>ng</literal> composantes. Chaque composante représente une surface.
+                        Elle doit être définie comme suit :
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>Une fonction Scilab.</para>
+                            <para>Sa séquence d'appel doit être
+                                <literal>surf(t, y)</literal>
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une liste.</para>
+                            <para>Elle doit se présenter comme suit :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+list(surf, x1, x2, ...)
+ ]]></programlisting>
+                            <para>où la séquence d'appel de la fonction
+                                <literal>surf</literal> est maintenant
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+r = surf(t, y, x1, x2, ...)
+ ]]></programlisting>
+                        </listitem>
+                        <listitem>
+                            <para>Une chaîne de caractères.</para>
+                            <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.
+                            </para>
+                            <para>En C, la séquence d'appel doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+void surf(int *ny, double *t, double y[], int *ng, double gout[])
+ ]]></programlisting>
+                            <para>En Fortran, elle doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+subroutine surf(ny, t, y, ng, gout)
+double precision t, y(*), gout(*)
+integer ny, ng
+ ]]></programlisting>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>info</term>
+                <listitem>
+                    <para>
+                        liste contenant <literal>14</literal> éléments. La valeur par défaut est
+                        <literal>list([], 0, [], [], [], 0, [], 0, [], 0, [], [], [], 1)</literal>.
+                    </para>
+                    <variablelist>
+                        <varlistentry>
+                            <term>info(1)</term>
+                            <listitem>
+                                <para>réel scalaire donnant le temps maximal pour lequel
+                                    <literal>g</literal> peut être évalué ou une matrice vide
+                                    <literal>[]</literal> si aucune limite de temps n'est imposée.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(2)</term>
+                            <listitem>
+                                <para>
+                                    drapeau indiquant si <literal>daskr</literal> retourne
+                                    ses valeurs intermédiaires calculées (<literal>= 1</literal>)
+                                    ou seulement les temps indiqués par l'utilisateur
+                                    (<literal>= 0</literal>).
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(3)</term>
+                            <listitem>
+                                <para>
+                                    vecteur de deux éléments donnant la définition
+                                    <literal>[ml,mu]</literal> de la matrice bande calculeé par
+                                    <literal>jac</literal>; <literal>r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"</literal>.
+                                    Si <literal>jac</literal> retourne une matrice pleine, fixer
+                                    <literal>info(3)=[]</literal>.
+                                    Inutile si
+                                    <literal>info(8)=1</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(4)</term>
+                            <listitem>
+                                <para>réel scalaire donnant la taille maximale du pas. Fixer
+                                    <literal>info(4)=[]</literal> si illimité.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(5)</term>
+                            <listitem>
+                                <para>réel scalaire donnant le pas initial. Fixer
+                                    <literal>info(5)=[]</literal> si non spécifié.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(6)</term>
+                            <listitem>
+                                <para>
+                                    fixer <literal>info(6)=1</literal> si la solution est
+                                    non-négative, sinon fixer
+                                    <literal>info(6)=0</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(7)</term>
+                            <listitem>
+                                <para>
+                                    si ydot0 est fixé tel que
+                                    <literal>g(t0, y0, ydot0) = 0</literal>, alors fixer
+                                    <literal>info(7)=[]</literal>. Sinon, fixer
+                                    <literal>info(7)=[+-1, ..., +-1]</literal>, avec
+                                    <literal>info(7)(i) = 1</literal> si y(i) est une variable différentielle et
+                                    <literal>info(7)(i) = -1</literal> si y(i) est une variable algébrique
+                                    (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(8)</term>
+                            <listitem>
+                                <para>
+                                    méthode directe / Krylov. Fixer <literal>info(8)=1</literal> et founrnir une routine <literal>psol</literal>
+                                    si vous souhaitez que le solveur utilise des itérations de Krylov, sinon (méthode directe) fixer
+                                    <literal>info(8)=0</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(9)</term>
+                            <listitem>
+                                <para>
+                                    paramètres de Krylov. Inutile si vous avez fixé
+                                    <literal>info(8)=0</literal>. Sinon, fixer
+                                    <literal>info(9)=[]</literal> ou
+                                    <literal>info(9)=[maxl kmp nrmax epli]</literal>, où :
+                                </para>
+                                <para>
+                                    - maxl = nombre maximal d'itérations de l'algorithme GMRes (par défaut
+                                    <literal>min(5, neq)</literal>),
+                                </para>
+                                <para>
+                                    - kmp = nombre de vecteurs sur lesquels l'orthogonalisation est faite dans GMRes
+                                    (par défaut maxl),
+                                </para>
+                                <para>
+                                    - nrmax = nombre maximal de redémarrages de GMRes par intération non-linéaire
+                                    (par défaut <literal>5</literal>),
+                                </para>
+                                <para>
+                                    - epli = constante du test de convergence de GMRes (par défaut <literal>0.05</literal>).
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(10)</term>
+                            <listitem>
+                                <para>
+                                    conditions initiales. A ignorer si
+                                    <literal>info(7)=[]</literal>. Fixer
+                                    <literal>info(10)=1</literal> si le solveur doit s'arrêter après
+                                    le calcul des valeurs initiales, sinon fixer
+                                    <literal>info(10)=0</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(11)</term>
+                            <listitem>
+                                <para>
+                                    routine pour le calcul et la factorisation LU du préconditionneur pour <literal>psol</literal>.
+                                    Inutile si <literal>info(8)=0</literal>. Fixer
+                                    <literal>info(11)=1</literal> et fournir une routine <literal>pjac</literal> si l'<link linkend="external">external</link>
+                                    <literal>psol</literal> doit utiliser une routine spécifique, sinon fixer
+                                    <literal>info(11)=0</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(12)</term>
+                            <listitem>
+                                <para>
+                                    si vous souhaitez contrôler l'erreur localement sur toutes les variables, fixez
+                                    <literal>info(12)=[]</literal>. Sinon, fixez
+                                    <literal>info(12)=[+-1, ..., +-1]</literal>, avec
+                                    <literal>info(12)(i) = 1</literal> si y(i) est une variable différentielle et
+                                    <literal>info(12)(i) = -1</literal> si y(i) est une variable algébrique
+                                    (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(13)</term>
+                            <listitem>
+                                <para>
+                                    paramètres heuristiques. Ignorer si
+                                    <literal>info(7)=[]</literal>. Sinon, fixer
+                                    <literal>info(13)=[]</literal> ou
+                                    <literal>info(13)=[mxnit mxnj mxnh lsoff stptol epinit]</literal>, où :
+                                </para>
+                                <para>
+                                    - mxnit = nombre maximal d'itérations de Newton par évaluation du préconditionneur (par défaut
+                                    <literal>5</literal> si <literal>info(8)=0</literal>, <literal>15</literal> sinon),
+                                </para>
+                                <para>
+                                    - mxnj = nombre maximal d'évaluations du préconditionneur (par défaut
+                                    <literal>6</literal> si <literal>info(8)=0</literal>, <literal>2</literal> sinon),
+                                </para>
+                                <para>
+                                    - mxnh = nombre maximal de valeurs artificielles du pas <literal>h</literal> à tenter si info(7) &#8800; [] (par défaut
+                                    <literal>5</literal>),
+                                </para>
+                                <para>
+                                    - lsoff = drapeau pour désactiver l'algorithme de recherche linéaire (lsoff = 0 pour activer, lsoff = 1 pour désactiver)
+                                    (par défaut <literal>0</literal>),
+                                </para>
+                                <para>
+                                    - stptol = pas minimal (dimmensionné) dans l'algorithme de recherche linéaire (par défaut <literal>(unit roundoff)^(2/3)</literal>),
+                                </para>
+                                <para>
+                                    - epinit = facteur déterminant dans le test de convergence de l'itération Newton (par défaut <literal>0.01</literal>).
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                        <varlistentry>
+                            <term>info(14)</term>
+                            <listitem>
+                                <para>
+                                    verbosité. Fixer <literal>info(14)=1</literal> pour une information minimale,
+                                    <literal>info(14)=2</literal> pour une information maximale, sinon fixer
+                                    <literal>info(14)=0</literal>.
+                                </para>
+                            </listitem>
+                        </varlistentry>
+                    </variablelist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>psol</term>
+                <listitem>
+                    <para>
+                        <link linkend="external">external</link> (fonction, liste ou chaîne de caractères).
+                        Résout un système linéraire <literal>P*x = b</literal>,
+                        où P est le préconditionneur LU-factorisé que <literal>pjac</literal>
+                        a calculé auparavant et stocké dans <literal>wp</literal> et <literal>iwp</literal>.
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>Une fonction Scilab.</para>
+                            <para>Sa séquence d'appel doit être
+                                <literal>[r, ier] = psol(wp, iwp, b)</literal> et doit retourner la solution du système dans
+                                <literal>r</literal> et un drapeau d'erreur <literal>ier</literal>.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une liste.</para>
+                            <para>Elle doit se présenter comme suit :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+list(psol, x1, x2, ...)
+ ]]></programlisting>
+                            <para>
+                                où la séquence d'appel de <literal>psol</literal> est
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+psol(wp, iwp, b, x1, x2, ...)
+ ]]></programlisting>
+                            <para>
+                                <literal>psol</literal> retourne toujours la solution dans <literal>r</literal>.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une chaîne de caractères.</para>
+                            <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab
+                            </para>
+                            <para>En C, la séquence d'appel doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+void psol (int*neq, double*t, double*y, double*ydot, double*savr,
+  double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
+ ]]></programlisting>
+                            où les vecteurs <literal>wp</literal> et <literal>iwp</literal> contiennent le préconditionneur LU-factorisé
+                            <literal>P</literal>, <literal>wp</literal> représentant les valeurs et
+                            <literal>iwp</literal> les pivots utilisés dans la factorisation.
+                            <para>En Fortran, elle doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
+  wp, iwp, b, eplin, ier, rpar, ipar)
+double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
+  b(*), eplin, rpar(*)
+integer neq, iwp(*), ier, ipar(*)
+ ]]></programlisting>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>pjac</term>
+                <listitem>
+                    <para>
+                        <link linkend="external">external</link> (fonction, liste ou chaîne de caractères). Calcule la valeur de
+                        <literal>dg/dy + cj*dg/dydot</literal> pour une valeur donnée du paramètre
+                        <literal>cj</literal> et la LU-factorise en deux vecteurs, réel et entier.
+                    </para>
+                    <itemizedlist>
+                        <listitem>
+                            <para>Une fonction Scilab.</para>
+                            <para>Sa séquence d'appel doit être
+                                <literal>[wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)</literal> et en retour,
+                                les vecteurs <literal>wp</literal> et <literal>iwp</literal>
+                                doivent contenir toutes les informations du préconditionneur factorisé.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une liste.</para>
+                            <para>Elle doit se présenter comme suit :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+list(pjac, x1, x2, ...)
+ ]]></programlisting>
+                            <para>
+                                où la séquence d'appel de <literal>pjac</literal> est
+                            </para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
+ ]]></programlisting>
+                            <para>
+                                <literal>pjac</literal> retourne toujours
+                                <literal>dg/dy + cj*dg/dydot</literal> comme fonction de
+                                <literal>(neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...)</literal>.
+                            </para>
+                        </listitem>
+                        <listitem>
+                            <para>Une chaîne de caractères.</para>
+                            <para>Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab
+                            </para>
+                            <para>En C, la séquence d'appel doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
+double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
+ ]]></programlisting>
+                            <para>En Fortran, elle doit être :</para>
+                            <programlisting role="no-scilab-exec"><![CDATA[
+subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
+  wk, h, cj, wp, iwp, ier, rpar, ipar)
+double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
+  wk(*), h, cj, wp(*), rpar(*)
+integer ires, neq, iwp(*), ier, ipar(*)
+ ]]></programlisting>
+                        </listitem>
+                    </itemizedlist>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>hd</term>
+                <listitem>
+                    <para>
+                        vecteur réel servant à stocker le contexte de
+                        <literal>daskr</literal> et reprendre l'intégration.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>r</term>
+                <listitem>
+                    <para>
+                        matrice réelle. Chaque colonne est le vecteur <literal>[t; x(t); xdot(t)]</literal> où
+                        <literal>t</literal> est l'indice en temps aulequel la solution a été calculée,
+                        <literal>x(t)</literal> est la valeur de la solution calculée,
+                        <literal>xdot(t)</literal> est la dérivée de la solution calculée.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>nn</term>
+                <listitem>
+                    <para>
+                        vecteur à deux entrées <literal>[times num]</literal>,
+                        <literal>times</literal> est la valeur du temps auquel la surface est traversée,
+                        <literal>num</literal> est le nombre de surfaces traversées.
+                    </para>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Description</title>
+        <para>Solution de l'équation différentielle implicite :</para>
+        <programlisting role="no-scilab-exec"><![CDATA[
+g(t, y, ydot) = 0
+y(t0) = y0  et  ydot(t0) = ydot0
+ ]]></programlisting>
+        <para>Retourne les temps de traversée de surface et le nombre de surfaces traversées
+            dans <literal>nn</literal>.
+        </para>
+        <para>Des exemples détaillés se trouvent dans SCI/modules/differential_equations/tests/unit_tests/daskr.tst</para>
+    </refsection>
+    <refsection>
+        <title>Exemples</title>
+        <programlisting role="example"><![CDATA[
+// dy/dt = ((2*log(y)+8)/t-5)*y,  y(1) = 1,  1 <= t <= 6
+// g1 = ((2*log(y)+8)/t-5)*y
+// g2 = log(y) - 2.2491
+y0 = 1; t = 2:6; t0 = 1; y0d = 3;
+atol = 1.d-6; rtol = 0; ng = 2;
+
+deff('[delta, ires] = res1(t, y, ydot)', 'ires = 0; delta = ydot-((2*log(y)+8)/t-5)*y')
+deff('[rts] = gr1(t, y)', 'rts = [((2*log(y)+8)/t-5)*y; log(y)-2.2491]')
+
+[yy, nn] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
+nn
+
+// Retourne nn = [2.4698972 2]
+ ]]></programlisting>
+    </refsection>
+    <refsection role="see also">
+        <title>Voir aussi</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="ode">ode</link>
+            </member>
+            <member>
+                <link linkend="dasrt">dasrt</link>
+            </member>
+            <member>
+                <link linkend="dassl">dassl</link>
+            </member>
+            <member>
+                <link linkend="impl">impl</link>
+            </member>
+            <member>
+                <link linkend="fort">fort</link>
+            </member>
+            <member>
+                <link linkend="link">link</link>
+            </member>
+            <member>
+                <link linkend="external">external</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>Historique</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>daskr ajouté</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 11f9f7c..2a5fb9a 100644 (file)
@@ -3,37 +3,46 @@
 #include "machine.h"
 #include "core_math.h"
 
-extern void C2F(dgesv)(int*, int*, double*, int*, int*, double*, int*, int*);
-extern double C2F(dlamch)(const char*);
-typedef void (*resfunc)(double*, double*, double*, double*, int*, double*, int*);
+extern void   C2F(dgefa) (double *A, int *lead_dim_A, int *n, int *ipivots, int *info);
+extern void   C2F(dgesl) (double *A, int *lead_dim_A, int *n, int *ipivots, double *B, int *job);
+extern double C2F(dlamch) (const char*);
+typedef void (*resfunc) (double*, double*, double*, double*, int*, double*, int*);
 
-void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *actual, double *actualP,
-                double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
-                int *ier, double *rpar, int *ipar)
+void C2F(pjac1) (resfunc res, int *ires, int *neq, double *tOld, double *actual, double *actualP,
+                 double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
+                 int *ier, double *rpar, int *ipar);
+void C2F(psol1) (int *neq, double *tOld, double *actual, double *actualP,
+                 double *savr, double *wk, double *cj, double *wght, double *wp,
+                 int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2);
+
+void C2F(pjac1) (resfunc res, int *ires, int *neq, double *tOld, double *actual, double *actualP,
+                 double *rewt, double *savr, double *wk, double *h, double *cj, double *wp, int *iwp,
+                 int *ier, double *rpar, int *ipar)
 {
     int i = 0;
     int j = 0;
+    int info = 0;
     int nrow = 0;
     double tx = 0;
     double del = 0;
     double delinv = 0;
     double ysave = 0;
     double ypsave = 0;
-    double* e = NULL;
+    double * e = NULL;
 
-    int neq = *nequations;
     double SQuround = sqrt(C2F(dlamch)("P"));
 
     tx = *tOld;
 
-    e = (double*) calloc(neq, sizeof(double));
+    e = (double *) calloc(*neq, sizeof(double));
 
-    for (i = 0 ; i < neq ; ++i)
+    for (i = 0 ; i < *neq ; ++i)
     {
-        del = Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
+        del =  Max(SQuround * Max(fabs(actual[i]), fabs(*h * actualP[i])), 1. / rewt[i]);
         del *= (*h * actualP[i] >= 0) ? 1 : -1;
-        ysave   =  actual[i];
-        ypsave  =  actualP[i];
+        del =  (actual[i] + del) - actual[i];
+        ysave  = actual[i];
+        ypsave = actualP[i];
         actual[i]  += del;
         actualP[i] += *cj * del;
         res(&tx, actual, actualP, e, ires, rpar, ipar);
@@ -46,57 +55,50 @@ void C2F(pjac1)(resfunc res, int *ires, int *nequations, double *tOld, double *a
         }
 
         delinv = 1. / del;
-        for (j = 0 ; j < neq ; j++)
+        for (j = 0 ; j < *neq ; j++)
         {
             wp[nrow + j] = (e[j] - savr[j]) * delinv;
             /* NaN test */
             if (ISNAN(wp[nrow + j]))
             {
                 *ier = -1;
-                free (e);
+                free(e);
                 return;
             }
-            iwp[nrow + j] = i + 1;
-            iwp[nrow + j + neq * neq] = j + 1;
         }
-        nrow       += neq;
+        nrow       += *neq;
         actual[i]  =  ysave;
         actualP[i] =  ypsave;
     }
 
+    /* Proceed to LU factorization of P. */
+    C2F(dgefa) (wp, neq, neq, iwp, &info);
+    if (info != 0)
+    {
+        *ier = -1;
+    }
+
     free(e);
 }
 
-void C2F(psol1)(int *nequations, double *tOld, double *actual, double *actualP,
-                double *savr, double *wk, double *cj, double *wght, double *wp,
-                int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
+void C2F(psol1) (int *neq, double *tOld, double *actual, double *actualP,
+                 double *savr, double *wk, double *cj, double *wght, double *wp,
+                 int *iwp, double *b, double *eplin, int *ier, double *dummy1, int *dummy2)
 {
-    int i = 0;
-    int nColB = 1;
-    int info = 0;
-    int *ipiv = NULL;
+    int i = 0, job = 0;
 
-    ipiv = (int *) malloc(*nequations * sizeof(int));
+    C2F(dgesl) (wp, neq, neq, iwp, b, &job);
 
-    C2F(dgesv) (nequations, &nColB, wp, nequations, ipiv, b, nequations, &info);
-
-    if (info == 0)
+    /* NaN test */
+    for (i = 0; i < *neq; ++i)
     {
-        /* NaN test */
-        for (i = 0; i < *nequations; ++i)
+        if (ISNAN(b[i]))
         {
-            if (ISNAN(b[i]))
-            {
-                *ier = -1;
-                free (ipiv);
-                return;
-            }
+            /* Indicate a recoverable error, meaning that the step will be retried with the same step size,
+               but with a call to 'jacpsol' to update necessary data, unless the Jacobian data is current,
+               if (b[i] - b[i] != 0) in which case the step will be retried with a smaller step size. */
+            *ier = -1;
+            return;
         }
     }
-    else
-    {
-        *ier = -1;
-    }
-
-    free (ipiv);
 }
index 10fef2a..284a5ef 100644 (file)
@@ -59,7 +59,7 @@ function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
         ypsave = ydot(i);
         y(i) = y(i) + del;
         ydot(i) = ydot(i) + cj*del;
-        [e ires]=res1(tx, y, ydot);
+        [e ires] = res1(tx, y, ydot);
         if ires < 0 then
             ires = -1;
             return;
index c11ea98..35c6696 100644 (file)
@@ -63,7 +63,7 @@ function [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)
         ypsave = ydot(i);
         y(i) = y(i) + del;
         ydot(i) = ydot(i) + cj*del;
-        [e ires]=res1(tx, y, ydot);
+        [e ires] = res1(tx, y, ydot);
         if ires < 0 then
             ires = -1;
             return;