ode corrected 94/15094/7
Cedric Delamarre [Mon, 25 Aug 2014 13:50:40 +0000 (15:50 +0200)]
test_run differential_equations odedc
test_run differential_equations matode
test_run differential_equations ode

Change-Id: I75f1fd73aadfac541749d037256191c27e503489

20 files changed:
scilab/modules/differential_equations/Makefile.am
scilab/modules/differential_equations/Makefile.in
scilab/modules/differential_equations/sci_gateway/cpp/sci_ode.cpp
scilab/modules/differential_equations/sci_gateway/fortran/Ex-ode.f [deleted file]
scilab/modules/differential_equations/src/c/Ex-daskr.c
scilab/modules/differential_equations/src/c/Ex-daskr.h
scilab/modules/differential_equations/src/c/Ex-ode.c [new file with mode: 0644]
scilab/modules/differential_equations/src/c/Ex-ode.h [new file with mode: 0644]
scilab/modules/differential_equations/src/c/Ex-odedc.h
scilab/modules/differential_equations/src/c/differential_equations.vcxproj
scilab/modules/differential_equations/src/c/differential_equations.vcxproj.filters
scilab/modules/differential_equations/src/c/elementary_functions_f_Import.def [new file with mode: 0644]
scilab/modules/differential_equations/src/cpp/differentialequationfunctions.cpp
scilab/modules/differential_equations/src/fortran/differential_equations_f.vfproj
scilab/modules/differential_equations/tests/unit_tests/matode.dia.ref
scilab/modules/differential_equations/tests/unit_tests/matode.tst
scilab/modules/differential_equations/tests/unit_tests/ode.tst
scilab/modules/differential_equations/tests/unit_tests/ode.unix.dia.ref
scilab/modules/differential_equations/tests/unit_tests/ode.win.dia.ref
scilab/modules/differential_equations/tests/unit_tests/odedc.dia.ref

index 36f3505..2673247 100644 (file)
@@ -15,6 +15,7 @@ DIFFERENTIAL_EQUATIONS_C_SOURCES = \
     src/c/arnol.c \
     src/c/rk4.c \
     src/c/Ex-odedc.c \
+    src/c/Ex-ode.c \
     src/c/Ex-daskr.c
 
 DIFFERENTIAL_EQUATIONS_FORTRAN_SOURCES = \
@@ -78,7 +79,6 @@ DIFFERENTIAL_EQUATIONS_FORTRAN_SOURCES = \
     src/fortran/writbufode.f
 
 GATEWAY_FORTRAN_SOURCES = \
-    sci_gateway/fortran/Ex-ode.f \
     sci_gateway/fortran/Ex-impl.f \
     sci_gateway/fortran/Ex-int2d.f \
     sci_gateway/fortran/Ex-int3d.f \
index a2df40d..d09c9a9 100644 (file)
@@ -168,6 +168,7 @@ am__objects_1 = src/c/libscidifferential_equations_algo_la-dassl.lo \
        src/c/libscidifferential_equations_algo_la-arnol.lo \
        src/c/libscidifferential_equations_algo_la-rk4.lo \
        src/c/libscidifferential_equations_algo_la-Ex-odedc.lo \
+       src/c/libscidifferential_equations_algo_la-Ex-ode.lo \
        src/c/libscidifferential_equations_algo_la-Ex-daskr.lo
 am__objects_2 = src/fortran/rscar1.lo src/fortran/bcomp.lo \
        src/fortran/lcomp.lo src/fortran/loren.lo src/fortran/prja.lo \
@@ -212,8 +213,8 @@ am__v_lt_1 =
 @MAINTAINER_MODE_TRUE@am_libscidifferential_equations_algo_la_rpath =
 libscidifferential_equations_la_DEPENDENCIES =  \
        libscidifferential_equations-algo.la
-am__objects_4 = sci_gateway/fortran/Ex-ode.lo \
-       sci_gateway/fortran/Ex-impl.lo sci_gateway/fortran/Ex-int2d.lo \
+am__objects_4 = sci_gateway/fortran/Ex-impl.lo \
+       sci_gateway/fortran/Ex-int2d.lo \
        sci_gateway/fortran/Ex-int3d.lo sci_gateway/fortran/Ex-intg.lo \
        sci_gateway/fortran/Ex-dasrt.lo \
        sci_gateway/fortran/Ex-dassl.lo \
@@ -630,6 +631,7 @@ DIFFERENTIAL_EQUATIONS_C_SOURCES = \
     src/c/arnol.c \
     src/c/rk4.c \
     src/c/Ex-odedc.c \
+    src/c/Ex-ode.c \
     src/c/Ex-daskr.c
 
 DIFFERENTIAL_EQUATIONS_FORTRAN_SOURCES = \
@@ -693,7 +695,6 @@ DIFFERENTIAL_EQUATIONS_FORTRAN_SOURCES = \
     src/fortran/writbufode.f
 
 GATEWAY_FORTRAN_SOURCES = \
-    sci_gateway/fortran/Ex-ode.f \
     sci_gateway/fortran/Ex-impl.f \
     sci_gateway/fortran/Ex-int2d.f \
     sci_gateway/fortran/Ex-int3d.f \
@@ -935,6 +936,8 @@ src/c/libscidifferential_equations_algo_la-rk4.lo:  \
        src/c/$(am__dirstamp) src/c/$(DEPDIR)/$(am__dirstamp)
 src/c/libscidifferential_equations_algo_la-Ex-odedc.lo:  \
        src/c/$(am__dirstamp) src/c/$(DEPDIR)/$(am__dirstamp)
+src/c/libscidifferential_equations_algo_la-Ex-ode.lo:  \
+       src/c/$(am__dirstamp) src/c/$(DEPDIR)/$(am__dirstamp)
 src/c/libscidifferential_equations_algo_la-Ex-daskr.lo:  \
        src/c/$(am__dirstamp) src/c/$(DEPDIR)/$(am__dirstamp)
 src/fortran/$(am__dirstamp):
@@ -1080,8 +1083,6 @@ sci_gateway/fortran/$(am__dirstamp):
 sci_gateway/fortran/$(DEPDIR)/$(am__dirstamp):
        @$(MKDIR_P) sci_gateway/fortran/$(DEPDIR)
        @: > sci_gateway/fortran/$(DEPDIR)/$(am__dirstamp)
-sci_gateway/fortran/Ex-ode.lo: sci_gateway/fortran/$(am__dirstamp) \
-       sci_gateway/fortran/$(DEPDIR)/$(am__dirstamp)
 sci_gateway/fortran/Ex-impl.lo: sci_gateway/fortran/$(am__dirstamp) \
        sci_gateway/fortran/$(DEPDIR)/$(am__dirstamp)
 sci_gateway/fortran/Ex-int2d.lo: sci_gateway/fortran/$(am__dirstamp) \
@@ -1171,6 +1172,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscidifferential_equations_la-sci_ode.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscidifferential_equations_la-sci_odedc.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-daskr.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-ode.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-odedc.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscidifferential_equations_algo_la-arnol.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscidifferential_equations_algo_la-dassl.Plo@am__quote@
@@ -1239,6 +1241,13 @@ src/c/libscidifferential_equations_algo_la-Ex-odedc.lo: src/c/Ex-odedc.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscidifferential_equations_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o src/c/libscidifferential_equations_algo_la-Ex-odedc.lo `test -f 'src/c/Ex-odedc.c' || echo '$(srcdir)/'`src/c/Ex-odedc.c
 
+src/c/libscidifferential_equations_algo_la-Ex-ode.lo: src/c/Ex-ode.c
+@am__fastdepCC_TRUE@   $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscidifferential_equations_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT src/c/libscidifferential_equations_algo_la-Ex-ode.lo -MD -MP -MF src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-ode.Tpo -c -o src/c/libscidifferential_equations_algo_la-Ex-ode.lo `test -f 'src/c/Ex-ode.c' || echo '$(srcdir)/'`src/c/Ex-ode.c
+@am__fastdepCC_TRUE@   $(AM_V_at)$(am__mv) src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-ode.Tpo src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-ode.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      $(AM_V_CC)source='src/c/Ex-ode.c' object='src/c/libscidifferential_equations_algo_la-Ex-ode.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscidifferential_equations_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o src/c/libscidifferential_equations_algo_la-Ex-ode.lo `test -f 'src/c/Ex-ode.c' || echo '$(srcdir)/'`src/c/Ex-ode.c
+
 src/c/libscidifferential_equations_algo_la-Ex-daskr.lo: src/c/Ex-daskr.c
 @am__fastdepCC_TRUE@   $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscidifferential_equations_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT src/c/libscidifferential_equations_algo_la-Ex-daskr.lo -MD -MP -MF src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-daskr.Tpo -c -o src/c/libscidifferential_equations_algo_la-Ex-daskr.lo `test -f 'src/c/Ex-daskr.c' || echo '$(srcdir)/'`src/c/Ex-daskr.c
 @am__fastdepCC_TRUE@   $(AM_V_at)$(am__mv) src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-daskr.Tpo src/c/$(DEPDIR)/libscidifferential_equations_algo_la-Ex-daskr.Plo
index ccf4a6f..c4147ec 100644 (file)
@@ -39,7 +39,7 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
 {
     // Methode
     types::String* pStrType     = NULL;
-    const wchar_t* wcsType            = L"lsoda";
+    const wchar_t* wcsType      = L"lsoda";
     int meth                    = 0;
 
     // y0
@@ -938,13 +938,35 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 rwork[4] = pDblOdeOptions->get(2);      // h0 =0
                 rwork[5] = pDblOdeOptions->get(3);      // hmax = %inf
                 rwork[6] = pDblOdeOptions->get(4);      // hmin = 0
-                iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
-                iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
-                iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
+
+                if (jt == 4 || jt == 5)
+                {
+                    iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
+                    iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
+                }
+
+                if (meth == 0 || meth == 3)
+                {
+                    // lsoda/lsodar
+                    iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
+                    iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
+                    iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
+                }
+                else // 1 or 2
+                {
+                    // lsode
+                    if (meth == 1)
+                    {
+                        iwork[4] = (int)pDblOdeOptions->get(7); // mxordn = 12
+                    }
+                    else // meth == 2
+                    {
+                        iwork[4] = (int)pDblOdeOptions->get(8); // mxords = 5
+                    }
+                }
+
                 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
                 iwork[6] = 10;  // mxhnil = 10 maximum number of messages printed per problem
-                iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
-                iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
             }
         }
     }
@@ -1125,18 +1147,21 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 break;
             }
 
-            if (meth == 3 && istate == 3 && getWarningMode())
+            if (meth == 3 && istate == 3)
             {
                 // istate == 3  means the integration was successful, and one or more
                 //              roots were found before satisfying the stop condition
                 //              specified by itask.  see jroot.
-
-                sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
-                for (int k = 0; k < pDblNg->get(0); k++)
+                if (getWarningMode())
                 {
-                    sciprint("\t%d", jroot[k]);
+                    sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
+                    for (int k = 0; k < pDblNg->get(0); k++)
+                    {
+                        sciprint("\t%d", jroot[k]);
+                    }
+                    sciprint("\n");
                 }
-                sciprint("\n");
+
                 break;
             }
         }
@@ -1213,7 +1238,7 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 else
                 {
         */
-        pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize());
+        pDblYOut = new types::Double(pDblY0->getRows(), pDblY0->getCols() * pDblT->getSize());
         //        }
         bool bBreak = false;
         for (int i = 0; i < pDblT->getSize(); i++)
@@ -1359,18 +1384,21 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 break;
             }
 
-            if (meth == 3 && istate == 3 && getWarningMode())
+            if (meth == 3 && istate == 3)
             {
                 // istate == 3  means the integration was successful, and one or more
                 //              roots were found before satisfying the stop condition
                 //              specified by itask.  see jroot.
 
-                sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
-                for (int k = 0; k < pDblNg->get(0); k++)
+                if (getWarningMode())
                 {
-                    sciprint("\t%d", jroot[k]);
+                    sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
+                    for (int k = 0; k < pDblNg->get(0); k++)
+                    {
+                        sciprint("\t%d", jroot[k]);
+                    }
+                    sciprint("\n");
                 }
-                sciprint("\n");
 
                 types::Double* pDblYOutTemp = pDblYOut;
                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
diff --git a/scilab/modules/differential_equations/sci_gateway/fortran/Ex-ode.f b/scilab/modules/differential_equations/sci_gateway/fortran/Ex-ode.f
deleted file mode 100644 (file)
index bad7b6d..0000000
+++ /dev/null
@@ -1,168 +0,0 @@
-c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-c Copyright (C) INRIA
-c ...
-c 
-c This file must be used under the terms of the CeCILL.
-c This source file is licensed as described in the file COPYING, which
-c you should have received as part of this distribution.  The terms
-c are also available at    
-c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
-c
-c-----------------------------------------------------------------
-c     Examples of argument functions for scilab function ode 
-c     The function provided here must be listed 
-c     in the file Flist in order to be used with scilab 
-c     fydot_list : contains the rhs function
-c     fjac_list   : contains the jacobian (if you want to give it)
-c                   ( example jex -> jacobian for fex)
-c     Examples :
-c       fex, jex : first example 
-c       fex2     : uses creadmat (or matz) to get parameters 
-c       fex3     : uses matptr to get parameters 
-c       fexab    : a more complex example with matpr
-c-----------------------------------------------------------------
-c
-      subroutine fex (neq, t, y, ydot)
-c     -------------------------------------------
-c     input variables neq, t, y
-c     neq = dimension of state vector y
-c     t = time
-c     y = state variable 
-c     output variable = ydot
-c
-c     argument routine must 
-c     load ydot(1) with d/dt ( y(1)(t) )
-c          ydot(2) with d/dt ( y(2)(t) )
-c          ...
-c     i.e. ydot vector of derivative of state y
-c
-c     Example:
-c     1/ call this fex routine: 
-c     ode([1;0;0],0,[0.4,4],'fex')
-c
-c     2/ May use dynamic link: 
-c     link('fex.o','fex')
-c     ode([1;0;0],0,[0.4,4],'fex')
-c
-c     3/ passing a parameter to fex routine:
-c     link('fex.o','fex')    
-c     y=ode([1;0;0],0,[0.4,4],list('fex',param))
-c     param can be retrieved in fex by:
-c     param(1)=y(neq+1) , param(2)=y(neq+2) etc 
-c     with this calling sequence y is a nc+nd+np vector
-c     where np=dimension of scilab variable param
-c
-      double precision t, y, ydot
-      dimension y(3), ydot(3)
-      ydot(1) = -.0400d+0*y(1) + 1.0d+4*y(2)*y(3)
-      ydot(3) = 3.0d+7*y(2)*y(2)
-      ydot(2) = -ydot(1) - ydot(3)
-      return
-      end
-
-
-      subroutine jex (neq, t, y, ml, mu, pd, nrpd)
-c     -------------------------------------------
-c     fex example continued : we provide here 
-c     the jacobian function 
-c     jacobian routine jex
-c     scilab ode
-c     ode([1;0;0],0,[0.4,4],'fex','jex')
-
-      double precision pd, t, y
-      dimension y(3), pd(nrpd,3)
-      pd(1,1) = -.040d+0
-      pd(1,2) = 1.0d+4*y(3)
-      pd(1,3) = 1.0d+4*y(2)
-      pd(2,1) = .040d+0
-      pd(2,3) = -pd(1,3)
-      pd(3,2) = 6.0d+7*y(2)
-      pd(2,2) = -pd(1,2) - pd(3,2)
-      return
-      end
-
-
-!      subroutine fex2(neq, t, y, ydot)
-!c     -------------------------------------------
-!c     exemple with a call to creadmat routine
-!c     param must be defined as a scilab variable
-!      double precision t, y, ydot, param
-!      dimension y(3), ydot(3), param(3)
-!
-!      call creadmat('param'//char(0),m,n,param)
-!c     *******************************
-!
-!      ydot(1) = param(1)*y(1) + param(2)*y(2)*y(3)
-!      ydot(3) = param(3)*y(2)*y(2)
-!      ydot(2) = -ydot(1) - ydot(3)
-!      return
-!      end
-
-
-!      subroutine fex3(neq, t, y, ydot)
-!c     -------------------------------------------
-!c     same example with call to matptr
-!c     param must be defined as a scilab variable
-!      double precision t, y, ydot
-!      dimension y(3), ydot(3)
-!
-!      include 'stack.h'
-!      call matptr('param'//char(0),m,n,lp)
-!      if(m.eq.-1) return
-!c     ***************************
-!c     param entries are in stk(lp),stk(lp+1),stk(lp+2)
-!c     m,n = dimensions of param = 3,1 (or 1,3 if row v.)
-!c     (note that vector param not used in this example)
-!      ydot(1) = stk(lp)*y(1) + stk(lp+1)*y(2)*y(3)
-!      ydot(3) = stk(lp+2)*y(2)*y(2)
-!      ydot(2) = -ydot(1) - ydot(3)
-!      return
-!      end
-
-
-!      subroutine fexab(neq, t, x, xdot)
-!c     -----------------------------------
-!c    Another example
-!c     xdot=A*xc+B*u
-!c     A and B real scilab matrices
-!c     u=u(t,x) given by function uinput below
-!c
-!c     [ User may extract this fexab.f routine 
-!c       and -->link('fexab.o','fexab') ]
-!c
-!c     -->A=..., B=... 
-!c     -->ode(x0,t0,t,'fexab')
-!c     
-!c     or customize this one and re-make scilab
-!c
-!c     input variables neq, t, x
-!c     neq = dimension of state vector x
-!c     t = time
-!c     x = state variable 
-!c     output variable = xdot
-!c
-!c     routine must 
-!c     load xdot(1) with d/dt ( x(1)(t) )
-!c          xdot(2) with d/dt ( x(2)(t) )
-!c          ...
-!c     i.e. xdot vector of derivative of state x
-!c
-!c
-!      double precision t,x(*),xdot(*)
-!c     set dimension of u here
-!      double precision u(1)
-!c
-!      include 'stack.h'
-!      call matptr('A'//char(0),m,n,la)
-!      call dmmul(stk(la),m,x,m,xdot,m,m,m,1)
-!      call matptr('B'//char(0),m,nb,lb)
-!      call uinput(t,x,u)
-!      call dgemm('n','n',m,1,nb,1.0d0,stk(lb),m,u,1,1.0d0,xdot,m)
-!      end
-!
-!      subroutine uinput(t,x,u)
-!      double precision t,x(*),u(*)
-!      u(1)=sin(3*t)
-!c     u(2)=...
-!      end
-
index a9de3f3..a242c41 100644 (file)
@@ -3,6 +3,7 @@
 #include "machine.h"
 #include "core_math.h"
 #include "Ex-daskr.h"
+#include "elem_common.h" // dlamch
 
 void 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,
index 512f147..5e4de8a 100644 (file)
@@ -11,7 +11,8 @@
 */
 /*--------------------------------------------------------------------------*/
 
-#include "elem_common.h" // dlamch
+#ifndef __EX_DASKR_H__
+#define __EX_DASKR_H__
 
 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);
@@ -26,3 +27,4 @@ void 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);
 
+#endif // __EX_DASKR_H__
diff --git a/scilab/modules/differential_equations/src/c/Ex-ode.c b/scilab/modules/differential_equations/src/c/Ex-ode.c
new file mode 100644 (file)
index 0000000..c9acfb0
--- /dev/null
@@ -0,0 +1,190 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) INRIA
+ * Copyright (C) 2014 - Scilab Enterprises - Cedric DELAMARRE
+ *
+ * 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
+ *
+ */
+/*--------------------------------------------------------------------------*/
+#include "Ex-ode.h"
+#include "api_scilab.h"
+#include "elem_common.h"
+
+//-----------------------------------------------------------------
+//     Examples of argument functions for scilab function ode
+//     The function provided here must be listed
+//     in the file Flist in order to be used with scilab
+//     fydot_list : contains the rhs function
+//     fjac_list   : contains the jacobian (if you want to give it)
+//                   ( example jex -> jacobian for fex)
+//     Examples :
+//       fex, jex : first example
+//       fex2     : uses creadmat (or matz) to get parameters
+//       fex3     : uses matptr to get parameters
+//       fexab    : a more complex example with matpr
+//-----------------------------------------------------------------
+
+extern C2F(dmmul)(double*, int*, double*, int*, double*, int*, int*, int*, int*);
+
+void fex(int* neq, double* t, double* y, double* ydot)
+{
+    //        input variables neq, t, y
+    //        neq = dimension of state vector y
+    //        t = time
+    //        y = state variable
+    //        output variable = ydot
+    //
+    //        argument routine must
+    //        load ydot(1) with d/dt ( y(1)(t) )
+    //          ydot(2) with d/dt ( y(2)(t) )
+    //          ...
+    //        i.e. ydot vector of derivative of state y
+    //
+    //        Example:
+    //        1/ call this fex routine:
+    //        ode([1;0;0],0,[0.4,4],'fex')
+    //
+    //        2/ May use dynamic link:
+    //        link('fex.o','fex')
+    //        ode([1;0;0],0,[0.4,4],'fex')
+    //
+    //        3/ passing a parameter to fex routine:
+    //        link('fex.o','fex')
+    //        y=ode([1;0;0],0,[0.4,4],list('fex',param))
+    //        param can be retrieved in fex by:
+    //        param(1)=y(neq+1) , param(2)=y(neq+2) etc
+    //        with this calling sequence y is a nc+nd+np vector
+    //        where np=dimension of scilab variable param
+
+    ydot[0] = -0.0400 * y[0] + 1e4 * y[1] * y[2];
+    ydot[2] = 3e7 * y[1] * y[1];
+    ydot[1] = -ydot[0] - ydot[2];
+}
+
+void jex(int* neq, double* t, double* y, double* ml, double* mu, double* pd, int* nrpd)
+{
+    //        fex example continued : we provide here
+    //        the jacobian function
+    //        jacobian routine jex
+    //        scilab ode
+    //        ode([1;0;0],0,[0.4,4],'fex','jex')
+    pd[0] = -0.040;
+    pd[3] = 1e4 * y[2];
+    pd[6] = 1e4 * y[1];
+    pd[1] = 0.040;
+    pd[7] = -pd[6];
+    pd[5] = 6e7 * y[1];
+    pd[4] = -pd[3] - pd[5];
+}
+
+
+void fex2(int* neq, double* t, double* y, double* ydot)
+{
+    //     exemple with a call to creadmat routine
+    //     param must be defined as a scilab variable
+
+    // pvApiCtx not used when we get a named variable
+    void* pvApiCtx = NULL;
+    int iRows = 0;
+    int iCols = 0;
+    double* pdblData = NULL;
+
+    readNamedMatrixOfDouble(pvApiCtx, "param", &iRows, &iCols, pdblData);
+    pdblData = (double*)malloc(iRows * iCols * sizeof(double));
+    readNamedMatrixOfDouble(pvApiCtx, "param", &iRows, &iCols, pdblData);
+
+    ydot[0] = pdblData[0] * y[0] + pdblData[1] * y[1] * y[2];
+    ydot[2] = pdblData[2] * y[1] * y[1];
+    ydot[1] = -ydot[0] - ydot[2];
+
+    free(pdblData);
+}
+
+void fex3(int* neq, double* t, double* y, double* ydot)
+{
+    //     same example with call to matptr
+    //     param must be defined as a scilab variable
+
+    // pvApiCtx not used when we get a named variable
+    void* pvApiCtx = NULL;
+    int iRows = 0;
+    int iCols = 0;
+    double* pdblData = NULL;
+
+    readNamedMatrixOfDouble(pvApiCtx, "param", &iRows, &iCols, pdblData);
+    pdblData = (double*)malloc(iRows * iCols * sizeof(double));
+    readNamedMatrixOfDouble(pvApiCtx, "param", &iRows, &iCols, pdblData);
+
+    if (iRows == -1)
+    {
+        return;
+    }
+
+    //     param entries are in stk(lp),stk(lp+1),stk(lp+2)
+    //     m,n = dimensions of param = 3,1 (or 1,3 if row v.)
+    //     (note that vector param not used in this example)
+    ydot[0] = pdblData[0] * y[0] + pdblData[1] * y[1] * y[2];
+    ydot[2] = pdblData[2] * y[1] * y[1];
+    ydot[1] = -ydot[0] - ydot[2];
+
+    free(pdblData);
+}
+
+void fexab(int* neq, double* t, double* x, double* xdot)
+{
+    //    Another example
+    //     xdot=A*xc+B*u
+    //     A and B real scilab matrices
+    //     u=u(t,x)
+    //
+    //     [ User may extract this fexab.f routine
+    //       and -->link('fexab.o','fexab') ]
+    //
+    //     -->A=..., B=...
+    //     -->ode(x0,t0,t,'fexab')
+    //
+    //     or customize this one and re-make scilab
+    //
+    //     input variables neq, t, x
+    //     neq = dimension of state vector x
+    //     t = time
+    //     x = state variable
+    //     output variable = xdot
+    //
+    //     routine must
+    //     load xdot(1) with d/dt ( x(1)(t) )
+    //          xdot(2) with d/dt ( x(2)(t) )
+    //          ...
+    //     i.e. xdot vector of derivative of state x
+
+    void* pvApiCtx = NULL;
+    char cN = 'n';
+    int iOne = 1;
+    int iRows = 0;
+    int iColsA = 0;
+    int iColsB = 0;
+    double pdblOne = 1.0;
+    double pdblU = 0;
+    double* pdblDataA = NULL;
+    double* pdblDataB = NULL;
+
+    readNamedMatrixOfDouble(pvApiCtx, "A", &iRows, &iColsA, pdblDataA);
+    pdblDataA = (double*)malloc(iRows * iColsA * sizeof(double));
+    readNamedMatrixOfDouble(pvApiCtx, "A", &iRows, &iColsA, pdblDataA);
+    C2F(dmmul)(pdblDataA, &iRows, x, &iRows, xdot, &iRows, &iRows, &iRows, &iOne);
+
+    readNamedMatrixOfDouble(pvApiCtx, "B", &iRows, &iColsB, pdblDataB);
+    pdblDataB = (double*)malloc(iRows * iColsB * sizeof(double));
+    readNamedMatrixOfDouble(pvApiCtx, "B", &iRows, &iColsB, pdblDataB);
+
+    pdblU = sin(3 * t[0]);
+    C2F(dgemm)(&cN, &cN, &iRows, &iOne, &iColsB, &pdblOne, pdblDataB, &iRows, &pdblU, &iOne, &pdblOne, xdot, &iRows);
+
+    free(pdblDataA);
+    free(pdblDataB);
+}
diff --git a/scilab/modules/differential_equations/src/c/Ex-ode.h b/scilab/modules/differential_equations/src/c/Ex-ode.h
new file mode 100644 (file)
index 0000000..36b183c
--- /dev/null
@@ -0,0 +1,23 @@
+/*
+* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+* Copyright (C) 2014 - Scilab Enterprises - Cedric DELAMARRE
+*
+* 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
+*
+*/
+/*--------------------------------------------------------------------------*/
+
+#ifndef __EX_ODE_H__
+#define __EX_ODE_H__
+
+void fex(int* neq, double* t, double* y, double* ydot);
+void jex(int* neq, double* t, double* y, double* ml, double* mu, double* pd, int* nrpd);
+void fex2(int* neq, double* t, double* y, double* ydot);
+void fex3(int* neq, double* t, double* y, double* ydot);
+void fexab(int* neq, double* t, double* x, double* xdot);
+
+#endif // __EX_ODE_H__
index 37a1471..ce46896 100644 (file)
@@ -11,6 +11,9 @@
 */
 /*--------------------------------------------------------------------------*/
 
+#ifndef __EX_ODEDC_H__
+#define __EX_ODEDC_H__
+
 void fexcd(int* jflag, int* nc, int* nd, double* t, double* y, double* ydp);
 void fcd(int* jflag, int* nc, int* nd, double* t, double* x, double* xdp);
 void fc(double t, double* xc, double* u, double* xdot);
@@ -27,3 +30,4 @@ void phis(int* jflag, int* nc, int* nd, double* t, double* x, double* xdp);
 void phit(int* jflag, int* nc, int* nd, double* t, double* x, double* xdp);
 void sbrc(double t, double* x, double* xdot);
 
+#endif // __EX_ODEDC_H__
index eeb8336..9f9ff41 100644 (file)
       <Message>Make dependencies</Message>
       <Command>lib /DEF:"$(ProjectDir)core_import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)differential_equations_f.lib" 1&gt;NUL 2&gt;NUL
-</Command>
+lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL</Command>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;elementary_functions_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <SubSystem>Windows</SubSystem>
@@ -117,10 +117,10 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
       <Message>Make dependencies</Message>
       <Command>lib /DEF:"$(ProjectDir)core_import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)differential_equations_f.lib" 1&gt;NUL 2&gt;NUL
-</Command>
+lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL</Command>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;elementary_functions_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <SubSystem>Windows</SubSystem>
@@ -145,10 +145,10 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
       <Message>Make dependencies</Message>
       <Command>lib /DEF:"$(ProjectDir)core_import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)differential_equations_f.lib" 1&gt;NUL 2&gt;NUL
-</Command>
+lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL</Command>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;elementary_functions_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <GenerateDebugInformation>false</GenerateDebugInformation>
       <SubSystem>Windows</SubSystem>
@@ -178,10 +178,10 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
       <Message>Make dependencies</Message>
       <Command>lib /DEF:"$(ProjectDir)core_import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)differential_equations_f.lib" 1&gt;NUL 2&gt;NUL
-</Command>
+lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL</Command>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/linpack_f.lib;core.lib;differential_equations_f.lib;elementary_functions_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <GenerateDebugInformation>false</GenerateDebugInformation>
       <SubSystem>Windows</SubSystem>
@@ -201,6 +201,7 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
     <ClCompile Include="dassl.c" />
     <ClCompile Include="DllmainDifferential_equations.c" />
     <ClCompile Include="Ex-daskr.c" />
+    <ClCompile Include="Ex-ode.c" />
     <ClCompile Include="Ex-odedc.c" />
     <ClCompile Include="feval.c" />
     <ClCompile Include="rk4.c" />
@@ -213,6 +214,7 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
     <ClInclude Include="arnol.h" />
     <ClInclude Include="..\..\includes\dynlib_differential_equations.h" />
     <ClInclude Include="Ex-daskr.h" />
+    <ClInclude Include="Ex-ode.h" />
     <ClInclude Include="Ex-odedc.h" />
     <ClInclude Include="feval.h" />
     <ClInclude Include="rk4.h" />
@@ -224,6 +226,7 @@ lib /DEF:"$(ProjectDir)differential_equations_f_Import.def" /SUBSYSTEM:WINDOWS /
     <None Include="..\..\differential_equations.iss" />
     <None Include="..\..\sci_gateway\differential_equations_gateway.xml" />
     <None Include="..\..\Makefile.am" />
+    <None Include="elementary_functions_f_Import.def" />
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="differential_equations.rc" />
index 18eccda..67ff839 100644 (file)
@@ -53,6 +53,9 @@
     <ClCompile Include="Ex-daskr.c">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="Ex-ode.c">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="arnol.h">
@@ -85,6 +88,9 @@
     <ClInclude Include="Ex-daskr.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="Ex-ode.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <None Include="differential_equations_f_Import.def">
     <None Include="..\..\locales\differential_equations.pot">
       <Filter>localization</Filter>
     </None>
+    <None Include="elementary_functions_f_Import.def">
+      <Filter>Libraries Dependencies\Imports</Filter>
+    </None>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="differential_equations.rc">
diff --git a/scilab/modules/differential_equations/src/c/elementary_functions_f_Import.def b/scilab/modules/differential_equations/src/c/elementary_functions_f_Import.def
new file mode 100644 (file)
index 0000000..424672e
--- /dev/null
@@ -0,0 +1,8 @@
+LIBRARY    elementary_functions_f.dll
+
+
+EXPORTS
+; ---------------------------------------
+;elementary_functions_f
+; ---------------------------------------
+dmmul_
\ No newline at end of file
index 72f265d..41456ab 100644 (file)
@@ -20,6 +20,7 @@ extern "C"
 #include "elem_common.h"
 #include "scifunctions.h"
 #include "Ex-odedc.h"
+#include "Ex-ode.h"
 #include "Ex-daskr.h"
 #include "localization.h"
 }
@@ -117,15 +118,15 @@ DifferentialEquationFunctions::DifferentialEquationFunctions(std::wstring caller
     if (callerName == L"ode")
     {
         m_staticFunctionMap[L"arnol"]   = (void*) C2F(arnol);
-        m_staticFunctionMap[L"fex"]     = (void*) C2F(fex);
-        //m_staticFunctionMap[L"fex2"]    = (void*) C2F(fex2);
-        //m_staticFunctionMap[L"fex3"]    = (void*) C2F(fex3);
-        //m_staticFunctionMap[L"fexab"]   = (void*) C2F(fexab);
+        m_staticFunctionMap[L"fex"]     = (void*) fex;
+        m_staticFunctionMap[L"fex2"]    = (void*) fex2;
+        m_staticFunctionMap[L"fex3"]    = (void*) fex3;
+        m_staticFunctionMap[L"fexab"]   = (void*) fexab;
         m_staticFunctionMap[L"loren"]   = (void*) C2F(loren);
         m_staticFunctionMap[L"bcomp"]   = (void*) C2F(bcomp);
         m_staticFunctionMap[L"lcomp"]   = (void*) C2F(lcomp);
 
-        m_staticFunctionMap[L"jex"]     = (void*) C2F(jex);
+        m_staticFunctionMap[L"jex"]     = (void*) jex;
     }
     else if (callerName == L"odedc")
     {
@@ -135,7 +136,7 @@ DifferentialEquationFunctions::DifferentialEquationFunctions(std::wstring caller
         m_staticFunctionMap[L"phis"]    = (void*) phis;
         m_staticFunctionMap[L"phit"]    = (void*) phit;
 
-        m_staticFunctionMap[L"jex"]     = (void*) C2F(jex);
+        m_staticFunctionMap[L"jex"]     = (void*) jex;
     }
     else if (callerName == L"intg")
     {
index c189863..e5a2af3 100644 (file)
@@ -77,6 +77,8 @@
                <File RelativePath=".\dqags.f"/>
                <File RelativePath=".\dqagse.f"/>
                <File RelativePath=".\dqelg.f"/>
+               <File RelativePath=".\dqk21.f"/>
+               <File RelativePath=".\dqpsrt.f"/>
                <File RelativePath=".\ewset.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-bvode.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-dasrt.f"/>
@@ -85,7 +87,6 @@
                <File RelativePath="..\..\sci_gateway\fortran\Ex-int2d.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-int3d.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-intg.f"/>
-               <File RelativePath="..\..\sci_gateway\fortran\Ex-ode.f"/>
                <File RelativePath=".\fnorm.f"/>
                <File RelativePath=".\greatr.f"/>
                <File RelativePath=".\hpdel.f"/>
                <File RelativePath=".\lsodi.f"/>
                <File RelativePath=".\lsrgk.f"/>
                <File RelativePath=".\odeint.f"/>
-               <File RelativePath=".\dqpsrt.f"/>
                <File RelativePath=".\prepj.f"/>
                <File RelativePath=".\prepji.f"/>
                <File RelativePath=".\prja.f"/>
-               <File RelativePath=".\dqk21.f"/>
                <File RelativePath=".\rchek.f"/>
                <File RelativePath=".\rchek2.f"/>
                <File RelativePath=".\rkf45.f"/>
index 3953508..5da75a7 100644 (file)
@@ -4,6 +4,7 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
+// <-- CLI SHELL MODE -->
 warning("off");
 funcprot(0);
 // Copyright INRIA
@@ -14,56 +15,68 @@ Leps=1.e-6;
 // on the interval from t = 0.0 to t = 4.e10, with initial conditions
 // y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
 //     definition of rhs
-deff('[ydot]=f(t,y)',...
-   'f1=-0.04*y(1)+1e4*y(2)*y(3),...
-    f3=3e7*y(2)*y(2),...
-    ydot=[f1;-f1-f3;f3]','n')
+function [ydot]=f(t,y)
+    f1=-0.04*y(1)+1e4*y(2)*y(3)
+    f3=3e7*y(2)*y(2)
+    ydot=[f1;-f1-f3;f3]
+endfunction
 //     jacobian
-deff('[jac]=j(t,y)',...
-'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
- jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
- jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);','n')
+function [jac]=j(t,y)
+    jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2)
+    jac(3,1)=0
+    jac(3,2)=6.e7*y(2)
+    jac(3,3)=0
+    jac(2,1)=-jac(1,1)
+    jac(2,2)=-jac(1,2)-jac(3,2)
+    jac(2,3)=-jac(1,3)
+endfunction
 //
-y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
+y0=[1;0;0];
+t0=0;
+t1=[0.4,4];
+nt=size(t1,"*");
 //    solution
 yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 //
 //  1. fortran called by fydot, without jacobian
-y1=ode(y0,t0,t1,'fex');
+y1=ode(y0,t0,t1,"fex");
 if max(y1-yref) > Leps then bugmes();quit;end
 //  2. fortran called by fydot, type given (stiff), no jacobian
-y2=ode('stiff',y0,t0,t1,'fex');
+y2=ode("stiff",y0,t0,t1,"fex");
 if max(y2-yref) > Leps then bugmes();quit;end
 //  3. fortran called by fydot , fjac, type given
-y3=ode('stiff',y0,t0,t1,'fex','jex');
+y3=ode("stiff",y0,t0,t1,"fex","jex");
 if max(y3-yref) > Leps then bugmes();quit;end
 //   hot restart
-[z,w,iw]=ode('stiff',y0,0,0.4,'fex','jex');
-z=ode('stiff',z,0.4,4,'fex','jex',w,iw);
+[z,w,iw]=ode("stiff",y0,0,0.4,"fex","jex");
+z=ode("stiff",z,0.4,4,"fex","jex",w,iw);
 if max(z-y3(:,2)) > %eps then bugmes();quit;end
-[y1,w,iw]=ode(y0,t0,t1(1),'fex');
-y2=ode(y0,t1(1),t1(2:nt),'fex',w,iw);
+[y1,w,iw]=ode(y0,t0,t1(1),"fex");
+y2=ode(y0,t1(1),t1(2:nt),"fex",w,iw);
 if max([y1 y2]-yref) > Leps then bugmes();quit;end
-[y1,w,iw]=ode(y0,t0,t1(1),'fex','jex');
-y2=ode(y0,t1(1),t1(2:nt),'fex','jex',w,iw);
+[y1,w,iw]=ode(y0,t0,t1(1),"fex","jex");
+y2=ode(y0,t1(1),t1(2:nt),"fex","jex",w,iw);
 if max([y1 y2]-yref) > Leps then bugmes();quit;end
 //   variation of tolerances
-atol=[0.001,0.0001,0.001];rtol=atol;
+atol=[0.001,0.0001,0.001];
+rtol=atol;
 //    externals
 //  4. type given , scilab lhs ,jacobian not passed
-y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
+y4=ode("stiff",y0,t0,t1(1),atol,rtol,f);
 if max(y4(:,1)-yref(:,1)) > 0.01 then bugmes();quit;end
 //  5. type non given, rhs and scilab jacobian
 y5=ode(y0,t0,t1,f,j);
 if max(y5-yref) > Leps then bugmes();quit;end
 //  6. type given (stiff),rhs and jacobian  by scilab
-y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
+y6=ode("stiff",y0,t0,t1,0.00001,0.00001,f,j);
 if (y6-yref) > 2*0.00001 then bugmes();quit;end
 //  7. matrix rhs, type given(adams),jacobian not passed
 //
 a=rand(3,3);ea=expm(a);
-deff('[ydot]=f(t,y)','ydot=a*y')
-t1=1;y=ode('adams',eye(a),t0,t1,f);
+function [ydot]=f(t,y)
+    ydot=a*y
+endfunction
+t1=1;y=ode("adams",eye(a),t0,t1,f);
 if max(ea-y) > Leps then bugmes();quit;end
 //
 //   DAE's
@@ -72,24 +85,27 @@ if max(ea-y) > Leps then bugmes();quit;end
 //       0.   = y1 + y2 + y3 - 1.
 //  y1(0) = 1.0, y2(0) = y3(0) = 0.
 // scilab macros
-deff('[r]=resid(t,y,s)','...
-r(1)=-0.04*y(1)+1.d4*y(2)*y(3)-s(1),...
-r(2)=0.04*y(1)-1.d4*y(2)*y(3)-3.d7*y(2)*y(2)-s(2),...
-r(3)=y(1)+y(2)+y(3)-1')
-deff('[p]=aplusp(t,y,p)','...
-p(1,1)=p(1,1)+1,...
-p(2,2)=p(2,2)+1')
-deff('[p]=dgbydy(t,y,s)','[-0.04,1.d4*y(3),1.d4*y(2);...
-0.04,-1.d4*y(3)-6.d7*y(2),-1.d4*y(2);...
-1,1,1]')
+function [r]=resid(t,y,s)
+    r(1)=-0.04*y(1)+1.d4*y(2)*y(3)-s(1)
+    r(2)=0.04*y(1)-1.d4*y(2)*y(3)-3.d7*y(2)*y(2)-s(2)
+    r(3)=y(1)+y(2)+y(3)-1
+endfunction
+function [p]=aplusp(t,y,p)
+    p(1,1)=p(1,1)+1
+    p(2,2)=p(2,2)+1
+endfunction
 //        %ODEOPTIONS tests
 //
-rand('seed',0);rand('normal');
+rand("seed",0);rand("normal");
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
-deff('y=f(t,x)','y=A*x')
-deff('J=j(t,x)','J=A')
+function y=f(t,x)
+    y=A*x
+endfunction
+function J=j(t,x)
+    J=A
+endfunction
 x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
-nt=size(t,'*');
+nt=size(t,"*");
 eAt=expm(A*t(nt));
 //        Test itask=%ODEOPTIONS(1)
 //itask=1  --->  usual call (t=vector of instants, solution at all t)
@@ -100,7 +116,6 @@ jacflag=2;ml=-1;mu=-1;
 xf=ode(x0,t0,t,f);   //lsoda
 if norm(xf(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
-
 if norm(xfj(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 //itask=2;   --->  solution at mesh points  ---> t=tmax
 //========================
@@ -118,12 +133,15 @@ if norm(x3(2:nx+1)-xlast,1)>Leps then bugmes();quit;end
 %ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
 tcrit=2.5;%ODEOPTIONS(2)=tcrit;
 chk=0;
-deff('y=fcrit(t,x)',['if t<=tcrit then'
-                      ' y=A*x;'
-                      'else'
-                      ' y=A*x;chk=resume(1);end'])
+function y=fcrit(t,x)
+    if t<=tcrit then
+        y=A*x;
+    else
+        y=A*x;
+        chk=resume(1);
+    end
+endfunction
 x42=ode(x0,t0,t,fcrit);
-
 if chk==1 then bugmes();quit;end
 [p,q]=size(x42);
 if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then bugmes();quit;end
@@ -144,25 +162,22 @@ if norm(x35-xf,1)>Leps then bugmes();quit;end
 //test of %ODEOPTIONS(6)=jacflag
 %ODEOPTIONS(6)=1;//Jacobian given
 %ODEOPTIONS(3:5)=[0 0 0];
-x61=ode('st',x0,t0,t,f,j);   //with Jacobian
+x61=ode("stiff",x0,t0,t,f,j);   //with Jacobian
 if norm(x61-xf,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=0; // jacobian nor called nor estimated
-x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
-
-x60=ode('st',x0,t0,t,f);    //Jacobian not used
+x60=ode("stiff",x0,t0,t,f,j);   //Jacobian not used (warning)
+x60=ode("stiff",x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=1;//Jacobian estimated
-x60=ode('st',x0,t0,t,f)  ;
-
+x60=ode("stiff",x0,t0,t,f);
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 //test of %ODEOPTIONS(6)=jacflag   (adams)
 %ODEOPTIONS(6)=1;//with given Jacobian
-x60=ode('ad',x0,t0,t,f,j) ;
+x60=ode("adams",x0,t0,t,f,j) ;
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=0;// jacobian nor called nor estimated
-x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
-
-x60=ode('ad',x0,t0,t,f);    //Jacobian not used
+x60=ode("adams",x0,t0,t,f,j);   //Jacobian not used (warning)
+x60=ode("adams",x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 // test lsoda
 %ODEOPTIONS(6)=2;// jacobian  estimated
@@ -170,7 +185,6 @@ x60=ode(x0,t0,t,f);
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 %ODEOPTIONS(6)=1;//Jacobian estimated
 x60=ode(x0,t0,t,f);
-
 if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 //   Banded Jacobian
 itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
@@ -185,12 +199,14 @@ for k=1:nx-2, A(k,k+2)=-2;end
 for k=1:nx-1, A(k+1,k)=-1;end
 for k=1:nx-2, A(k+2,k)=2;end
 for k=1:nx-3, A(k+3,k)=-3;end
-deff('xd=f(t,x)','xd=A*x')
+function xd=f(t,x)
+    xd=A*x
+endfunction
 ml=3;mu=2;
 %ODEOPTIONS(11:12)=[ml,mu];
 for i=1:nx;
     for j=1:nx;
-if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
+        if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
 end;end;
 // J is a ml+mu+1 x ny matrix.
 // Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
@@ -198,18 +214,23 @@ end;end;
 // Column 2 of J is made of mu-1 zeros followed by df1/dx2, df2/dx2, ...
 // etc...
 %ODEOPTIONS(6)=1;%ODEOPTIONS(11:12)=[-1,-1];
-deff('jj=j1(t,x)','jj=A')
-xnotband=ode('st',x0,t0,t,f,j1);
+function jj=j1(t,x)
+    jj=A
+endfunction
+xnotband=ode("stiff",x0,t0,t,f,j1);
 %ODEOPTIONS(6)=4;//banded jacobian external given
 %ODEOPTIONS(11:12)=[3,2];
-deff('jj=j(t,x)','jj=J')
-xband=ode('st',x0,t0,t,f,j);
+function jj=j(t,x)
+    jj=J
+endfunction
+xband=ode("stiff",x0,t0,t,f,j);
 if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 %ODEOPTIONS(6)=5;//banded jacobian evaluated
 %ODEOPTIONS(11:12)=[3,2];
-deff('jj=j(t,x)','jj=J')
-xband=ode('st',x0,t0,t,f,j);
-
+function jj=j(t,x)
+    jj=J
+endfunction
+xband=ode("stiff",x0,t0,t,f,j);
 if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 //            Test of %ODEOPTIONS(7)
 //%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
@@ -233,50 +254,60 @@ ww=ode(x0,t0,t,f);norm(wref-ww,1);
 if norm(wref-ode(x0,t0,t,f),1)>Leps then bugmes();quit;end
 //using stiff method
 %ODEOPTIONS(9)=0;
-wref=ode('st',x0,t0,t,f);
+wref=ode("stiff",x0,t0,t,f);
 %ODEOPTIONS(9)=5;
-if norm(wref-ode('st',x0,t0,t,f),1) >Leps then bugmes();quit;end //=0
+if norm(wref-ode("stiff",x0,t0,t,f),1) >Leps then bugmes();quit;end //=0
 %ODEOPTIONS(9)=4;
-if norm(wref-ode('st',x0,t0,t,f),1)  >Leps then bugmes();quit;end  //small
+if norm(wref-ode("stiff",x0,t0,t,f),1)  >Leps then bugmes();quit;end  //small
 //using nonstiff method
 %ODEOPTIONS(8)=0;
-wref=ode('ad',x0,t0,t,f);
+wref=ode("adams",x0,t0,t,f);
 %ODEOPTIONS(8)=12;
-if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end  //=0
+if norm(wref-ode("adams",x0,t0,t,f),1) >Leps then bugmes();quit;end  //=0
 %ODEOPTIONS(8)=5;
-if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end   //small
+if norm(wref-ode("adams",x0,t0,t,f),1) >Leps then bugmes();quit;end   //small
 //mixed
 %ODEOPTIONS(8:9)=[5,12];
 wref=ode(x0,t0,t,f);
 %ODEOPTIONS(8:9)=[4,10];
 if norm(ode(x0,t0,t,f)-wref,1)>Leps then bugmes();quit;end   //small
 A=diag([-10,-0.01,-1]);
-deff('uu=u(t)','uu=sin(t)');
+function uu=u(t)
+    uu=sin(t)
+endfunction
 B=rand(3,1);
-deff('y=f(t,x)','y=A*x+B*u(t)')
+function y=f(t,x)
+    y=A*x+B*u(t)
+endfunction
 %ODEOPTIONS(1)=2;
-yy1=ode('stiff',[1;1;1],0,1,f);
-yy2=ode('stiff',[1;1;1],0,2,f);
+yy1=ode("stiff",[1;1;1],0,1,f);
+yy2=ode("stiff",[1;1;1],0,2,f);
 %ODEOPTIONS(1)=3;
-yy1=ode('stiff',[1;1;1],0,1,f);
-yy2=ode('stiff',[1;1;1],0,2,f);
+yy1=ode("stiff",[1;1;1],0,1,f);
+yy2=ode("stiff",[1;1;1],0,2,f);
 clear %ODEOPTIONS;
-rand('seed',0);rand('normal');
+rand("seed",0);rand("normal");
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
-deff('y=f(t,x)','y=A*x')
-deff('J=j(t,x)','J=A')
+function y=f(t,x)
+    y=A*x
+endfunction
+function J=j(t,x)
+    J=A
+endfunction
 //%ODEOPTIONS(1)=1;
-y2=ode('stiff',ones(nx,1),0,2,f,j);
-[y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
-y2p=ode('stiff',y1,1,2,f,j,w,iw);
-y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
+y2=ode("stiff",ones(nx,1),0,2,f,j);
+[y1,w,iw]=ode("stiff",ones(nx,1),0,1,f,j);
+y2p=ode("stiff",y1,1,2,f,j,w,iw);
+y12=ode("stiff",ones(nx,1),0,[1,2],f,j);
 norm(y12(:,2)-y2p);
-yaf=ode('adams',ones(nx,1),0,2,f,j);
-yaj=ode('adams',ones(nx,1),0,2,f,j);
-ysf=ode('stiff',ones(nx,1),0,2,f,j);
-ysj=ode('stiff',ones(nx,1),0,2,f,j);
-deff('xd=f(t,x)','xd=A*x+B*sin(3*t)')
+yaf=ode("adams",ones(nx,1),0,2,f,j);
+yaj=ode("adams",ones(nx,1),0,2,f,j);
+ysf=ode("stiff",ones(nx,1),0,2,f,j);
+ysj=ode("stiff",ones(nx,1),0,2,f,j);
+function xd=f(t,x)
+    xd=A*x+B*sin(3*t)
+endfunction
 A=rand(10,10)-4.5*eye();B=rand(10,1);
 x=ode(ones(10,1),0,[1,2,3],f);
 //link('fexab.o','fexab')
-if norm(x-ode(ones(10,1),0,[1,2,3],'fexab'),1) > 1.d-9 then bugmes();quit;end
+if norm(x-ode(ones(10,1),0,[1,2,3],"fexab"),1) > 1.d-9 then bugmes();quit;end
index b36c3dc..3900e6c 100644 (file)
@@ -35,56 +35,56 @@ function [jac]=j(t,y)
     jac(2,2)=-jac(1,2)-jac(3,2)
     jac(2,3)=-jac(1,3)
 endfunction
-//    
+//
 y0=[1;0;0];
 t0=0;
 t1=[0.4,4];
-nt=size(t1,'*');
-//    solution 
+nt=size(t1,"*");
+//    solution
 yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 //
 //  1. fortran called by fydot, without jacobian
-y1=ode(y0,t0,t1,'fex');
+y1=ode(y0,t0,t1,"fex");
 if max(y1-yref) > Leps then pause,end
 //  2. fortran called by fydot, type given (stiff), no jacobian
-y2=ode('stiff',y0,t0,t1,'fex');
+y2=ode("stiff",y0,t0,t1,"fex");
 if max(y2-yref) > Leps then pause,end
 //  3. fortran called by fydot , fjac, type given
-y3=ode('stiff',y0,t0,t1,'fex','jex');
+y3=ode("stiff",y0,t0,t1,"fex","jex");
 if max(y3-yref) > Leps then pause,end
 //   hot restart
-[z,w,iw]=ode('stiff',y0,0,0.4,'fex','jex');
-z=ode('stiff',z,0.4,4,'fex','jex',w,iw);
+[z,w,iw]=ode("stiff",y0,0,0.4,"fex","jex");
+z=ode("stiff",z,0.4,4,"fex","jex",w,iw);
 if max(z-y3(:,2)) > %eps then pause,end
 
-[y1,w,iw]=ode(y0,t0,t1(1),'fex');
-y2=ode(y0,t1(1),t1(2:nt),'fex',w,iw);
+[y1,w,iw]=ode(y0,t0,t1(1),"fex");
+y2=ode(y0,t1(1),t1(2:nt),"fex",w,iw);
 if max([y1 y2]-yref) > Leps then pause,end
 
-[y1,w,iw]=ode(y0,t0,t1(1),'fex','jex');
-y2=ode(y0,t1(1),t1(2:nt),'fex','jex',w,iw);
+[y1,w,iw]=ode(y0,t0,t1(1),"fex","jex");
+y2=ode(y0,t1(1),t1(2:nt),"fex","jex",w,iw);
 if max([y1 y2]-yref) > Leps then pause,end
 
 //   variation of tolerances
 atol=[0.001,0.0001,0.001];
 rtol=atol;
-//    externals 
+//    externals
 //  4. type given , scilab lhs ,jacobian not passed
-y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
+y4=ode("stiff",y0,t0,t1(1),atol,rtol,f);
 if max(y4(:,1)-yref(:,1)) > 0.01 then pause,end
 //  5. type non given, rhs and scilab jacobian
 y5=ode(y0,t0,t1,f,j);
 if max(y5-yref) > Leps then pause,end
 //  6. type given (stiff),rhs and jacobian  by scilab
-y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
+y6=ode("stiff",y0,t0,t1,0.00001,0.00001,f,j);
 if (y6-yref) > 2*0.00001 then pause,end
 //  7. matrix rhs, type given(adams),jacobian not passed
-// 
+//
 a=rand(3,3);ea=expm(a);
 function [ydot]=f(t,y)
     ydot=a*y
 endfunction
-t1=1;y=ode('adams',eye(a),t0,t1,f);
+t1=1;y=ode("adams",eye(a),t0,t1,f);
 if max(ea-y) > Leps then pause,end
 //
 //   DAE's
@@ -105,8 +105,8 @@ function [p]=aplusp(t,y,p)
 endfunction
 
 //        %ODEOPTIONS tests
-//       
-rand('seed',0);rand('normal');
+//
+rand("seed",0);rand("normal");
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
 function y=f(t,x)
     y=A*x
@@ -117,7 +117,7 @@ function J=j(t,x)
 endfunction
 
 x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
-nt=size(t,'*');
+nt=size(t,"*");
 eAt=expm(A*t(nt));
 
 //        Test itask=%ODEOPTIONS(1)
@@ -145,7 +145,7 @@ if norm(xlast-expm(A*xft(1,q))*x0)>Leps then pause,end
 //itask=3;   ---> solution at first mesh point beyond t=tmax
 %ODEOPTIONS(1)=3;
 x3=ode(x0,t0,tmax,f);
-if norm(x3(2:nx+1)-xlast,1)>Leps then pause,end 
+if norm(x3(2:nx+1)-xlast,1)>Leps then pause,end
 
 //itask=4; test with %ODEOPTIONS(2)=tcrit
 %ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
@@ -163,7 +163,7 @@ endfunction
 x42=ode(x0,t0,t,fcrit);
 if chk==1 then pause,end
 [p,q]=size(x42);
-if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then pause,end    
+if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then pause,end
 
 //itask=5; test with %ODEOPTIONS(2)=tcrit
 %ODEOPTIONS(1)=5;  //---> computation at mesh points and t>tcrit are not called
@@ -184,29 +184,29 @@ if norm(x35-xf,1)>Leps then pause,end
 //test of %ODEOPTIONS(6)=jacflag
 %ODEOPTIONS(6)=1;//Jacobian given
 %ODEOPTIONS(3:5)=[0 0 0];
-x61=ode('stiff',x0,t0,t,f,j);   //with Jacobian
+x61=ode("stiff",x0,t0,t,f,j);   //with Jacobian
 if norm(x61-xf,1)>10*Leps then pause,end
 
 
 %ODEOPTIONS(6)=0; // jacobian nor called nor estimated
-x60=ode('stiff',x0,t0,t,f,j);   //Jacobian not used (warning)
-x60=ode('stiff',x0,t0,t,f);    //Jacobian not used
+x60=ode("stiff",x0,t0,t,f,j);   //Jacobian not used (warning)
+x60=ode("stiff",x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then pause,end
 
 
 %ODEOPTIONS(6)=1;//Jacobian estimated
-x60=ode('stiff',x0,t0,t,f);
+x60=ode("stiff",x0,t0,t,f);
 if norm(x60-x61,1)>10*Leps then pause,end
 
 //test of %ODEOPTIONS(6)=jacflag   (adams)
 %ODEOPTIONS(6)=1;//with given Jacobian
-x60=ode('adams',x0,t0,t,f,j) ;  
+x60=ode("adams",x0,t0,t,f,j) ;
 if norm(x60-x61,1)>10*Leps then pause,end
 
 
 %ODEOPTIONS(6)=0;// jacobian nor called nor estimated
-x60=ode('adams',x0,t0,t,f,j);   //Jacobian not used (warning)
-x60=ode('adams',x0,t0,t,f);    //Jacobian not used
+x60=ode("adams",x0,t0,t,f,j);   //Jacobian not used (warning)
+x60=ode("adams",x0,t0,t,f);    //Jacobian not used
 if norm(x60-x61,1)>10*Leps then pause,end
 
 // test lsoda
@@ -241,7 +241,7 @@ ml=3;mu=2;
 %ODEOPTIONS(11:12)=[ml,mu];
 for i=1:nx;
     for j=1:nx;
-if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
+        if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
 end;end;
 // J is a ml+mu+1 x ny matrix.
 // Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
@@ -253,7 +253,7 @@ function jj=j1(t,x)
     jj=A
 endfunction
 
-xnotband=ode('stiff',x0,t0,t,f,j1);
+xnotband=ode("stiff",x0,t0,t,f,j1);
 
 
 %ODEOPTIONS(6)=4;//banded jacobian external given
@@ -262,7 +262,7 @@ function jj=j(t,x)
     jj=J
 endfunction
 
-xband=ode('stiff',x0,t0,t,f,j);
+xband=ode("stiff",x0,t0,t,f,j);
 if norm(xnotband-xband,1)>Leps then pause,end
 
 %ODEOPTIONS(6)=5;//banded jacobian evaluated
@@ -270,11 +270,11 @@ if norm(xnotband-xband,1)>Leps then pause,end
 function jj=j(t,x)
     jj=J
 endfunction
-xband=ode('stiff',x0,t0,t,f,j);
+xband=ode("stiff",x0,t0,t,f,j);
 if norm(xnotband-xband,1)>Leps then pause,end
 
 
-//            Test of %ODEOPTIONS(7)  
+//            Test of %ODEOPTIONS(7)
 //%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
 itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 //provisional values as default
@@ -283,7 +283,7 @@ jacflag=2;ml=-1;mu=-1;
 %ODEOPTIONS(7)=10;
 //ode(x0,t0,t,f);  // ---> Non convergence
 
-//            Test of %ODEOPTIONS(8:9)  
+//            Test of %ODEOPTIONS(8:9)
 //%ODEOPTIONS(8:9)=[maxordn,maxords] ---> maximum order for nonstiff and stiff
 itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 //provisional values as default
@@ -301,23 +301,23 @@ if norm(wref-ode(x0,t0,t,f),1)>Leps then pause,end
 //using stiff method
 
 %ODEOPTIONS(9)=0;
-wref=ode('stiff',x0,t0,t,f);
+wref=ode("stiff",x0,t0,t,f);
 %ODEOPTIONS(9)=5;
-if norm(wref-ode('stiff',x0,t0,t,f),1) >Leps then pause,end //=0
+if norm(wref-ode("stiff",x0,t0,t,f),1) >Leps then pause,end //=0
 %ODEOPTIONS(9)=4;
-if norm(wref-ode('stiff',x0,t0,t,f),1)  >Leps then pause,end  //small
+if norm(wref-ode("stiff",x0,t0,t,f),1)  >Leps then pause,end  //small
 
 
 //using nonstiff method
 
 %ODEOPTIONS(8)=0;
-wref=ode('ad',x0,t0,t,f);
+wref=ode("adams",x0,t0,t,f);
 %ODEOPTIONS(8)=12;
-if norm(wref-ode('adams',x0,t0,t,f),1) >Leps then pause,end  //=0
+if norm(wref-ode("adams",x0,t0,t,f),1) >Leps then pause,end  //=0
 %ODEOPTIONS(8)=5;
-if norm(wref-ode('adams',x0,t0,t,f),1) >Leps then pause,end   //small
+if norm(wref-ode("adams",x0,t0,t,f),1) >Leps then pause,end   //small
 
-//mixed 
+//mixed
 %ODEOPTIONS(8:9)=[5,12];
 wref=ode(x0,t0,t,f);
 %ODEOPTIONS(8:9)=[4,10];
@@ -334,14 +334,14 @@ function y=f(t,x)
 endfunction
 
 %ODEOPTIONS(1)=2;
-yy1=ode('stiff',[1;1;1],0,1,f);
-yy2=ode('stiff',[1;1;1],0,2,f);
+yy1=ode("stiff",[1;1;1],0,1,f);
+yy2=ode("stiff",[1;1;1],0,2,f);
 %ODEOPTIONS(1)=3;
-yy1=ode('stiff',[1;1;1],0,1,f);
-yy2=ode('stiff',[1;1;1],0,2,f);
+yy1=ode("stiff",[1;1;1],0,1,f);
+yy2=ode("stiff",[1;1;1],0,2,f);
 
 clear %ODEOPTIONS;
-rand('seed',0);rand('normal');
+rand("seed",0);rand("normal");
 nx=20;A=rand(nx,nx);A=A-4.5*eye();
 function y=f(t,x)
     y=A*x
@@ -350,15 +350,15 @@ function J=j(t,x)
     J=A
 endfunction
 //%ODEOPTIONS(1)=1;
-y2=ode('stiff',ones(nx,1),0,2,f,j);
-[y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
-y2p=ode('stiff',y1,1,2,f,j,w,iw);
-y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
+y2=ode("stiff",ones(nx,1),0,2,f,j);
+[y1,w,iw]=ode("stiff",ones(nx,1),0,1,f,j);
+y2p=ode("stiff",y1,1,2,f,j,w,iw);
+y12=ode("stiff",ones(nx,1),0,[1,2],f,j);
 norm(y12(:,2)-y2p);
-yaf=ode('adams',ones(nx,1),0,2,f,j);
-yaj=ode('adams',ones(nx,1),0,2,f,j);
-ysf=ode('stiff',ones(nx,1),0,2,f,j);
-ysj=ode('stiff',ones(nx,1),0,2,f,j);
+yaf=ode("adams",ones(nx,1),0,2,f,j);
+yaj=ode("adams",ones(nx,1),0,2,f,j);
+ysf=ode("stiff",ones(nx,1),0,2,f,j);
+ysj=ode("stiff",ones(nx,1),0,2,f,j);
 
 function xd=f(t,x)
     xd=A*x+B*sin(3*t)
@@ -367,5 +367,5 @@ endfunction
 A=rand(10,10)-4.5*eye();B=rand(10,1);
 x=ode(ones(10,1),0,[1,2,3],f);
 //link('fexab.o','fexab')
-if norm(x-ode(ones(10,1),0,[1,2,3],'fexab'),1) > 1.d-9 then pause,end
+if norm(x-ode(ones(10,1),0,[1,2,3],"fexab"),1) > 1.d-9 then pause,end
 
index a788a95..d21c442 100644 (file)
@@ -8,17 +8,17 @@
 // =============================================================================
 
 // <-- CLI SHELL MODE -->
-
+warning("off")
 ilib_verbose(0);
 
 version = getversion("scilab");
 
- // to check that ode works
- // ---------- Simple one dimension ODE (Scilab function external)
- // dy/dt=y^2-y sin(t)+cos(t), y(0)=0
- function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
- y0=0;t0=0;t=0:0.1:%pi;
- y=ode(y0,t0,t,f);
+// to check that ode works
+// ---------- Simple one dimension ODE (Scilab function external)
+// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
+function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
+y0=0;t0=0;t=0:0.1:%pi;
+y=ode(y0,t0,t,f);
 assert_checkalmostequal(size(y), [1 32] , %eps, [], "matrix");
 clear y;
 clear t;
@@ -26,25 +26,25 @@ clear t;
 // create functions
 cd TMPDIR;
 
-CC=['void fex1(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = -0.04*y[0] + 1.0e+4*y[1]*y[2];'
-    '   ydot[2] = 3.0e+7*y[1]*y[1];'
-    '   ydot[1] = -ydot[0] - ydot[2];'
-    '}'];
-mputl(CC,TMPDIR+'/fex1.c');
-ilib_for_link('fex1','fex1.c',[],'c');
+CC=["void fex1(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = -0.04*y[0] + 1.0e+4*y[1]*y[2];"
+"   ydot[2] = 3.0e+7*y[1]*y[1];"
+"   ydot[1] = -ydot[0] - ydot[2];"
+"}"];
+mputl(CC,TMPDIR+"/fex1.c");
+ilib_for_link("fex1","fex1.c",[],"c");
 exec loader.sce;
 
-C=[ 'void fex2(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = y[4]*y[0] + y[5]*y[1]*y[2];'
-    '   ydot[2] = y[3]*y[1]*y[1];'
-    '   ydot[1] = -ydot[0] - ydot[2];'
-    '}'];
+C=[ "void fex2(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = y[4]*y[0] + y[5]*y[1]*y[2];"
+"   ydot[2] = y[3]*y[1]*y[1];"
+"   ydot[1] = -ydot[0] - ydot[2];"
+"}"];
 
-mputl(C,TMPDIR+'/fex2.c');
-ilib_for_link('fex2','fex2.c',[],'c');
+mputl(C,TMPDIR+"/fex2.c");
+ilib_for_link("fex2","fex2.c",[],"c");
 exec loader.sce;
 clear f;
 function ydot = f(t,yin)
@@ -104,10 +104,10 @@ resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
 
 // f as a string (dynamic link function)
-res  = ode(y, t,tout, rtol, atol, 'fex1');
+res  = ode(y, t,tout, rtol, atol, "fex1");
 // f as a list(string,...
 if version(1) > 5 then
-    res2 = ode(y, t, tout, rtol, atol, list('fex2', 3.0d+7, -0.04, 1.0d+4));
+    res2 = ode(y, t, tout, rtol, atol, list("fex2", 3.0d+7, -0.04, 1.0d+4));
 end
 // f as a macro
 res3 = ode(y, t, tout, rtol, atol, f);
@@ -116,7 +116,7 @@ res4 = ode(y, t, tout, rtol, atol, list(f1, 3.0d+7, -0.04, 1.0d+4));
 args = [ 3.0d+7, -0.04, 1.0d+4 ];
 res5 = ode(y, t, tout, rtol, atol, list(f2, args));
 // f as a string (static link function)
-res6 = ode(y, t,tout, rtol, atol, 'fex');
+res6 = ode(y, t,tout, rtol, atol, "fex");
 
 // check results
 assert_checkalmostequal(resDoc, res, 2d-7, [], "matrix"); // There are a little diff between resDoc and res
@@ -151,48 +151,48 @@ assert_checkalmostequal(yout2(3*12+1:3*13), yout1, %eps, [], "matrix");
 //for i=1:12, assert_checkequal(poly(res(:,i), "s", "coeff"), res7(i), %eps); end
 
 //*************************** function Jac and lsode ********************************/
-CC=['void jac(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)'
-    '{'
-    '  j[0] = y[6]; '
-    '  j[1] = y[7]; '
-    '  j[2] = y[8]; '
-    '  j[3] = y[9]; '
-    '}'];
-
-mputl(CC,TMPDIR+'/jac.c');
-ilib_for_link('jac','jac.c',[],'c');
+CC=["void jac(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
+"{"
+"  j[0] = y[6]; "
+"  j[1] = y[7]; "
+"  j[2] = y[8]; "
+"  j[3] = y[9]; "
+"}"];
+
+mputl(CC,TMPDIR+"/jac.c");
+ilib_for_link("jac","jac.c",[],"c");
 exec loader.sce;
 
-CC=['void jac2(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)'
-    '{'
-    '  j[0] = 10; '
-    '  j[1] = 0; '
-    '  j[2] = 0; '
-    '  j[3] = -1; '
-    '}'];
-mputl(CC,TMPDIR+'/jac2.c');
-ilib_for_link('jac2','jac2.c',[],'c');
+CC=["void jac2(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
+"{"
+"  j[0] = 10; "
+"  j[1] = 0; "
+"  j[2] = 0; "
+"  j[3] = -1; "
+"}"];
+mputl(CC,TMPDIR+"/jac2.c");
+ilib_for_link("jac2","jac2.c",[],"c");
 exec loader.sce;
 
 // ydot = A * y
-C=[ 'void fext(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = y[2]*y[0] + y[4]*y[1];'
-    '   ydot[1] = y[3]*y[0] + y[5]*y[1];'
-    '}'];
-
-mputl(C,TMPDIR+'/fext.c');
-ilib_for_link('fext','fext.c',[],'c');
+C=[ "void fext(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = y[2]*y[0] + y[4]*y[1];"
+"   ydot[1] = y[3]*y[0] + y[5]*y[1];"
+"}"];
+
+mputl(C,TMPDIR+"/fext.c");
+ilib_for_link("fext","fext.c",[],"c");
 exec loader.sce;
 
-C=[ 'void fext2(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = 10*y[0] + 0*y[1];'
-    '   ydot[1] = 0*y[0] + (-1)*y[1];'
-    '}'];
+C=[ "void fext2(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = 10*y[0] + 0*y[1];"
+"   ydot[1] = 0*y[0] + (-1)*y[1];"
+"}"];
 
-mputl(C,TMPDIR+'/fext2.c');
-ilib_for_link('fext2','fext2.c',[],'c');
+mputl(C,TMPDIR+"/fext2.c");
+ilib_for_link("fext2","fext2.c",[],"c");
 exec loader.sce;
 clear f;
 
@@ -211,11 +211,11 @@ t=1;
 
 res  = ode("stiff", y0, t0, t, f, Jacobian);
 if version(1) > 5 then
-    res1 = ode("stiff", y0, t0, t, list('fext', 10, 0, 0,-1), list('jac', A));
+    res1 = ode("stiff", y0, t0, t, list("fext", 10, 0, 0,-1), list("jac", A));
 end
-res2 = ode("stiff", y0, t0, t, 'fext2', 'jac2');
-res3 = ode("stiff", y0, t0, t, f, 'jac2');
-res4 = ode("stiff", y0, t0, t, 'fext2', Jacobian);
+res2 = ode("stiff", y0, t0, t, "fext2", "jac2");
+res3 = ode("stiff", y0, t0, t, f, "jac2");
+res4 = ode("stiff", y0, t0, t, "fext2", Jacobian);
 
 assert_checkalmostequal(res, expm(A*t)*y0, 1.0D-7, [], "matrix");
 if version(1) > 5 then
@@ -252,15 +252,15 @@ assert_checkalmostequal(y, y1, %eps, [], "matrix");
 y0=1;
 ng=1;
 
-C=[ 'void fextern(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    'int i = 0;'
-    'for( i = 0; i < *neq; i++)'
-    '   ydot[i] = y[i];'
-    '}'];
+C=[ "void fextern(int* neq, double* t, double* y, double* ydot)"
+"{"
+"int i = 0;"
+"for( i = 0; i < *neq; i++)"
+"   ydot[i] = y[i];"
+"}"];
 
-mputl(C,TMPDIR+'/fextern.c');
-ilib_for_link('fextern','fextern.c',[],'c');
+mputl(C,TMPDIR+"/fextern.c");
+ilib_for_link("fextern","fextern.c",[],"c");
 exec loader.sce;
 clear f;
 function ydot=f(t,y)
@@ -270,14 +270,14 @@ endfunction
 clear g;
 // check rd result
 function z=g(t,y)
-z=y-2;
+    z=y-2;
 endfunction
 [y,rd]=ode("root",y0,0,2,f,ng,g);
 assert_checkequal( rd(2) <> 1 , %f);
 
 clear g;
 function z=g(t,y)
-z=y-[2;2;33];
+    z=y-[2;2;33];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,3,g);
 assert_checkequal( rd(2) <> 1 , %f);
@@ -285,7 +285,7 @@ assert_checkequal( rd(3) <> 2 , %f);
 
 clear g;
 function z=g(t,y)
-z=y-[2;2;2];
+    z=y-[2;2;2];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,3,g);
 assert_checkequal( rd(2) <> 1 , %f);
@@ -294,7 +294,7 @@ assert_checkequal( rd(4) <> 3 , %f);
 
 clear g;
 function z=g(t,y)
-z=y-[2;6;2;2];
+    z=y-[2;6;2;2];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,4,g);
 assert_checkequal( rd(2) <> 1 , %f);
@@ -336,14 +336,14 @@ resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
 //   at t =  4.0000e+10
 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
 
-G=[ 'void gex(int* neq, double* t, double* y, int* ng, double* gout)'
-    '{'
-        'gout[0] = y[0] - 1.0e-4;'
-        'gout[1] = y[2] - 1.0e-2;'
-    '}'];
+G=[ "void gex(int* neq, double* t, double* y, int* ng, double* gout)"
+"{"
+"gout[0] = y[0] - 1.0e-4;"
+"gout[1] = y[2] - 1.0e-2;"
+"}"];
 
-mputl(G,TMPDIR+'/gex.c');
-ilib_for_link('gex','gex.c',[],'c');
+mputl(G,TMPDIR+"/gex.c");
+ilib_for_link("gex","gex.c",[],"c");
 exec loader.sce;
 
 y(1) = 1;
@@ -351,9 +351,9 @@ y(2) = 0;
 y(3) = 0;
 t0   = 0;
 
-[yout1,rd1,w,iw] = ode("root", y, t0, tout, 'fex1', 2, 'gex');
+[yout1,rd1,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex");
 assert_checkalmostequal(rd1(1), 2.64d-01, 1d-4);
-[yout2,rd2,w,iw] = ode("root", y, t0, tout, 'fex1', 2, 'gex', w, iw);
+[yout2,rd2,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex", w, iw);
 assert_checkalmostequal(rd2(1), 2.0795776d+07, 4d-5);
 err = execstr("[yout3,rd,w,iw] = ode(""root"", y, t0, tout, ""fex1"", 2, ""gex"", w, iw);","errcatch");
 assert_checkequal( err == 0 , %f);
@@ -385,3 +385,4 @@ rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094
 assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");
+warning("on")
index e19a95a..63c2529 100644 (file)
@@ -7,6 +7,7 @@
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 // <-- CLI SHELL MODE -->
+warning("off")
 ilib_verbose(0);
 version = getversion("scilab");
  // to check that ode works
@@ -309,6 +310,9 @@ assert_checkalmostequal(rd1(1), 2.64d-01, 1d-4);
 [yout2,rd2,w,iw] = ode("root", y, t0, tout, 'fex1', 2, 'gex', w, iw);
 assert_checkalmostequal(rd2(1), 2.0795776d+07, 4d-5);
 err = execstr("[yout3,rd,w,iw] = ode(""root"", y, t0, tout, ""fex1"", 2, ""gex"", w, iw);","errcatch");
+      t n est pas entre tcur - hu (= r1) et tcur (=r2)
+      where r1 is :   0.2033607066275D+08   and r2 :   0.2109164231072D+08      
+Illegal input detected (see printed message).
 assert_checkequal( err == 0 , %f);
 // check results
 assert_checkalmostequal(resDocRoot(:,1), yout1, 2.0D-8, [], "matrix");
@@ -331,3 +335,4 @@ rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094
 assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");
+warning("on")
index e19a95a..c12a0ff 100644 (file)
@@ -7,37 +7,38 @@
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 // <-- CLI SHELL MODE -->
+warning("off")
 ilib_verbose(0);
 version = getversion("scilab");
- // to check that ode works
- // ---------- Simple one dimension ODE (Scilab function external)
- // dy/dt=y^2-y sin(t)+cos(t), y(0)=0
- function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
- y0=0;t0=0;t=0:0.1:%pi;
- y=ode(y0,t0,t,f);
+// to check that ode works
+// ---------- Simple one dimension ODE (Scilab function external)
+// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
+function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
+y0=0;t0=0;t=0:0.1:%pi;
+y=ode(y0,t0,t,f);
 assert_checkalmostequal(size(y), [1 32] , %eps, [], "matrix");
 clear y;
 clear t;
 //*************************** function F and lsoda ********************************/
 // create functions
 cd TMPDIR;
-CC=['void fex1(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = -0.04*y[0] + 1.0e+4*y[1]*y[2];'
-    '   ydot[2] = 3.0e+7*y[1]*y[1];'
-    '   ydot[1] = -ydot[0] - ydot[2];'
-    '}'];
-mputl(CC,TMPDIR+'/fex1.c');
-ilib_for_link('fex1','fex1.c',[],'c');
+CC=["void fex1(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = -0.04*y[0] + 1.0e+4*y[1]*y[2];"
+"   ydot[2] = 3.0e+7*y[1]*y[1];"
+"   ydot[1] = -ydot[0] - ydot[2];"
+"}"];
+mputl(CC,TMPDIR+"/fex1.c");
+ilib_for_link("fex1","fex1.c",[],"c");
 exec loader.sce;
-C=[ 'void fex2(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = y[4]*y[0] + y[5]*y[1]*y[2];'
-    '   ydot[2] = y[3]*y[1]*y[1];'
-    '   ydot[1] = -ydot[0] - ydot[2];'
-    '}'];
-mputl(C,TMPDIR+'/fex2.c');
-ilib_for_link('fex2','fex2.c',[],'c');
+C=[ "void fex2(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = y[4]*y[0] + y[5]*y[1]*y[2];"
+"   ydot[2] = y[3]*y[1]*y[1];"
+"   ydot[1] = -ydot[0] - ydot[2];"
+"}"];
+mputl(C,TMPDIR+"/fex2.c");
+ilib_for_link("fex2","fex2.c",[],"c");
 exec loader.sce;
 clear f;
 function ydot = f(t,yin)
@@ -92,10 +93,10 @@ resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
 //   at t =  4.0000e+10
 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
 // f as a string (dynamic link function)
-res  = ode(y, t,tout, rtol, atol, 'fex1');
+res  = ode(y, t,tout, rtol, atol, "fex1");
 // f as a list(string,...
 if version(1) > 5 then
-    res2 = ode(y, t, tout, rtol, atol, list('fex2', 3.0d+7, -0.04, 1.0d+4));
+    res2 = ode(y, t, tout, rtol, atol, list("fex2", 3.0d+7, -0.04, 1.0d+4));
 end
 // f as a macro
 res3 = ode(y, t, tout, rtol, atol, f);
@@ -104,7 +105,7 @@ res4 = ode(y, t, tout, rtol, atol, list(f1, 3.0d+7, -0.04, 1.0d+4));
 args = [ 3.0d+7, -0.04, 1.0d+4 ];
 res5 = ode(y, t, tout, rtol, atol, list(f2, args));
 // f as a string (static link function)
-res6 = ode(y, t,tout, rtol, atol, 'fex');
+res6 = ode(y, t,tout, rtol, atol, "fex");
 // check results
 assert_checkalmostequal(resDoc, res, 2d-7, [], "matrix"); // There are a little diff between resDoc and res
 if version(1) > 5 then
@@ -132,42 +133,42 @@ assert_checkalmostequal(yout2(3*12+1:3*13), yout1, %eps, [], "matrix");
 //res7 = ode(yy, t, tout, rtol, atol, list('fex2', 3.0d+7, -0.04, 1.0d+4));
 //for i=1:12, assert_checkequal(poly(res(:,i), "s", "coeff"), res7(i), %eps); end
 //*************************** function Jac and lsode ********************************/
-CC=['void jac(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)'
-    '{'
-    '  j[0] = y[6]; '
-    '  j[1] = y[7]; '
-    '  j[2] = y[8]; '
-    '  j[3] = y[9]; '
-    '}'];
-mputl(CC,TMPDIR+'/jac.c');
-ilib_for_link('jac','jac.c',[],'c');
+CC=["void jac(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
+"{"
+"  j[0] = y[6]; "
+"  j[1] = y[7]; "
+"  j[2] = y[8]; "
+"  j[3] = y[9]; "
+"}"];
+mputl(CC,TMPDIR+"/jac.c");
+ilib_for_link("jac","jac.c",[],"c");
 exec loader.sce;
-CC=['void jac2(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)'
-    '{'
-    '  j[0] = 10; '
-    '  j[1] = 0; '
-    '  j[2] = 0; '
-    '  j[3] = -1; '
-    '}'];
-mputl(CC,TMPDIR+'/jac2.c');
-ilib_for_link('jac2','jac2.c',[],'c');
+CC=["void jac2(int* neq, double* t, double* y, int* mu, int* ml, double* j, int nj)"
+"{"
+"  j[0] = 10; "
+"  j[1] = 0; "
+"  j[2] = 0; "
+"  j[3] = -1; "
+"}"];
+mputl(CC,TMPDIR+"/jac2.c");
+ilib_for_link("jac2","jac2.c",[],"c");
 exec loader.sce;
 // ydot = A * y
-C=[ 'void fext(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = y[2]*y[0] + y[4]*y[1];'
-    '   ydot[1] = y[3]*y[0] + y[5]*y[1];'
-    '}'];
-mputl(C,TMPDIR+'/fext.c');
-ilib_for_link('fext','fext.c',[],'c');
+C=[ "void fext(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = y[2]*y[0] + y[4]*y[1];"
+"   ydot[1] = y[3]*y[0] + y[5]*y[1];"
+"}"];
+mputl(C,TMPDIR+"/fext.c");
+ilib_for_link("fext","fext.c",[],"c");
 exec loader.sce;
-C=[ 'void fext2(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    '   ydot[0] = 10*y[0] + 0*y[1];'
-    '   ydot[1] = 0*y[0] + (-1)*y[1];'
-    '}'];
-mputl(C,TMPDIR+'/fext2.c');
-ilib_for_link('fext2','fext2.c',[],'c');
+C=[ "void fext2(int* neq, double* t, double* y, double* ydot)"
+"{"
+"   ydot[0] = 10*y[0] + 0*y[1];"
+"   ydot[1] = 0*y[0] + (-1)*y[1];"
+"}"];
+mputl(C,TMPDIR+"/fext2.c");
+ilib_for_link("fext2","fext2.c",[],"c");
 exec loader.sce;
 clear f;
 function ydot=f(t, y)
@@ -182,11 +183,11 @@ t0=0;
 t=1;
 res  = ode("stiff", y0, t0, t, f, Jacobian);
 if version(1) > 5 then
-    res1 = ode("stiff", y0, t0, t, list('fext', 10, 0, 0,-1), list('jac', A));
+    res1 = ode("stiff", y0, t0, t, list("fext", 10, 0, 0,-1), list("jac", A));
 end
-res2 = ode("stiff", y0, t0, t, 'fext2', 'jac2');
-res3 = ode("stiff", y0, t0, t, f, 'jac2');
-res4 = ode("stiff", y0, t0, t, 'fext2', Jacobian);
+res2 = ode("stiff", y0, t0, t, "fext2", "jac2");
+res3 = ode("stiff", y0, t0, t, f, "jac2");
+res4 = ode("stiff", y0, t0, t, "fext2", Jacobian);
 assert_checkalmostequal(res, expm(A*t)*y0, 1.0D-7, [], "matrix");
 if version(1) > 5 then
     assert_checkalmostequal(res, res1, %eps, [], "matrix");
@@ -216,14 +217,14 @@ assert_checkalmostequal(y, y1, %eps, [], "matrix");
 //*************************** root ********************************/
 y0=1;
 ng=1;
-C=[ 'void fextern(int* neq, double* t, double* y, double* ydot)'
-    '{'
-    'int i = 0;'
-    'for( i = 0; i < *neq; i++)'
-    '   ydot[i] = y[i];'
-    '}'];
-mputl(C,TMPDIR+'/fextern.c');
-ilib_for_link('fextern','fextern.c',[],'c');
+C=[ "void fextern(int* neq, double* t, double* y, double* ydot)"
+"{"
+"int i = 0;"
+"for( i = 0; i < *neq; i++)"
+"   ydot[i] = y[i];"
+"}"];
+mputl(C,TMPDIR+"/fextern.c");
+ilib_for_link("fextern","fextern.c",[],"c");
 exec loader.sce;
 clear f;
 function ydot=f(t,y)
@@ -232,20 +233,20 @@ endfunction
 clear g;
 // check rd result
 function z=g(t,y)
-z=y-2;
+    z=y-2;
 endfunction
 [y,rd]=ode("root",y0,0,2,f,ng,g);
 assert_checkequal( rd(2) <> 1 , %f);
 clear g;
 function z=g(t,y)
-z=y-[2;2;33];
+    z=y-[2;2;33];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,3,g);
 assert_checkequal( rd(2) <> 1 , %f);
 assert_checkequal( rd(3) <> 2 , %f);
 clear g;
 function z=g(t,y)
-z=y-[2;2;2];
+    z=y-[2;2;2];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,3,g);
 assert_checkequal( rd(2) <> 1 , %f);
@@ -253,7 +254,7 @@ assert_checkequal( rd(3) <> 2 , %f);
 assert_checkequal( rd(4) <> 3 , %f);
 clear g;
 function z=g(t,y)
-z=y-[2;6;2;2];
+    z=y-[2;6;2;2];
 endfunction
 [y,rd]=ode("root",y0,0,2,f,4,g);
 assert_checkequal( rd(2) <> 1 , %f);
@@ -292,23 +293,26 @@ resDoc(:,10) = [ 5.283401e-06 ; 2.113371e-11 ; 9.999947e-01 ];
 resDoc(:,11) = [ 4.659031e-07 ; 1.863613e-12 ; 9.999995e-01 ];
 //   at t =  4.0000e+10
 resDoc(:,12) = [ 1.404280e-08 ; 5.617126e-14 ; 1.000000e+00 ];
-G=[ 'void gex(int* neq, double* t, double* y, int* ng, double* gout)'
-    '{'
-        'gout[0] = y[0] - 1.0e-4;'
-        'gout[1] = y[2] - 1.0e-2;'
-    '}'];
-mputl(G,TMPDIR+'/gex.c');
-ilib_for_link('gex','gex.c',[],'c');
+G=[ "void gex(int* neq, double* t, double* y, int* ng, double* gout)"
+"{"
+"gout[0] = y[0] - 1.0e-4;"
+"gout[1] = y[2] - 1.0e-2;"
+"}"];
+mputl(G,TMPDIR+"/gex.c");
+ilib_for_link("gex","gex.c",[],"c");
 exec loader.sce;
 y(1) = 1;
 y(2) = 0;
 y(3) = 0;
 t0   = 0;
-[yout1,rd1,w,iw] = ode("root", y, t0, tout, 'fex1', 2, 'gex');
+[yout1,rd1,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex");
 assert_checkalmostequal(rd1(1), 2.64d-01, 1d-4);
-[yout2,rd2,w,iw] = ode("root", y, t0, tout, 'fex1', 2, 'gex', w, iw);
+[yout2,rd2,w,iw] = ode("root", y, t0, tout, "fex1", 2, "gex", w, iw);
 assert_checkalmostequal(rd2(1), 2.0795776d+07, 4d-5);
 err = execstr("[yout3,rd,w,iw] = ode(""root"", y, t0, tout, ""fex1"", 2, ""gex"", w, iw);","errcatch");
+      t n est pas entre tcur - hu (= r1) et tcur (=r2)
+      where r1 is :   0.2033607066275D+08   and r2 :   0.2109164231072D+08      
+Illegal input detected (see printed message).
 assert_checkequal( err == 0 , %f);
 // check results
 assert_checkalmostequal(resDocRoot(:,1), yout1, 2.0D-8, [], "matrix");
@@ -331,3 +335,4 @@ rkfRes(18:32) = [    0.99166478935902302    0.97384760697150774    0.94630006094
 assert_checkalmostequal(rkRes, rk, %eps * 20, [], "matrix");
 assert_checkalmostequal(rkfRes, rkf, %eps, [], "matrix");
 assert_checkalmostequal(rkf, fixx, %eps, [], "matrix");
+warning("on")
index 0902fee..42d4991 100644 (file)
-
 // =============================================================================
-
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-
 // Copyright (C) 2007-2008 - INRIA
-
 //
-
 //  This file is distributed under the same license as the Scilab package.
-
 // =============================================================================
-
-
 // <-- CLI SHELL MODE -->
-
-
 // Copyright INRIA
-
 function xcd=f(t,xc,xd,iflag)
-
     if iflag==0 then
-
         xcd=fc(t,xc,e(t)-hd(t,xd));
-
     else
-
         xcd=fd(xd,hc(t,xc));
-
     end
-
 endfunction
-
-
 A=[-10,2,3;4,-10,6;7,8,-10];
-
 B=[1;1;1];
-
 C=[1,1,1];
-
 Ad=[1/2,1;0,1/20];
-
 Bd=[1;1];
-
 Cd=[1,1];
-
-
 function st=e(t)
-
     st=sin(3*t)
-
 endfunction
-
-
 function xdot=fc(t,x,u)
-
     xdot=A*x+B*u
-
 endfunction
-
-
 function y=hc(t,x)
-
     y=C*x
-
 endfunction
-
-
 function xp=fd(x,y)
-
     xp=Ad*x + Bd*y
-
 endfunction
-
-
 function u=hd(t,x)
-
     u=Cd*x
-
 endfunction
-
-
 h=0.1;
-
 t0=0;
-
 t=0:0.1:2;
-
 x0c=[0;0;0];
-
 x0d=[0;0];
-
 nd=2;
-
 xcd=odedc([x0c;x0d],nd,h,t0,t,f);
-
 if norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd'),1) > 1.0d-5 then bugmes();quit;end
-
 //(see fydot2.f)
-
 if norm(xcd-odedc([x0c;x0d],nd,h,t0,t,'fcd1'),1) >1.d-5 then bugmes();quit;end
-
-
 function xcd=phis(t,xc,xd,iflag)
-
     if iflag==0 then
-
         xcd=A*xc+B*xd;
-
     else
-
         xcd=1-xd;
-
     end
-
 endfunction
-
-
 t=0:0.1:30;
-
 xcd=odedc([x0c;0],1,1,t0,t,phis);
-
 //link('phis.o','phis')
-
 xcd2=odedc([x0c;0],1,1,t0,t,'phis');
-
-
 if norm(xcd-xcd2,1) > 1.d-5 then bugmes();quit;end
-
 function xd=ff(t,x)
-
     xd=A*x+B*u
-
 endfunction
-
-
 u=1/2;
-
 xn=ode(x0c,t0,t,ff);
-
 //plot2d([t',t',t',t'],[(xcd(1,:))',(xcd(2,:))',(xcd(3,:))',(xcd(4,:))'])
-
-
 function xcd=phit(t,xc,xd,iflag)
-
     if iflag==0 then
-
         xcd=[A*xc(1:3,:)+B*xc(4);xd];
-
     else
-
         xcd=-xd;
-
     end
-
 endfunction
-
-
 xcdt=odedc([x0c;1;1],1,1,t0,t,phit);
-
 //link('phit.o','phit')
-
 xcdt2=odedc([x0c;1;1],1,1,t0,t,'phit');
-
 if norm(xcdt-xcdt2,1) > 1.d-10 then bugmes();quit;end
-
 //plot2d([t',t',t',t'],[(xcdt(1,:))',(xcdt(2,:))',(xcdt(3,:))',(xcdt(4,:))'])
-
-
 xcdt3=odedc('adams',[x0c;1;1],1,1,t0,t,'phit');
-
 if norm(xcdt3-xcdt2,1) > 0.0001 then bugmes();quit;end
-
-
 xcdt4=odedc('fix',[x0c;1;1],1,1,t0,t,'phit');
-
 if norm(xcdt4-xcdt2,1) > 0.001 then bugmes();quit;end
-
-
 xcdt5=odedc('stiff',[x0c;1;1],1,1,t0,t,'phit');
-
 if norm(xcdt5-xcdt2,1) > 0.0001 then bugmes();quit;end
-