fix non regression tests in differential_equations module.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_ode.cpp
index 093f990..d7234c9 100644 (file)
@@ -2,11 +2,14 @@
 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 * Copyright (C) 2011 - DIGITEO - 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
+ * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
 *
 */
 /*--------------------------------------------------------------------------*/
@@ -20,6 +23,7 @@
 #include "differentialequationfunctions.hxx"
 #include "runvisitor.hxx"
 #include "context.hxx"
+#include "checkodeerror.hxx"
 
 extern "C"
 {
@@ -27,10 +31,9 @@ extern "C"
 #include "Scierror.h"
 #include "scifunctions.h"
 #include "elem_common.h"
-#include "warningmode.h"
+#include "configvariable_interface.h"
 #include "sciprint.h"
 #include "common_structure.h"
-#include "checkodeerror.h"
 #include "sci_malloc.h"
 }
 /*--------------------------------------------------------------------------*/
@@ -39,7 +42,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
@@ -85,6 +88,10 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
     // For root methode
     int* jroot = NULL;
 
+    // error message catched
+    std::wostringstream os;
+    bool bCatch = false;
+
     // *** check the minimal number of input args. ***
     if (in.size() < 4)
     {
@@ -233,10 +240,10 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
     pDblT = in[iPos]->getAs<types::Double>();
 
     // get next inputs
-    DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"ode");
-    DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
-    deFunctionsManager->setOdeYRows(pDblY0->getRows());
-    deFunctionsManager->setOdeYCols(pDblY0->getCols());
+    DifferentialEquationFunctions deFunctionsManager(L"ode");
+    DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
+    deFunctionsManager.setOdeYRows(pDblY0->getRows());
+    deFunctionsManager.setOdeYCols(pDblY0->getCols());
 
     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
     *YSize = pDblY0->getSize();
@@ -336,17 +343,17 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
             if (bFuncF == false)
             {
-                deFunctionsManager->setFFunction(pCall);
+                deFunctionsManager.setFFunction(pCall);
                 bFuncF = true;
             }
             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
             {
-                deFunctionsManager->setJacFunction(pCall);
+                deFunctionsManager.setJacFunction(pCall);
                 bFuncJac = true;
             }
             else if (bFuncG == false && meth == 3)
             {
-                deFunctionsManager->setGFunction(pCall);
+                deFunctionsManager.setGFunction(pCall);
                 bFuncG = true;
             }
             else
@@ -365,17 +372,17 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
 
             if (bFuncF == false)
             {
-                bOK = deFunctionsManager->setFFunction(pStr);
+                bOK = deFunctionsManager.setFFunction(pStr);
                 bFuncF = true;
             }
             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
             {
-                bOK = deFunctionsManager->setJacFunction(pStr);
+                bOK = deFunctionsManager.setJacFunction(pStr);
                 bFuncJac = true;
             }
             else if (bFuncG == false && meth == 3)
             {
-                bOK = deFunctionsManager->setGFunction(pStr);
+                bOK = deFunctionsManager.setGFunction(pStr);
                 bFuncG = true;
             }
             else
@@ -428,13 +435,13 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 if (bFuncF == false)
                 {
                     bFuncF = true;
-                    bOK = deFunctionsManager->setFFunction(pStr);
+                    bOK = deFunctionsManager.setFFunction(pStr);
                     sizeOfpdYData = *YSize;
                 }
                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
                 {
                     bFuncJac = true;
-                    bOK = deFunctionsManager->setJacFunction(pStr);
+                    bOK = deFunctionsManager.setJacFunction(pStr);
                     if (sizeOfpdYData == 0)
                     {
                         sizeOfpdYData = *YSize;
@@ -443,7 +450,7 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 else if (bFuncG == false && meth == 3)
                 {
                     bFuncG = true;
-                    bOK = deFunctionsManager->setGFunction(pStr);
+                    bOK = deFunctionsManager.setGFunction(pStr);
                     if (sizeOfpdYData == 0)
                     {
                         sizeOfpdYData = *YSize;
@@ -505,28 +512,28 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 if (bFuncF == false)
                 {
                     bFuncF = true;
-                    deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
+                    deFunctionsManager.setFFunction(pList->get(0)->getAs<types::Callable>());
                     for (int iter = 1; iter < pList->getSize(); iter++)
                     {
-                        deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
+                        deFunctionsManager.setFArgs(pList->get(iter)->getAs<types::InternalType>());
                     }
                 }
                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
                 {
                     bFuncJac = true;
-                    deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
+                    deFunctionsManager.setJacFunction(pList->get(0)->getAs<types::Callable>());
                     for (int iter = 1; iter < pList->getSize(); iter++)
                     {
-                        deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
+                        deFunctionsManager.setJacArgs(pList->get(iter)->getAs<types::InternalType>());
                     }
                 }
                 else if (bFuncG == false && meth == 3)
                 {
                     bFuncG = true;
-                    deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
+                    deFunctionsManager.setGFunction(pList->get(0)->getAs<types::Callable>());
                     for (int iter = 1; iter < pList->getSize(); iter++)
                     {
-                        deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
+                        deFunctionsManager.setGArgs(pList->get(iter)->getAs<types::InternalType>());
                     }
                 }
             }
@@ -739,7 +746,7 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
     {
         if (getWarningMode())
         {
-            sciprint(_("%s: Warning: Wrong value for maximun stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), "ode", 12, 5);
+            sciprint(_("%s: Warning: Wrong value for maximum stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), "ode", 12, 5);
         }
 
         mxordn = 12;
@@ -800,8 +807,8 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
             lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
             lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
 
-            rworkSize   = max(lrn, lrs);
-            iworkSize   = 20 + *YSize;
+            rworkSize = std::max(lrn, lrs);
+            iworkSize = 20 + *YSize;
 
             dStructTabSize += 241;
             iStructTabSize += 50;
@@ -938,13 +945,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
             }
         }
     }
@@ -1073,12 +1102,10 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
                 }
             }
-            catch (ast::ScilabError &e)
+            catch (ast::InternalError &ie)
             {
-                char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
-                sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
-                Scierror(999, pstrMsg);
-                FREE(pstrMsg);
+                os << ie.GetErrorMessage();
+                bCatch = true;
                 err = 1;
             }
 
@@ -1110,6 +1137,17 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 {
                     FREE(rtol);
                 }
+
+                if (bCatch)
+                {
+                                       char sError[bsiz];
+                                       os_sprintf(sError, "%ls: An error occured in '%s' subroutine.\n", L"ode", strMeth.c_str());
+                                       wchar_t* szError = to_wide_string(sError);
+                    os << szError;
+                    throw ast::InternalError(os.str());
+                                       FREE(szError);
+                }
+
                 return types::Function::Error;
             }
 
@@ -1120,23 +1158,26 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
             {
                 if (getWarningMode())
                 {
-                    sciprint(_("Integration was stoped at t = %lf.\n"), t0);
+                    sciprint(_("Integration was stopped at t = %lf.\n"), t0);
                 }
                 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;
             }
         }
@@ -1150,11 +1191,11 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                     // first lime is the time of result stored in second line.
                     int size = 2 * iSizeList;
                     int* iRanks = (int*)MALLOC(size * sizeof(int));
-
+                    int iMaxRank = pPolyY0->getMaxRank();
                     for(int i = 0; i < iSizeList; i++)
                     {
                         iRanks[i * 2] = 1; // time rank
-                        iRanks[i * 2 + 1] = pPolyY0->getMaxRank(); // result rank
+                        iRanks[i * 2 + 1] = iMaxRank; // result rank
                     }
 
                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
@@ -1200,9 +1241,10 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 {
                     int size = pDblT->getSize();
                     int* iRanks = (int*)MALLOC(size * sizeof(int));
+                    int iMaxRank = pPolyY0->getMaxRank();
                     for(int i = 0; i < size; i++)
                     {
-                        iRanks[i] = pPolyY0->getMaxRank();
+                        iRanks[i] = iMaxRank;
                     }
 
                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
@@ -1212,7 +1254,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++)
@@ -1283,15 +1325,12 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
                 }
             }
-            catch (ast::ScilabError &e)
+            catch (ast::InternalError &ie)
             {
-                char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
-                sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
-                Scierror(999, pstrMsg);
-                FREE(pstrMsg);
+                os << ie.GetErrorMessage();
+                bCatch = true;
                 err = 1;
             }
-
             // FREE allocated data
             if (err == 1) // error case
             {
@@ -1320,6 +1359,17 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
                 {
                     FREE(rtol);
                 }
+
+                if (bCatch)
+                {
+                                       char sError[bsiz];
+                                       os_sprintf(sError, "%ls: An error occured in '%s' subroutine.\n", L"ode", strMeth.c_str());
+                                       wchar_t* szError = to_wide_string(sError);
+                    os << szError;
+                    throw ast::InternalError(os.str());
+                                       FREE(szError);
+                }
+
                 return types::Function::Error;
             }
 
@@ -1344,7 +1394,7 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
             {
                 if (getWarningMode())
                 {
-                    sciprint(_("Integration was stoped at t = %lf.\n"), t0);
+                    sciprint(_("Integration was stopped at t = %lf.\n"), t0);
                 }
 
                 types::Double* pDblYOutTemp = pDblYOut;
@@ -1358,18 +1408,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);
@@ -1465,15 +1518,23 @@ types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, type
             }
         }
 
-        types::Double* pDblRd = new types::Double(1, sizeOfRd);
-        //rd: The first entry contains the stopping time.
-        pDblRd->set(0, C2F(lsr001).tlast);
-        for (int i = 0; i < pDblNg->get(0); i++)
+        types::Double* pDblRd = NULL;
+        if (sizeOfRd == 1) // Not root found, return empty matrix
         {
-            if (jroot[i])
+            pDblRd = types::Double::Empty();
+        }
+        else
+        {
+            pDblRd = new types::Double(1, sizeOfRd);
+            //rd: The first entry contains the stopping time.
+            pDblRd->set(0, C2F(lsr001).tlast);
+            for (int i = 0; i < pDblNg->get(0); i++)
             {
-                k++;
-                pDblRd->set(k, (double)i + 1);
+                if (jroot[i])
+                {
+                    k++;
+                    pDblRd->set(k, (double)i + 1);
+                }
             }
         }
         out.push_back(pDblRd); // rd