Xcos solvers : ODEPACK's 'lsodar' implementation 43/9743/53
Clément DAVID [Tue, 2 Oct 2012 15:25:33 +0000 (17:25 +0200)]
Xcos can now use 'lsodar' solver from Scilab's 'ODEPACK'. Added tests & doc

Change-Id: Ic527681cf6be07baf920653f3a0bb8f9ae75e321

20 files changed:
scilab/CHANGES_5.4.X
scilab/modules/scicos/Makefile.am
scilab/modules/scicos/Makefile.in
scilab/modules/scicos/src/c/differential_equations_f_Import.def [new file with mode: 0644]
scilab/modules/scicos/src/c/lsodar.c [new file with mode: 0644]
scilab/modules/scicos/src/c/lsodar.h [new file with mode: 0644]
scilab/modules/scicos/src/c/scicos.c
scilab/modules/scicos/src/c/scicos.vcxproj
scilab/modules/scicos/src/c/scicos.vcxproj.filters
scilab/modules/scicos/src/scicos_sundials/src/cvode/cvode_impl.h
scilab/modules/scicos/tests/unit_tests/LSodar.dia.ref [new file with mode: 0644]
scilab/modules/scicos/tests/unit_tests/LSodar.tst [new file with mode: 0644]
scilab/modules/xcos/examples/solvers/integLSodar.sce [new file with mode: 0644]
scilab/modules/xcos/examples/solvers/integRK.sce
scilab/modules/xcos/examples/solvers/integSundials.sce [moved from scilab/modules/xcos/examples/solvers/integ.sce with 95% similarity]
scilab/modules/xcos/help/en_US/solvers/CVode.xml
scilab/modules/xcos/help/en_US/solvers/IDA.xml
scilab/modules/xcos/help/en_US/solvers/LSodar.xml [new file with mode: 0644]
scilab/modules/xcos/help/en_US/solvers/Runge-Kutta.xml
scilab/modules/xcos/tests/unit_tests/LSodar_test.zcos [new file with mode: 0644]

index 56a8390..cc875ff 100644 (file)
@@ -53,6 +53,8 @@ Xcos
 
 * Update the Palette icons and the rendering of some blocks.
 
+* Add a new ODE solver : LSodar. Automatically switches methods to efficiently solve both stiff and nonstiff problems. Includes a rootfinding feature.
+
 * Bug #10040 fixed - VARIABLE_DELAY documentation does not fully describe the
                      initial value behavioral.
 
index 942d9cb..49326b5 100644 (file)
@@ -42,7 +42,8 @@ src/c/createblklist.c \
 src/c/var2sci.c \
 src/c/copyvarfromlistentry.c \
 src/c/MlistGetFieldNumber.c \
-src/c/extractblklist.c
+src/c/extractblklist.c \
+src/c/lsodar.c
 
 SCICOS_FORTRAN_SOURCES = src/fortran/scitovv.f \
 src/fortran/vvtosci.f \
index 7009fcf..173ebd3 100644 (file)
@@ -169,10 +169,10 @@ am__libsciscicos_algo_la_SOURCES_DIST = src/c/tree.c \
        src/c/il_state.c src/c/il_sim.c src/c/createblklist.c \
        src/c/var2sci.c src/c/copyvarfromlistentry.c \
        src/c/MlistGetFieldNumber.c src/c/extractblklist.c \
-       src/fortran/scitovv.f src/fortran/vvtosci.f \
-       src/fortran/skipvars.f src/fortran/coselm.f \
-       src/fortran/scitod.f src/fortran/sctree.f src/fortran/scierr.f \
-       src/fortran/ftree2.f src/fortran/ftree3.f \
+       src/c/lsodar.c src/fortran/scitovv.f \
+       src/fortran/vvtosci.f src/fortran/skipvars.f \
+       src/fortran/coselm.f src/fortran/scitod.f src/fortran/sctree.f \
+       src/fortran/scierr.f src/fortran/ftree2.f src/fortran/ftree3.f \
        src/fortran/list2vars.f src/fortran/ftree4.f \
        src/fortran/scitoi.f src/fortran/scifunc.f
 @XCOS_TRUE@am__objects_1 = libsciscicos_algo_la-tree.lo \
@@ -191,7 +191,8 @@ am__libsciscicos_algo_la_SOURCES_DIST = src/c/tree.c \
 @XCOS_TRUE@    libsciscicos_algo_la-var2sci.lo \
 @XCOS_TRUE@    libsciscicos_algo_la-copyvarfromlistentry.lo \
 @XCOS_TRUE@    libsciscicos_algo_la-MlistGetFieldNumber.lo \
-@XCOS_TRUE@    libsciscicos_algo_la-extractblklist.lo
+@XCOS_TRUE@    libsciscicos_algo_la-extractblklist.lo \
+@XCOS_TRUE@    libsciscicos_algo_la-lsodar.lo
 @XCOS_TRUE@am__objects_2 = scitovv.lo vvtosci.lo skipvars.lo coselm.lo \
 @XCOS_TRUE@    scitod.lo sctree.lo scierr.lo ftree2.lo ftree3.lo \
 @XCOS_TRUE@    list2vars.lo ftree4.lo scitoi.lo scifunc.lo
@@ -872,7 +873,8 @@ HELP_CHAPTERLANG = en_US fr_FR pt_BR
 @XCOS_TRUE@src/c/var2sci.c \
 @XCOS_TRUE@src/c/copyvarfromlistentry.c \
 @XCOS_TRUE@src/c/MlistGetFieldNumber.c \
-@XCOS_TRUE@src/c/extractblklist.c
+@XCOS_TRUE@src/c/extractblklist.c \
+@XCOS_TRUE@src/c/lsodar.c
 
 @XCOS_TRUE@SCICOS_FORTRAN_SOURCES = src/fortran/scitovv.f \
 @XCOS_TRUE@src/fortran/vvtosci.f \
@@ -1210,6 +1212,7 @@ distclean-compile:
        -rm -f *.tab.c
 
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libsciscicos_algo_la-MlistGetFieldNumber.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libsciscicos_algo_la-lsodar.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libsciscicos_algo_la-copyvarfromlistentry.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libsciscicos_algo_la-createblklist.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libsciscicos_algo_la-extractblklist.Plo@am__quote@
@@ -1417,6 +1420,13 @@ libsciscicos_algo_la-extractblklist.lo: src/c/extractblklist.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(LIBTOOL)  --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libsciscicos_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libsciscicos_algo_la-extractblklist.lo `test -f 'src/c/extractblklist.c' || echo '$(srcdir)/'`src/c/extractblklist.c
 
+libsciscicos_algo_la-lsodar.lo: src/c/lsodar.c
+@am__fastdepCC_TRUE@   $(LIBTOOL)  --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libsciscicos_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libsciscicos_algo_la-lsodar.lo -MD -MP -MF $(DEPDIR)/libsciscicos_algo_la-lsodar.Tpo -c -o libsciscicos_algo_la-lsodar.lo `test -f 'src/c/lsodar.c' || echo '$(srcdir)/'`src/c/lsodar.c
+@am__fastdepCC_TRUE@   $(am__mv) $(DEPDIR)/libsciscicos_algo_la-lsodar.Tpo $(DEPDIR)/libsciscicos_algo_la-lsodar.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/lsodar.c' object='libsciscicos_algo_la-lsodar.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL)  --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libsciscicos_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libsciscicos_algo_la-lsodar.lo `test -f 'src/c/lsodar.c' || echo '$(srcdir)/'`src/c/lsodar.c
+
 libsciscicos_la-noscicos.lo: src/c/noscicos/noscicos.c
 @am__fastdepCC_TRUE@   $(LIBTOOL)  --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libsciscicos_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libsciscicos_la-noscicos.lo -MD -MP -MF $(DEPDIR)/libsciscicos_la-noscicos.Tpo -c -o libsciscicos_la-noscicos.lo `test -f 'src/c/noscicos/noscicos.c' || echo '$(srcdir)/'`src/c/noscicos/noscicos.c
 @am__fastdepCC_TRUE@   $(am__mv) $(DEPDIR)/libsciscicos_la-noscicos.Tpo $(DEPDIR)/libsciscicos_la-noscicos.Plo
diff --git a/scilab/modules/scicos/src/c/differential_equations_f_Import.def b/scilab/modules/scicos/src/c/differential_equations_f_Import.def
new file mode 100644 (file)
index 0000000..93e16db
--- /dev/null
@@ -0,0 +1,3 @@
+LIBRARY    differential_equations_f.dll
+EXPORTS
+       lsodar_
diff --git a/scilab/modules/scicos/src/c/lsodar.c b/scilab/modules/scicos/src/c/lsodar.c
new file mode 100644 (file)
index 0000000..9ab0603
--- /dev/null
@@ -0,0 +1,564 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2012 - 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
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <machine.h>
+
+#include "lsodar.h"
+
+/* =============================
+ *
+ *            lsodar
+ *
+ * =============================
+ *
+ * Actual solving function, from 'ODEPACK' in 'differential_equations' module
+ */
+
+extern void C2F(lsodar) (LSRhsFn f, int *neq, realtype *y, realtype *t, realtype *tout, int *itol, realtype *reltol, realtype *abstol, enum iTask_t *itask, int *istate, int *iopt, struct rWork_t *rwork, int *lrw, int *iwork, int *liw,  int *jacobian, int *jacType, LSRootFn grblk, int *ng, int *jroot);
+
+/* =============================
+ *
+ *         LSodarCreate
+ *
+ * =============================
+ *
+ * LSodarCreate creates an internal memory block for a problem to be solved by LSODAR.
+ * If successful, LSodarCreate returns a pointer to the problem memory.
+ * This pointer should be passed to LSodarMalloc.
+ * If an initialization error occurs,
+ * LSodarCreate prints an error message to standard err and returns NULL.
+ */
+
+void * LSodarCreate (int * neq, int ng)
+{
+    int lRn, lRs, lIw, lRw;
+
+    /* Allocate the problem memory space */
+    LSodarMem lsodar_mem = NULL;
+    lsodar_mem = (LSodarMem) malloc(sizeof(struct LSodarMemRec));
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, 0, "LSODAR", "LSodarCreate", MSGCV_CVMEM_FAIL);
+        return (NULL);
+    }
+
+    /* Set the 'rwork' and 'iwork' workspaces lengths */
+    lRn = 20 + 16 * (*neq) + 3 * ng;
+    lRs = 22 + ((*neq) * (9 + (*neq))) + 3 * ng;
+    lIw = 20 + (*neq);
+    lRw = max(lRn, lRs);
+
+    /* Copy the variables into the problem memory space */
+    lsodar_mem->nEquations = neq;
+    lsodar_mem->iState = 1;
+    lsodar_mem->iOpt = 0;
+    lsodar_mem->rwork = NULL;
+    lsodar_mem->lrw = lRw;
+    lsodar_mem->iwork = NULL;
+    lsodar_mem->liw = lIw;
+    lsodar_mem->jacobian = 0;
+    lsodar_mem->jacType = 2;
+    lsodar_mem->g_fun = NULL;
+    lsodar_mem->ng_fun = ng;
+    lsodar_mem->jroot = NULL;
+
+    return ((void *) lsodar_mem);
+}
+
+/* Shortcuts to problem memory space parameters */
+# define func    ls_mem->func
+# define nEq     ls_mem->nEquations
+# define yVec    ls_mem->yVector
+# define tStart  ls_mem->tStart
+# define tEnd    ls_mem->tEnd
+# define iTol    ls_mem->iTol
+# define relTol  ls_mem->relTol
+# define absTol  ls_mem->absTol
+# define iState  ls_mem->iState
+# define iOpt    ls_mem->iOpt
+# define rwork   ls_mem->rwork
+# define lrw     ls_mem->lrw
+# define iwork   ls_mem->iwork
+# define liw     ls_mem->liw
+# define jac     ls_mem->jacobian
+# define jacType ls_mem->jacType
+# define g_fun   ls_mem->g_fun
+# define ng_fun  ls_mem->ng_fun
+# define jroot   ls_mem->jroot
+
+/* =============================
+ *
+ *         LSodarMalloc
+ *
+ * =============================
+ *
+ * LSodarMalloc allocates and initializes memory for a problem.
+ * All problem inputs are checked for errors. If any error occurs during initialization,
+ * it is reported to the file whose file pointer is errfp and an error flag is returned.
+ * Otherwise, it returns CV_SUCCESS.
+ */
+
+int LSodarMalloc (void * lsodar_mem, LSRhsFn f, realtype t0, N_Vector y, int itol, realtype reltol, void * abstol)
+{
+    /* Check the input arguments */
+
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarMalloc", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (f == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarMalloc", MSGCV_NULL_F);
+        return (CV_ILL_INPUT);
+    }
+    if (y == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarMalloc", MSGCV_NULL_Y0);
+        return (CV_ILL_INPUT);
+    }
+    if (reltol < 0.)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarMalloc", MSGCV_BAD_RELTOL);
+        return (CV_ILL_INPUT);
+    }
+    if (*((realtype *) abstol) < 0.)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarMalloc", MSGCV_BAD_ABSTOL);
+        return (CV_ILL_INPUT);
+    }
+
+    /* Copy the arguments into the problem memory space */
+    func   =  f;
+    yVec   =  NV_DATA_S(y);
+    tStart =  t0;
+    iTol   =  itol;
+    relTol =  reltol;
+    absTol =  *((realtype *) abstol);
+
+    /* Allocate rwork and iwork workspaces and set them to zero.
+       Their size is lrw and liw, respectively */
+    rwork = (struct rWork_t *) calloc(lrw, sizeof(realtype));
+    iwork = calloc(liw, sizeof(int));
+
+    return (CV_SUCCESS);
+}
+
+/* =============================
+ *
+ *         LSodarReInit
+ *
+ * =============================
+ *
+ * LSodarReInit re-initializes LSODAR's memory for a problem,
+ * assuming it has already been allocated in a prior LSodarMalloc call.
+ * All problem specification inputs are checked for errors.
+ * If any error occurs during initialization, it is reported to the file whose file pointer is errfp.
+ * The return value is CV_SUCCESS = 0 if no errors occurred, or a negative value otherwise.
+ */
+
+int LSodarReInit (void * lsodar_mem, LSRhsFn f, realtype tOld, N_Vector y, int itol, realtype reltol, void * abstol)
+{
+    double rwork0, rwork5;
+
+    /* Check the input arguments */
+
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarReInit", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (y == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarReInit", MSGCV_NULL_Y0);
+        return (CV_ILL_INPUT);
+    }
+    if (f == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarReInit", MSGCV_NULL_F);
+        return (CV_ILL_INPUT);
+    }
+    if (reltol < 0.)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarReInit", MSGCV_BAD_RELTOL);
+        return (CV_ILL_INPUT);
+    }
+    if (*((realtype *) abstol) < 0.)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarReInit", MSGCV_BAD_ABSTOL);
+        return (CV_ILL_INPUT);
+    }
+
+    /* Reset the problem memory space variables to the arguments */
+    func   =  f;
+    *nEq   =  NV_LENGTH_S(y);
+    yVec   =  NV_DATA_S(y);
+    tStart =  tOld;
+    iTol   =  itol;
+    relTol =  reltol;
+    absTol =  *((realtype *) abstol);
+    iState =  1;
+
+    /* Reinitialize rwork and iwork, leave rwork->tcrit and rwork->hmax unchanged, containing tcrit and hmax */
+    rwork0 = rwork->tcrit;
+    rwork5 = rwork->hmax;
+    memset(rwork, 0, lrw);
+    memset(iwork, 0, liw);
+    rwork->tcrit = rwork0;
+    rwork->hmax = rwork5;
+
+    return (CV_SUCCESS);
+}
+
+/* =============================
+ *
+ *         LSodarRootInit
+ *
+ * =============================
+ *
+ * LSodarRootInit initializes a rootfinding problem to be solved during the integration of the ODE system.
+ * It loads the root function pointer and the number of root functions, and allocates workspace memory.
+ * The return value is CV_SUCCESS = 0 if no errors occurred, or a negative value otherwise.
+ */
+
+int LSodarRootInit (void * lsodar_mem, int ng, LSRootFn g, void *gdata)
+{
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarRootInit", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (g == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodarRootInit", MSGCV_NULL_G);
+        return (CV_MEM_NULL);
+    }
+
+    g_fun = g;
+    ng_fun = ng;
+
+    /* Allocate jroot and set it to zero */
+    if (ng > 0)
+    {
+        jroot = calloc(ng, sizeof(int));
+    }
+
+    return (CV_SUCCESS);
+}
+
+
+/* =============================
+ *
+ *       LSodarSetMaxStep
+ *
+ * =============================
+ *
+ * Sets the maximum step size, stocked in rwork->hmax.
+ * Sets iOpt to 1 for rwork->hmax to be taken in consideration by lsodar().
+ */
+
+int LSodarSetMaxStep (void * lsodar_mem, realtype hMax)
+{
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarSetMaxStep", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (iOpt == 0)
+    {
+        iOpt = 1;
+    }
+    rwork->hmax = hMax;
+
+    return (CV_SUCCESS);
+}
+
+/* =============================
+ *
+ *       LSodarSetStopTime
+ *
+ * =============================
+ *
+ * Specifies the time beyond which the integration is not to proceed, stocked in rwork->tcrit.
+ * Sets iOpt to 1 for rwork->tcrit to be taken in consideration by lsodar().
+ */
+
+int LSodarSetStopTime (void * lsodar_mem, realtype tCrit)
+{
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarSetStopTime", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (iOpt == 0)
+    {
+        iOpt = 1;
+    }
+    rwork->tcrit = tCrit;
+
+    return (CV_SUCCESS);
+}
+
+/* =============================
+ *
+ *            LSodar
+ *
+ * =============================
+ *
+ * This routine is the main driver of LSODAR.
+ *
+ * It integrates and looks for roots over a time interval defined by the user.
+ *
+ * The first time that LSodar is called for a successfully initialized problem,
+ * it computes a tentative initial step size h.
+ *
+ * LSodar supports five modes, specified by itask: LS_NORMAL, LS_ONE_STEP, LS_MESH_STEP, LS_NORMAL_TSTOP, and LS_ONE_STEP_TSTOP.
+ *
+ * In the LS_NORMAL mode, the solver steps until it reaches or passes tout and then interpolates to obtain y(tOut).
+ * In the LS_ONE_STEP mode, it takes one internal step and returns.
+ * LS_MESH_STEP means stop at the first internal mesh point at or beyond t = tout and return.
+ * LS_NORMAL_TSTOP and LS_ONE_STEP_TSTOP are similar to LS_NORMAL and LS_ONE_STEP, respectively,
+ * but the integration never proceeds past tcrit (which must have been defined through a call to LSodarSetStopTime).
+ *
+ * It returns CV_ROOT_RETURN if a root was detected, CV_SUCCESS if the integration went correctly,
+ * or a corresponding error flag.
+ */
+
+int LSodar (void * lsodar_mem, realtype tOut, N_Vector yOut, realtype * tOld, enum iTask_t itask)
+{
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodar", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    if (yOut == NULL)
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodar", MSGCV_YOUT_NULL);
+        return (CV_ILL_INPUT);
+    }
+    if ((itask != LS_NORMAL) && (itask != LS_ONE_STEP) && (itask != LS_MESH_STEP) &&
+            (itask != LS_NORMAL_TSTOP) && (itask != CV_ONE_STEP_TSTOP))
+    {
+        LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodar", MSGCV_BAD_ITASK);
+        return(CV_ILL_INPUT);
+    }
+
+    /* Retrieve nEq if it has changed, use a copy of the solution vector and stock the simulation times */
+    *nEq = NV_LENGTH_S(yOut);
+    yVec = NV_DATA_S(yOut);
+    tStart = *tOld;
+    tEnd = tOut;
+
+    /* Launch the simulation with the memory space parameters.
+       lsodar() will update yVec, iState, rwork, iwork and jroot */
+    C2F(lsodar) (func, nEq, yVec, &tStart, &tEnd, &iTol, &relTol, &absTol, &itask, &iState, &iOpt, rwork, &lrw, iwork, &liw, &jac, &jacType, g_fun, &ng_fun, jroot);
+
+    /* Increment the start times */
+    *tOld = tOut;
+    tStart = tEnd;
+
+    /* lsodar() stocked the completion status in iState; return accordingly  */
+    switch (iState)
+    {
+        case 3:
+            return (CV_ROOT_RETURN);
+        case -1:
+            LSProcessError(ls_mem, CV_TOO_MUCH_WORK, "LSODAR", "LSodar", MSGCV_MAX_STEPS, tStart);
+            return (CV_TOO_MUCH_WORK);
+        case -2:
+            LSProcessError(ls_mem, CV_TOO_MUCH_ACC, "LSODAR", "LSodar", MSGCV_TOO_MUCH_ACC, tStart);
+            return (CV_TOO_MUCH_ACC);
+        case -3:
+            LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodar", MSGCV_BAD_INPUT, tStart);
+            return (CV_ILL_INPUT);
+        case -4:
+            LSProcessError(ls_mem, CV_ERR_FAILURE, "LSODAR", "LSodar", MSGCV_ERR_FAILS, tStart);
+            return (CV_ERR_FAILURE);
+        case -5:
+            LSProcessError(ls_mem, CV_CONV_FAILURE, "LSODAR", "LSodar", MSGCV_CONV_FAILS, tStart);
+            return (CV_CONV_FAILURE);
+        case -6:
+            LSProcessError(ls_mem, CV_ILL_INPUT, "LSODAR", "LSodar", MSGCV_EWT_NOW_BAD, tStart);
+            return (CV_ILL_INPUT);
+        default:
+            return (CV_SUCCESS);
+    }
+}
+
+/* =============================
+ *
+ *       LSodarGetRootInfo
+ *
+ * =============================
+ *
+ * Updates rootsfound[] to the computed roots stocked in jroot[].
+ */
+
+int LSodarGetRootInfo (void * lsodar_mem, int * rootsfound)
+{
+    LSodarMem ls_mem = NULL;
+    if (lsodar_mem == NULL)
+    {
+        LSProcessError(NULL, CV_MEM_NULL, "LSODAR", "LSodarGetRootInfo", MSGCV_NO_MEM);
+        return (CV_MEM_NULL);
+    }
+    ls_mem = (LSodarMem) lsodar_mem;
+
+    /* Copy jroot to rootsfound */
+    memcpy(rootsfound, jroot, ng_fun * sizeof(int));
+
+    return (CV_SUCCESS);
+}
+
+/* =============================
+ *
+ *            LSFree
+ *
+ * =============================
+ *
+ * This routine frees the problem memory allocated by LSodarMalloc.
+ */
+
+void LSodarFree (void ** lsodar_mem)
+{
+    LSodarMem ls_mem = NULL;
+    if (*lsodar_mem == NULL)
+    {
+        return;
+    }
+    ls_mem = (LSodarMem) (*lsodar_mem);
+
+    /* Free the inner vectors */
+    LSFreeVectors (ls_mem);
+
+    free (*lsodar_mem);
+    *lsodar_mem = NULL;
+}
+
+/* =============================
+ *
+ *         LSFreeVectors
+ *
+ * =============================
+ *
+ * Frees the problem memory space vectors.
+ */
+
+void LSFreeVectors (LSodarMem ls_mem)
+{
+    /* rwork, iwork and jroot have been allocated; free them */
+    free (rwork);
+    free (iwork);
+    free (jroot);
+}
+
+#define ehfun    LSErrHandler
+#define eh_data  (void *) lsodar_mem
+
+/* =============================
+ *
+ *         LSProcessError
+ *
+ * =============================
+ *
+ * Error handling function.
+ */
+
+void LSProcessError (LSodarMem lsodar_mem, int error_code, const char *module, const char *fname, const char *msgfmt, ...)
+{
+    va_list ap;
+    char msg[256];
+
+    /* Initialize the argument pointer variable
+       (msgfmt is the last required argument to LSProcessError) */
+    va_start(ap, msgfmt);
+
+    if (lsodar_mem == NULL)      /* We write to stderr */
+    {
+#ifndef NO_FPRINTF_OUTPUT
+        fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);
+        fprintf(stderr, msgfmt);
+        fprintf(stderr, "\n\n");
+#endif
+    }
+    else                     /* We can call ehfun */
+    {
+        /* Compose the message */
+        vsprintf(msg, msgfmt, ap);
+
+        /* Call ehfun */
+        ehfun(error_code, module, fname, msg, eh_data);
+    }
+
+    /* Finalize argument processing */
+    va_end(ap);
+
+    return;
+}
+
+#define errfp    stderr
+
+/* =============================
+ *
+ *          LSErrHandler
+ *
+ * =============================
+ *
+ * Default error handling function.
+ */
+
+void LSErrHandler (int error_code, const char *module, const char *function, char *msg, void *data)
+{
+    char err_type[10];
+
+    if (error_code == CV_WARNING)
+    {
+        sprintf(err_type, "WARNING");
+    }
+    else
+    {
+        sprintf(err_type, "ERROR");
+    }
+
+#ifndef NO_FPRINTF_OUTPUT
+    if (errfp != NULL)
+    {
+        fprintf(errfp, "\n[%s %s]  %s\n", module, err_type, function);
+        fprintf(errfp, "  %s\n\n", msg);
+    }
+#endif
+
+    return;
+}
+
diff --git a/scilab/modules/scicos/src/c/lsodar.h b/scilab/modules/scicos/src/c/lsodar.h
new file mode 100644 (file)
index 0000000..0939c44
--- /dev/null
@@ -0,0 +1,112 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2012 - 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
+ *
+ */
+
+#ifndef _LSODAR_H
+#define _LSODAR_H
+
+#include "sundials_extension.h"
+#include "sundials/sundials_types.h" // Definition of types 'realtype' and 'booleantype'
+#include "nvector/nvector_serial.h"  // Type 'N_Vector'
+#include "../scicos_sundials/src/cvode/cvode_impl.h" // Error handling
+
+#ifndef max
+#define max(A,B) ((A>B) ? A:B)  // 'max()' function
+#endif
+
+// realtype workspace
+struct rWork_t
+{
+    realtype tcrit;
+    realtype rwork2;
+    realtype rwork3;
+    realtype rwork4;
+    realtype h0;
+    realtype hmax;
+    realtype hmin;
+    realtype rwork[1];
+};
+
+// Derivative computation and Root functions
+typedef int (*LSRhsFn) (int * neq, realtype * t, realtype * y, realtype * rwork);
+typedef int (*LSRootFn) (int * neq, realtype * t, realtype * y, int * ng, realtype * rwork);
+
+// Potential tasks
+enum iTask_t
+{
+    LS_NORMAL = 1,
+    LS_ONE_STEP = 2,
+    LS_MESH_STEP = 3,
+    LS_NORMAL_TSTOP = 4,
+    LS_ONE_STEP_TSTOP = 5
+};
+
+// LSodar problem memory structure
+typedef struct LSodarMemRec
+{
+    LSRhsFn func;
+    int * nEquations;
+    realtype * yVector;
+    realtype tStart;
+    realtype tEnd;
+    int iTol;
+    realtype relTol;
+    realtype absTol;
+    int iState;
+    int iOpt;
+    struct rWork_t * rwork;
+    int lrw;
+    int * iwork;
+    int liw;
+    int jacobian;
+    int jacType;
+    LSRootFn g_fun;
+    int ng_fun;
+    int * jroot;
+} *LSodarMem;
+
+// Creating the problem
+void * LSodarCreate (int * neq, int ng);
+
+// Allocating the problem
+int LSodarMalloc (void * lsodar_mem, LSRhsFn f, realtype t0, N_Vector y, int itol, realtype reltol, void * abstol);
+
+// Reinitializing the problem
+int LSodarReInit (void * lsodar_mem, LSRhsFn f, realtype tOld, N_Vector y, int itol, realtype reltol, void * abstol);
+
+// Initializing the root-finding problem
+int LSodarRootInit (void * lsodar_mem, int ng, LSRootFn g, void *gdata);
+
+// Specifying the maximum step size
+int LSodarSetMaxStep (void * lsodar_mem, realtype hmax);
+
+// Specifying the time beyond which the integration is not to proceed
+int LSodarSetStopTime (void * lsodar_mem, realtype tcrit);
+
+// Solving the problem
+int LSodar (void * lsodar_mem, realtype tOut, N_Vector yVec, realtype * tOld, enum iTask_t itask);
+
+// Update rootsfound to the computed jroots
+int LSodarGetRootInfo (void * lsodar_mem, int * rootsfound);
+
+// Freeing the problem memory allocated by lsodarMalloc
+void LSodarFree (void ** lsodar_mem);
+
+// Freeing the lsodar vectors allocated in lsodarAllocVectors
+void LSFreeVectors (LSodarMem lsodar_mem);
+
+// Error handling function
+void LSProcessError (LSodarMem lsodar_mem, int error_code, const char *module, const char *fname, const char *msgfmt, ...);
+
+// Default error handling function
+void LSErrHandler (int error_code, const char *module, const char *function, char *msg, void *data);
+
+#endif
index 27ca045..f84ef7f 100644 (file)
@@ -74,6 +74,8 @@
 #include "sciblk4.h"
 #include "dynlib_scicos.h"
 
+#include "lsodar.h"           /* prototypes for lsodar fcts. and consts. */
+
 #if defined(linux) && defined(__i386__)
 #include "setPrecisionFPU.h"
 #endif
@@ -96,7 +98,7 @@ SCICOS_IMPEXP SCSPTR_struct C2F(scsptr);
 /*--------------------------------------------------------------------------*/
 
 #define freeall                                        \
-       if (*neq>0) CVodeFree(&cvode_mem);              \
+       if (*neq>0) if (C2F(cmsolver).solver) (CVodeFree(&cvode_mem)); else (LSodarFree(&cvode_mem));           \
        if (*neq>0) N_VDestroy_Serial(y);               \
        if ( ng>0 ) FREE(jroot);                        \
        if ( ng>0 ) FREE(zcros);
@@ -227,6 +229,8 @@ static int CallKinsol(double *told);
 static int simblk(realtype t, N_Vector yy, N_Vector yp, void *f_data);
 static int grblkdaskr(realtype t, N_Vector yy, N_Vector yp, realtype *gout, void *g_data);
 static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data);
+static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res);
+static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngroot, realtype * res);
 static void addevs(double t, int *evtnb, int *ierr1);
 static int synchro_g_nev(ScicosImport *scs_imp, double *g, int kf, int *ierr);
 static void Multp(double *A, double *B, double *R, int ra, int rb, int ca, int cb);
@@ -799,19 +803,23 @@ int C2F(scicos)(double *x_in, int *xptr_in, double *z__,
     {
 
         /*     integration */
-        if (C2F(cmsolver).solver == 0)        /*  CVODE: Method: BDF,   Nonlinear solver= NEWTON     */
+        if (C2F(cmsolver).solver == 0)        /*  LSODAR: Method: DYNAMIC, Nonlinear solver= DYNAMIC */
         {
             cossim(t0);
         }
-        else if (C2F(cmsolver).solver == 1)   /*  CVODE: Method: BDF,   Nonlinear solver= FUNCTIONAL */
+        else if (C2F(cmsolver).solver == 1)   /*  CVODE: Method: BDF,   Nonlinear solver= NEWTON     */
         {
             cossim(t0);
         }
-        else if (C2F(cmsolver).solver == 2)   /*  CVODE: Method: ADAMS, Nonlinear solver= NEWTON     */
+        else if (C2F(cmsolver).solver == 2)   /*  CVODE: Method: BDF,   Nonlinear solver= FUNCTIONAL */
         {
             cossim(t0);
         }
-        else if (C2F(cmsolver).solver == 3)   /*  CVODE: Method: ADAMS, Nonlinear solver= FUNCTIONAL */
+        else if (C2F(cmsolver).solver == 3)   /*  CVODE: Method: ADAMS, Nonlinear solver= NEWTON     */
+        {
+            cossim(t0);
+        }
+        else if (C2F(cmsolver).solver == 4)   /*  CVODE: Method: ADAMS, Nonlinear solver= FUNCTIONAL */
         {
             cossim(t0);
         }
@@ -1352,15 +1360,18 @@ static void cossim(double *told)
         switch (C2F(cmsolver).solver)
         {
             case 0:
-                cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
+                cvode_mem = LSodarCreate(neq, ng); /* Create the lsodar problem */
                 break;
             case 1:
-                cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
+                cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
                 break;
             case 2:
-                cvode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
+                cvode_mem = CVodeCreate(CV_BDF, CV_FUNCTIONAL);
                 break;
             case 3:
+                cvode_mem = CVodeCreate(CV_ADAMS, CV_NEWTON);
+                break;
+            case 4:
                 cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
                 break;
         }
@@ -1376,6 +1387,11 @@ static void cossim(double *told)
             return;
         }
 
+        if (!C2F(cmsolver).solver)
+        {
+            flag = LSodarMalloc(cvode_mem, simblklsodar, T0, y, CV_SS, reltol, &abstol);
+        }
+        else
         flag = CVodeMalloc(cvode_mem, simblk, T0, y, CV_SS, reltol, &abstol);
         if (check_flag(&flag, "CVodeMalloc", 1))
         {
@@ -1384,6 +1400,11 @@ static void cossim(double *told)
             return;
         }
 
+        if (!C2F(cmsolver).solver)
+        {
+            flag = LSodarRootInit(cvode_mem, ng, grblklsodar, NULL);
+        }
+        else
         flag = CVodeRootInit(cvode_mem, ng, grblk, NULL);
         if (check_flag(&flag, "CVodeRootInit", 1))
         {
@@ -1392,6 +1413,7 @@ static void cossim(double *told)
             return;
         }
 
+        if (C2F(cmsolver).solver)
         /* Call CVDense to specify the CVDENSE dense linear solver */
         flag = CVDense(cvode_mem, *neq);
         if (check_flag(&flag, "CVDense", 1))
@@ -1403,6 +1425,11 @@ static void cossim(double *told)
 
         if (hmax > 0)
         {
+            if (!C2F(cmsolver).solver)
+            {
+                flag = LSodarSetMaxStep(cvode_mem, (realtype) hmax);
+            }
+            else
             flag = CVodeSetMaxStep(cvode_mem, (realtype) hmax);
             if (check_flag(&flag, "CVodeSetMaxStep", 1))
             {
@@ -1566,6 +1593,11 @@ L30:
 
                 if (hot == 0) /* hot==0 : cold restart*/
                 {
+                    if (!C2F(cmsolver).solver)
+                    {
+                        flag = LSodarSetStopTime(cvode_mem, (realtype)tstop);
+                    }
+                    else
                     flag = CVodeSetStopTime(cvode_mem, (realtype)tstop);  /* Setting the stop time*/
                     if (check_flag(&flag, "CVodeSetStopTime", 1))
                     {
@@ -1574,6 +1606,11 @@ L30:
                         return;
                     }
 
+                    if (!C2F(cmsolver).solver)
+                    {
+                        flag = LSodarReInit(cvode_mem, simblklsodar, (realtype)(*told), y, CV_SS, reltol, &abstol);
+                    }
+                    else
                     flag = CVodeReInit(cvode_mem, simblk, (realtype)(*told), y, CV_SS, reltol, &abstol);
                     if (check_flag(&flag, "CVodeReInit", 1))
                     {
@@ -1620,6 +1657,11 @@ L30:
                 if (Discrete_Jump == 0) /* if there was a dzero, its event should be activated*/
                 {
                     phase = 2;
+                    if (!C2F(cmsolver).solver)
+                    {
+                        flag = LSodar(cvode_mem, t, y, told, LS_NORMAL);
+                    }
+                    else
                     flag = CVode(cvode_mem, t, y, told, CV_NORMAL_TSTOP);
                     if (*ierr != 0)
                     {
@@ -1682,6 +1724,11 @@ L30:
                     hot = 0;
                     if (Discrete_Jump == 0)
                     {
+                        if (!C2F(cmsolver).solver)
+                        {
+                            flagr = LSodarGetRootInfo(cvode_mem, jroot);
+                        }
+                        else
                         flagr = CVodeGetRootInfo(cvode_mem, jroot);
                         if (check_flag(&flagr, "CVodeGetRootInfo", 1))
                         {
@@ -3427,6 +3474,69 @@ static int grblk(realtype t, N_Vector yy, realtype *gout, void *g_data)
     return 0;
 } /* grblk */
 /*--------------------------------------------------------------------------*/
+/* simblklsodar */
+static int simblklsodar(int * nequations, realtype * tOld, realtype * actual, realtype * res)
+{
+    double tx = 0.;
+    int i = 0, nantest = 0;
+
+    tx = (double) *tOld;
+
+    for (i = 0; i < *nequations; ++i) res[i] = 0; /* à la place de "C2F(dset)(neq, &c_b14,xcdot , &c__1);"*/
+    C2F(ierode).iero = 0;
+    *ierr = 0;
+    odoit(&tx, actual, res, res);
+    C2F(ierode).iero = *ierr;
+
+    if (*ierr == 0)
+    {
+        nantest = 0;
+        for (i = 0; i < *nequations; i++) /* NaN checking */
+        {
+            if ((res[i] - res[i] != 0))
+            {
+                sciprint(_("\nWarning: The computing function #%d returns a NaN/Inf"), i);
+                nantest = 1;
+                break;
+            }
+        }
+        if (nantest == 1) return 349; /* recoverable error; */
+    }
+
+    return (abs(*ierr)); /* ierr>0 recoverable error; ierr>0 unrecoverable error; ierr=0: ok*/
+
+} /* simblklsodar */
+/*--------------------------------------------------------------------------*/
+/* grblklsodar */
+static int grblklsodar(int * nequations, realtype * tOld, realtype * actual, int * ngc, realtype * res)
+{
+    double tx = 0.;
+    int jj = 0, nantest = 0;
+
+   tx = (double) *tOld;
+
+    C2F(ierode).iero = 0;
+    *ierr = 0;
+
+    zdoit(&tx, actual, actual, res);
+
+    if (*ierr == 0)
+    {
+        nantest = 0;
+        for (jj = 0; jj < *ngc; jj++)
+            if (res[jj] - res[jj] != 0)
+            {
+                sciprint(_("\nWarning: The zero_crossing function #%d returns a NaN/Inf"), jj);
+                nantest = 1;
+                break;
+            } /* NaN checking */
+        if (nantest == 1) return 350; /* recoverable error; */
+    }
+    C2F(ierode).iero = *ierr;
+
+    return 0;
+} /* grblklsodar */
+/*--------------------------------------------------------------------------*/
 /* simblkdaskr */
 static int simblkdaskr(realtype tres, N_Vector yy, N_Vector yp, N_Vector resval, void *rdata)
 {
index 713d81b..c10a492 100644 (file)
@@ -94,10 +94,11 @@ lib /DEF:"$(ProjectDir)Graphics_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platfo
 lib /DEF:"$(ProjectDir)Scicos_blocks_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scicos_blocks_f.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)Localization_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scilocalization.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL
-lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.lib" 1&gt;NUL 2&gt;NUL</Command>
+lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.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>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <IgnoreSpecificDefaultLibraries>
       </IgnoreSpecificDefaultLibraries>
@@ -136,10 +137,11 @@ lib /DEF:"$(ProjectDir)Graphics_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platfo
 lib /DEF:"$(ProjectDir)Scicos_blocks_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scicos_blocks_f.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)Localization_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scilocalization.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL
-lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.lib" 1&gt;NUL 2&gt;NUL</Command>
+lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.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>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <IgnoreSpecificDefaultLibraries>
       </IgnoreSpecificDefaultLibraries>
@@ -178,10 +180,11 @@ lib /DEF:"$(ProjectDir)Graphics_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platfo
 lib /DEF:"$(ProjectDir)Scicos_blocks_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scicos_blocks_f.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)Localization_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scilocalization.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL
-lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.lib" 1&gt;NUL 2&gt;NUL</Command>
+lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.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>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <IgnoreSpecificDefaultLibraries>
       </IgnoreSpecificDefaultLibraries>
@@ -225,10 +228,11 @@ lib /DEF:"$(ProjectDir)Graphics_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platfo
 lib /DEF:"$(ProjectDir)Scicos_blocks_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scicos_blocks_f.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)Localization_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)scilocalization.lib" 1&gt;NUL 2&gt;NUL
 lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)elementary_functions_f.lib" 1&gt;NUL 2&gt;NUL
-lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.lib" 1&gt;NUL 2&gt;NUL</Command>
+lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform) /OUT:"$(ProjectDir)core_f.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>
     </PreLinkEvent>
     <Link>
-      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>../../../../bin/blasplus.lib;core.lib;output_stream.lib;string.lib;dynamic_link.lib;graphics.lib;scicos_blocks_f.lib;scicos_f.lib;scilocalization.lib;elementary_functions_f.lib;core_f.lib;differential_equations_f.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <IgnoreSpecificDefaultLibraries>
       </IgnoreSpecificDefaultLibraries>
@@ -257,6 +261,7 @@ lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform
     <ClCompile Include="il_state.c" />
     <ClCompile Include="import.c" />
     <ClCompile Include="MlistGetFieldNumber.c" />
+    <ClCompile Include="lsodar.c" />
     <ClCompile Include="..\..\sci_gateway\c\sci_buildouttb.c" />
     <ClCompile Include="..\..\sci_gateway\c\sci_callblk.c" />
     <ClCompile Include="..\..\sci_gateway\c\sci_coserror.c" />
@@ -309,6 +314,7 @@ lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform
     <ClInclude Include="il_state.h" />
     <ClInclude Include="..\..\includes\import.h" />
     <ClInclude Include="MlistGetFieldNumber.h" />
+    <ClInclude Include="lsodar.h" />
     <ClInclude Include="..\..\includes\scicos-def.h" />
     <ClInclude Include="..\..\includes\scicos.h" />
     <ClInclude Include="..\..\includes\scicos_free.h" />
@@ -322,6 +328,7 @@ lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform
   <ItemGroup>
     <None Include="..\..\locales\scicos.pot" />
     <None Include="Core_f_Import.def" />
+    <None Include="differential_equations_f_Import.def" />
     <None Include="Dynamic_link_Import.def" />
     <None Include="elementary_functions_f_Import.def" />
     <None Include="Graphics_Import.def" />
@@ -375,4 +382,4 @@ lib /DEF:"$(ProjectDir)core_f_Import.def" /SUBSYSTEM:WINDOWS /MACHINE:$(Platform
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
   </ImportGroup>
-</Project>
\ No newline at end of file
+</Project>
index 068b9f4..26fce24 100644 (file)
@@ -56,6 +56,9 @@
     <ClCompile Include="MlistGetFieldNumber.c">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="lsodar.c">
+      <Filter>Source Files</Filter>
+    </ClCompile>
     <ClCompile Include="..\..\sci_gateway\c\sci_buildouttb.c">
       <Filter>Source Files</Filter>
     </ClCompile>
     <ClInclude Include="MlistGetFieldNumber.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="lsodar.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
     <ClInclude Include="..\..\includes\scicos-def.h">
       <Filter>Header Files</Filter>
     </ClInclude>
     <None Include="Core_f_Import.def">
       <Filter>Libraries Dependencies\Imports</Filter>
     </None>
+    <None Include="differential_equations_f_Import.def">
+      <Filter>Libraries Dependencies\Imports</Filter>
+    </None>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="scicos.rc">
       <Filter>Resource Files</Filter>
     </ResourceCompile>
   </ItemGroup>
-</Project>
\ No newline at end of file
+</Project>
index dbed09f..499ffe0 100644 (file)
@@ -473,6 +473,7 @@ extern "C" {
 
   /* CVode Error Messages */
 
+#define MSGCV_BAD_INPUT "One of the arguments is illegal"
 #define MSGCV_LSOLVE_NULL "The linear solver's solve routine is NULL."
 #define MSGCV_YOUT_NULL "yout = NULL illegal."
 #define MSGCV_TRET_NULL "tret = NULL illegal."
diff --git a/scilab/modules/scicos/tests/unit_tests/LSodar.dia.ref b/scilab/modules/scicos/tests/unit_tests/LSodar.dia.ref
new file mode 100644 (file)
index 0000000..83584fa
--- /dev/null
@@ -0,0 +1,32 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- ENGLISH IMPOSED -->
+// Execute with exec("SCI/modules/scicos/tests/unit_tests/LSodar.tst");
+//  or test_run('scicos', 'LSodar', ['no_check_error_output']);
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/LSodar_test.zcos"));
+// Set solver to LSodar + run LSodar + save results
+scs_m.props.tol(6) = 0;                                                                                        // Set solver to LSodar
+try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;        // Run LSodar
+lsodarval = res.values;                // Results
+time = res.time;                       // Time
+// Set solver to CVode BDF/Newton + run + save results
+scs_m.props.tol(6) = 1;
+try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+cvval = res.values;
+// Compare results
+compa = abs(lsodarval-cvval);
+// Extract mean, standard deviation, maximum
+mea = mean(compa);
+[maxi, indexMaxi] = max(compa);
+stdeviation = st_deviation(compa);
+// Verifying closeness of the results
+assert_checktrue(maxi <= 10^-8);
+assert_checktrue(mea <= 10^-8);
+assert_checktrue(stdeviation <= 10^-8);
diff --git a/scilab/modules/scicos/tests/unit_tests/LSodar.tst b/scilab/modules/scicos/tests/unit_tests/LSodar.tst
new file mode 100644 (file)
index 0000000..39d2539
--- /dev/null
@@ -0,0 +1,40 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- ENGLISH IMPOSED -->
+
+// Execute with exec("SCI/modules/scicos/tests/unit_tests/LSodar.tst");
+//  or test_run('scicos', 'LSodar', ['no_check_error_output']);
+
+// Import diagram
+loadScicos();
+loadXcosLibs();
+assert_checktrue(importXcosDiagram("SCI/modules/xcos/tests/unit_tests/LSodar_test.zcos"));
+
+// Set solver to LSodar + run LSodar + save results
+scs_m.props.tol(6) = 0;                                         // Set solver to LSodar
+try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end; // Run LSodar
+lsodarval = res.values;  // Results
+time = res.time;         // Time
+
+// Set solver to CVode BDF/Newton + run + save results
+scs_m.props.tol(6) = 1;
+try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+cvval = res.values;
+
+// Compare results
+compa = abs(lsodarval-cvval);
+
+// Extract mean, standard deviation, maximum
+mea = mean(compa);
+[maxi, indexMaxi] = max(compa);
+stdeviation = st_deviation(compa);
+
+// Verifying closeness of the results
+assert_checktrue(maxi <= 10^-8);
+assert_checktrue(mea <= 10^-8);
+assert_checktrue(stdeviation <= 10^-8);
diff --git a/scilab/modules/xcos/examples/solvers/integLSodar.sce b/scilab/modules/xcos/examples/solvers/integLSodar.sce
new file mode 100644 (file)
index 0000000..b2dd54f
--- /dev/null
@@ -0,0 +1,31 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - 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
+
+// Run with exec("SCI/modules/xcos/examples/solvers/integLSodar.sce");
+
+// Import the diagram and augment the ending time
+loadScicos();
+loadXcosLibs();
+importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
+scs_m.props.tf = 30000;
+
+solverName = ["LSodar", "BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional"];
+
+for solver = 0:4
+
+ // Select the solver
+ scs_m.props.tol(6) = solver;
+
+ // Start the timer, launch the simulation and display time
+ tic();
+ try scicos_simulate(scs_m, 'nw'); catch disp(lasterror()); end;
+ t = toc();
+ disp(t, "Time for " + solverName(solver+1) + " :");
+
+end
index 44fbb4e..a8e93f3 100644 (file)
@@ -15,15 +15,15 @@ loadXcosLibs();
 importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
 scs_m.props.tf = 30000;
 
-solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functionnal", "Runge-Kutta"];
+solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional", "Runge-Kutta"];
 
 for solver=0:4
 
  // Select the solver
- scs_m.props.tol(6) = solver;
+ scs_m.props.tol(6) = solver+1;
 
  // Set max step size if Runge-Kutta
- if (solver == 4) scs_m.props.tol(7) = 0.01;
+ if ((solver+1) == 5) scs_m.props.tol(7) = 0.01;
 
  // Start the timer, launch the simulation and display time
  tic();
@@ -15,12 +15,12 @@ loadXcosLibs();
 importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
 scs_m.props.tf = 30000;
 
-solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functionnal"];
+solverName=["BDF/Newton", "BDF/Functional", "Adams/Newton", "Adams/Functional"];
 
 for solver=0:3
 
  // Select the solver
- scs_m.props.tol(6) = solver;
+ scs_m.props.tol(6) = solver+1;
 
  // Start the timer, launch the simulation and display time
  tic();
index 3e7c5ac..dac1a61 100644 (file)
 loadScicos();
 loadXcosLibs();
 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
+scs_m.props.tol(6) = 1;
 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
 ]]></scilab:image>
         </para>
@@ -206,7 +207,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       scs_m.props.tf = 10000;
 
       // Select the solver BDF / Newton
-      scs_m.props.tol(6) = 0;
+      scs_m.props.tol(6) = 1;
 
       // Start the timer, launch the simulation and display time
       tic();
@@ -255,6 +256,9 @@ Time for Adams / Functional :
                 <link linkend="IDA">IDA</link>
             </member>
             <member>
+                <link linkend="LSodar">LSodar</link>
+            </member>
+            <member>
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
index 0bdaa98..b814b56 100644 (file)
@@ -164,6 +164,9 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
                 <link linkend="CVode">CVode</link>
             </member>
             <member>
+                <link linkend="LSodar">LSodar</link>
+            </member>
+            <member>
                 <link linkend="RK">Runge-Kutta 4(5)</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/help/en_US/solvers/LSodar.xml b/scilab/modules/xcos/help/en_US/solvers/LSodar.xml
new file mode 100644 (file)
index 0000000..f223379
--- /dev/null
@@ -0,0 +1,205 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) Scilab Enterprises - 2012 - 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:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en_US" xml:id="LSodar">
+    <refnamediv>
+        <refname>LSodar</refname>
+        <refpurpose>
+            LSODAR (short for Livermore Solver for Ordinary Differential equations, with Automatic method switching for stiff and nonstiff problems, and with Root-finding) is a numerical solver providing an efficient and stable method to solve Ordinary Differential Equations (ODEs) Initial Value Problems. Called by <link linkend="xcos">xcos</link>.
+        </refpurpose>
+    </refnamediv>
+    <refsection>
+        <title>Description</title>
+        <para>
+            LSODAR (short for Livermore Solver for Ordinary Differential equations, with Automatic method switching for stiff and nonstiff problems, and with Root-finding) is a numerical solver providing an efficient and stable variable-size step method to solve Initial Value Problems of the form :
+        </para>
+        <para>
+            <latex>
+                \begin{eqnarray}
+                \dot{y} = f(t,y), \hspace{3 mm} y(t_0) = y_0, \hspace{3 mm} y \in R^N
+                \end{eqnarray}
+            </latex>
+        </para>
+        <para>
+            LSodar is similar to <emphasis>CVode</emphasis> in many ways :
+            <itemizedlist>
+                <listitem>
+                    It uses variable-size steps,
+                </listitem>
+                <listitem>
+                    It can potentially use <emphasis>BDF</emphasis> and <emphasis>Adams</emphasis> integration, methods,
+                </listitem>
+                <listitem>
+                    <emphasis>BDF</emphasis> and <emphasis>Adams</emphasis> being implicit stable methods, <emphasis>LSodar</emphasis> is suitable for stiff and nonstiff problems,
+                </listitem>
+                <listitem>
+                    They both look for roots over the integration interval.
+                </listitem>
+            </itemizedlist>
+        </para>
+        <para>
+            The main difference though is that <emphasis>LSodar</emphasis> is <emphasis role="bold">fully automated</emphasis>, and chooses between <emphasis>BDF</emphasis> and <emphasis>Adams</emphasis> itself, by checking for stiffness at every step.
+        </para>
+        <para>
+            If the step is considered stiff, then <emphasis>BDF</emphasis> (with max order set to 5) is used and the Modified Newton method <emphasis>'Chord'</emphasis> iteration is selected.
+        </para>
+        <para>
+            Otherwise, the program uses <emphasis>Adams</emphasis> integration (with max order set to 12) and <emphasis>Functional</emphasis> iterations.
+        </para>
+        <para>
+            The stiffness detection is done by step size attempts with both methods.
+        </para>
+        <para>
+            First, if we are in <emphasis>Adams</emphasis> mode and the order is greater than 5, then we assume the problem is nonstiff and proceed with <emphasis>Adams</emphasis>.
+        </para>
+        <para>
+            The first twenty steps use <emphasis>Adams / Functional</emphasis> method.
+            Then <emphasis>LSodar</emphasis> computes the ideal step size of both methods. If the step size advantage is at least <emphasis>ratio = 5</emphasis>, then the current method switches (<emphasis>Adams / Functional</emphasis> to <emphasis>BDF / Chord Newton</emphasis> or vice versa).
+        </para>
+        <para>
+            After every switch, <emphasis>LSodar</emphasis> takes twenty steps, then starts comparing the step sizes at every step.
+        </para>
+        <para>
+            Such strategy induces a minor overhead computational cost if the problem stiffness is known, but is very effective on problems that require differentiate precision. For instance, discontinuities-sensitive problems.
+        </para>
+        <para>
+            Concerning precision, the two integration/iteration methods being close to <emphasis>CVode</emphasis>'s, the results are very similar.
+        </para>
+    </refsection>
+    <refsection>
+        <title>Examples</title>
+        <para>
+            <link type="scilab" linkend="scilab.xcos/xcos/examples/solvers/ODE_Example.xcos">
+                <inlinemediaobject>
+                    <imageobject>
+                        <imagedata align="center" fileref="../../../examples/solvers/ODE_Example.xcos" valign="middle"/>
+                    </imageobject>
+                </inlinemediaobject>
+            </link>
+            <scilab:image><![CDATA[
+loadScicos();
+loadXcosLibs();
+importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
+scs_m.props.tol(6) = 0;
+try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+]]></scilab:image>
+        </para>
+        <para>
+            The integral block returns its continuous state, we can evaluate it with LSodar by running the example :
+        </para>
+        <para>
+            <programlisting language="example"><![CDATA[
+      // Import the diagram and set the ending time
+      loadScicos();
+      loadXcosLibs();
+      importXcosDiagram("SCI/modules/xcos/examples/solvers/ODE_Example.xcos");
+      scs_m.props.tf = 10000;
+
+      // Select the LSodar solver
+      scs_m.props.tol(6) = 0;
+
+      // Start the timer, launch the simulation and display time
+      tic();
+      try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
+      t = toc();
+      disp(t, "Time for LSodar :");
+      ]]></programlisting>
+        </para>
+        <para>
+            The Scilab console displays :
+            <screen><![CDATA[
+Time for LSodar :
+ 30.584
+            ]]></screen>
+        </para>
+        <para>
+            Now, in the following script, we compare the time difference between the methods by running the example with the five solvers in turn :
+            <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integLSodar.sce">
+                Open the script
+            </link>
+        </para>
+        <para>
+            <screen><![CDATA[
+Time for LSodar :
+ 5.001
+
+Time for BDF / Newton :
+ 12.943
+
+Time for BDF / Functional :
+ 12.415
+
+Time for Adams / Functional :
+ 8.651
+
+Time for Adams / Newton :
+ 7.962
+            ]]></screen>
+        </para>
+        <para>
+            These results show that on a nonstiff problem, for the same precision required, LSodar is significantly faster. Other tests prove the proximity of the results. Indeed, we find that the solution difference order between LSodar and CVode is close to the order of the highest tolerance ( 
+            <emphasis>
+                y<subscript>lsodar</subscript> - y<subscript>cvode</subscript>
+            </emphasis>
+            &#8776; <emphasis>max(reltol, abstol)</emphasis> ).
+        </para>
+        <para>
+            Variable step-size ODE solvers are not appropriate for deterministic real-time applications because the computational overhead of taking a time step varies over the course of an application.
+        </para>
+    </refsection>
+    <refsection>
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="CVode">CVode</link>
+            </member>
+            <member>
+                <link linkend="IDA">IDA</link>
+            </member>
+            <member>
+                <link linkend="RK">Runge-Kutta 4(5)</link>
+            </member>
+            <member>
+                <link linkend="ode">ode</link>
+            </member>
+            <member>
+                <link linkend="ode_discrete">ode_discrete</link>
+            </member>
+            <member>
+                <link linkend="ode_root">ode_root</link>
+            </member>
+            <member>
+                <link linkend="odedc">odedc</link>
+            </member>
+            <member>
+                <link linkend="impl">impl</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>Bibliography</title>
+        <para>
+            ACM SIGNUM Newsletter, Volume 15, Issue 4, December 1980, Pages 10-11 <ulink url="http://dl.acm.org/citation.cfm?id=1218052.1218054&amp;coll=DL&amp;dl=GUIDE&amp;CFID=213153447&amp;CFTOKEN=26634355">LSode - LSodi</ulink>
+        </para>
+        <para>
+            <ulink url="https://computation.llnl.gov/casc/sundials/documentation/documentation.html">Sundials Documentation</ulink>
+        </para>
+    </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.4.1</revnumber>
+                <revdescription>LSodar solver added</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 939a8c2..39c2b22 100644 (file)
             <emphasis>
                 10<superscript>-3</superscript>
             </emphasis>
-            , rounding errors sometimes come into play as we approach 
-            <emphasis>
+            , rounding errors sometimes come into play as we approach <emphasis>
                 4*10<superscript>-4</superscript>
             </emphasis>
             . Indeed, the interval splitting cannot be done properly and we get capricious results.
 loadScicos();
 loadXcosLibs();
 importXcosDiagram(SCI + "/modules/xcos/examples/solvers/ODE_Example.xcos");
-scs_m.props.tol(6) = 4;
+scs_m.props.tol(6) = 5;
 scs_m.props.tol(7) = 10^-2;
 try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
 ]]></scilab:image>
@@ -165,7 +164,7 @@ try xcos_simulate(scs_m, 4); catch disp(lasterror()); end;
       scs_m.props.tf = 10000;
 
       // Select the solver Runge-Kutta and set the precision
-      scs_m.props.tol(6) = 4;
+      scs_m.props.tol(6) = 5;
       scs_m.props.tol(7) = 10^-2;
 
       // Start the timer, launch the simulation and display time
@@ -183,7 +182,7 @@ Time for Runge-Kutta :
             ]]></screen>
         </para>
         <para>
-            Now, in the following script, we compare the time difference between the methods by running the example with the five solvers in turn :
+            Now, in the following script, we compare the time difference between Runge-Kutta and Sundials by running the example with the five solvers in turn :
             <link type="scilab" linkend ="scilab.scinotes/xcos/examples/solvers/integRK.sce">
                 Open the script
             </link>
@@ -223,6 +222,9 @@ Time for Runge-Kutta :
                 <link linkend="IDA">IDA</link>
             </member>
             <member>
+                <link linkend="LSodar">LSodar</link>
+            </member>
+            <member>
                 <link linkend="ode">ode</link>
             </member>
             <member>
diff --git a/scilab/modules/xcos/tests/unit_tests/LSodar_test.zcos b/scilab/modules/xcos/tests/unit_tests/LSodar_test.zcos
new file mode 100644 (file)
index 0000000..e464047
Binary files /dev/null and b/scilab/modules/xcos/tests/unit_tests/LSodar_test.zcos differ