scilab/modules/differential_equations/sci_gateway/cpp/sci_ode.cpp

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
   4 *
   5 * This file must be used under the terms of the CeCILL.
   6 * This source file is licensed as described in the file COPYING, which
   7 * you should have received as part of this distribution.  The terms
   8 * are also available at
   9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  10 *
  11 */
  12 /*--------------------------------------------------------------------------*/
  13
  14 #include "differential_equations_gw.hxx"
  15 #include "function.hxx"
  16 #include "double.hxx"
  17 #include "string.hxx"
  18 #include "list.hxx"
  19 #include "callable.hxx"
  20 #include "differentialequationfunctions.hxx"
  21 #include "runvisitor.hxx"
  22 #include "context.hxx"
  23 #include "checkodeerror.hxx"
  24
  25 extern "C"
  26 {
  27 #include "localization.h"
  28 #include "Scierror.h"
  29 #include "scifunctions.h"
  30 #include "elem_common.h"
  31 #include "configvariable_interface.h"
  32 #include "sciprint.h"
  33 #include "common_structure.h"
  34 #include "sci_malloc.h"
  35 }
  36 /*--------------------------------------------------------------------------*/
  37
  38 types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, types::typed_list &out)
  39 {
  40     // Methode
  41     types::String* pStrType     = NULL;
  42     const wchar_t* wcsType      = L"lsoda";
  43     int meth                    = 0;
  44
  45     // y0
  46     //    types::Polynom* pPolyY0     = NULL;
  47     types::Double* pDblY0       = NULL;
  48     double* pdYData             = NULL; // contain y0 following by all args data in list case.
  49     int sizeOfpdYData           = 0;
  50
  51     // Other input args
  52     types::Double* pDblT0       = NULL;
  53     types::Double* pDblT        = NULL;
  54     types::Double* pDblRtol     = NULL;
  55     types::Double* pDblAtol     = NULL;
  56     types::Double* pDblNg       = NULL;
  57     types::Double* pDblW        = NULL;
  58     types::Double* pDblIw       = NULL;
  59
  60     // %ODEOPTIONS
  61     types::Double* pDblOdeOptions = NULL;
  62
  63     // Result
  64     types::Double* pDblYOut     = NULL;
  65     //    types::Double* pDblTOut     = NULL;
  66     //    types::Polynom* pPolyYOut   = NULL;
  67
  68     // Indicate if the function is given.
  69     bool bFuncF     = false;
  70     bool bFuncJac   = false;
  71     bool bFuncG     = false;
  72
  73     int iPos        = 0; // Position in types::typed_list in
  74
  75     int sizeOfYSize = 1;
  76     int* YSize      = NULL;    // YSize(1) = size of y0,
  77     // YSize(n) = size of Args(n) in list case.
  78
  79     C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
  80     // 1 means print, 0 means no printing.
  81     C2F(eh0001).lunit = 6;  // 6 = stdout
  82
  83     int one = 1; // used in C2F(dcopy)
  84
  85     // For root methode
  86     int* jroot = NULL;
  87
  88     // *** check the minimal number of input args. ***
  89     if (in.size() < 4)
  90     {
  91         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 4);
  92         return types::Function::Error;
  93     }
  94
  95     // *** Get the methode. ***
  96     if (in[0]->isString())
  97     {
  98         pStrType = in[0]->getAs<types::String>();
  99         wcsType = pStrType->get(0);
 100         iPos++;
 101     }
 102
 103     if (iPos)
 104     {
 105         if (wcscmp(wcsType, L"adams") == 0)
 106         {
 107             meth = 1;
 108         }
 109         else if (wcscmp(wcsType, L"stiff") == 0)
 110         {
 111             meth = 2;
 112         }
 113         else if (wcscmp(wcsType, L"root") == 0)
 114         {
 115             meth = 3;
 116         }
 117         else if (wcscmp(wcsType, L"discrete") == 0)
 118         {
 119             meth = 4;
 120         }
 121         else if (wcscmp(wcsType, L"rk") == 0)
 122         {
 123             meth = 5;
 124         }
 125         else if (wcscmp(wcsType, L"rkf") == 0)
 126         {
 127             meth = 6;
 128         }
 129         else if (wcscmp(wcsType, L"fix") == 0)
 130         {
 131             meth = 7;
 132         }
 133         else
 134         {
 135             Scierror(999, _("%s: Wrong value for input argument #%d: It must be one of the following strings: adams, stiff, rk, rkf, fix, root or discrete.\n"), "ode", 1);
 136             return types::Function::Error;
 137         }
 138     }
 139
 140     // *** check number of output args according the methode. ***
 141     if (meth < 3)
 142     {
 143         if (_iRetCount != 1 && _iRetCount != 3)
 144         {
 145             Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "ode", 1, 3);
 146             return types::Function::Error;
 147         }
 148     }
 149     else if (meth == 3)
 150     {
 151         if (_iRetCount == 3 || _iRetCount > 4)
 152         {
 153             Scierror(78, _("%s: Wrong number of output argument(s): %d, %d or %d expected.\n"), "ode", 1, 2, 4);
 154             return types::Function::Error;
 155         }
 156     }
 157     else // meth > 3
 158     {
 159         if (_iRetCount != 1)
 160         {
 161             Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ode", 1);
 162             return types::Function::Error;
 163         }
 164     }
 165
 166     // *** check type of input args and get it. ***
 167     // y0
 168     if (in[iPos]->isDouble())
 169     {
 170         pDblY0 = in[iPos]->getAs<types::Double>();
 171         if (pDblY0->isComplex())
 172         {
 173             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ode", iPos + 1);
 174             return types::Function::Error;
 175         }
 176     }
 177     /*
 178         else if(in[iPos]->isPoly())
 179         {
 180             pPolyY0 = in[iPos]->getAs<types::Polynom>();
 181             if(pPolyY0->isComplex())
 182             {
 183                 Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "ode", iPos+1);
 184                 return types::Function::Error;
 185             }
 186
 187             if(pPolyY0->isScalar() == false)
 188             {
 189                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar polynom expected.\n"), "ode", iPos+1);
 190                 return types::Function::Error;
 191             }
 192
 193             double* dbl = (double*)MALLOC(pPolyY0->getCoef()->getSize() * sizeof(double));
 194             vTransposeRealMatrix(   pPolyY0->getCoef()->get(),
 195                                     pPolyY0->getCoef()->getRows(),
 196                                     pPolyY0->getCoef()->getCols(),
 197                                     dbl);
 198
 199             pDblY0 = new types::Double(pPolyY0->getCoef()->getCols(), pPolyY0->getCoef()->getRows());
 200             pDblY0->set(dbl);
 201         }
 202     */
 203     else
 204     {
 205         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
 206         return types::Function::Error;
 207     }
 208
 209     // t0
 210     iPos++;
 211     if (in[iPos]->isDouble() == false)
 212     {
 213         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
 214         return types::Function::Error;
 215     }
 216
 217     pDblT0 = in[iPos]->getAs<types::Double>();
 218
 219     if (pDblT0->isScalar() == false)
 220     {
 221         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
 222         return types::Function::Error;
 223     }
 224
 225     // t
 226     iPos++;
 227     if (in[iPos]->isDouble() == false)
 228     {
 229         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
 230         return types::Function::Error;
 231     }
 232
 233     pDblT = in[iPos]->getAs<types::Double>();
 234
 235     // get next inputs
 236     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"ode");
 237     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
 238     deFunctionsManager->setOdeYRows(pDblY0->getRows());
 239     deFunctionsManager->setOdeYCols(pDblY0->getCols());
 240
 241     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
 242     *YSize = pDblY0->getSize();
 243     pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
 244     C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
 245
 246     if (meth == 4)
 247     {
 248         if (in.size() != 5)
 249         {
 250             Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 5);
 251             DifferentialEquation::removeDifferentialEquationFunctions();
 252             FREE(pdYData);
 253             FREE(YSize);
 254             return types::Function::Error;
 255         }
 256
 257         if (in[4]->isCallable() == false && in[4]->isString() == false && in[4]->isList() == false)
 258         {
 259             Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", 5);
 260             DifferentialEquation::removeDifferentialEquationFunctions();
 261             FREE(pdYData);
 262             FREE(YSize);
 263             return types::Function::Error;
 264         }
 265     }
 266
 267     for (iPos++; iPos < (int)in.size(); iPos++)
 268     {
 269         if (in[iPos]->isDouble())
 270         {
 271             if (pDblRtol == NULL && bFuncF == false)
 272             {
 273                 pDblRtol = in[iPos]->getAs<types::Double>();
 274                 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
 275                 {
 276                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
 277                     DifferentialEquation::removeDifferentialEquationFunctions();
 278                     FREE(pdYData);
 279                     FREE(YSize);
 280                     return types::Function::Error;
 281                 }
 282             }
 283             else if (pDblAtol == NULL && bFuncF == false)
 284             {
 285                 pDblAtol = in[iPos]->getAs<types::Double>();
 286                 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
 287                 {
 288                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
 289                     DifferentialEquation::removeDifferentialEquationFunctions();
 290                     FREE(pdYData);
 291                     FREE(YSize);
 292                     return types::Function::Error;
 293                 }
 294             }
 295             else if (pDblNg == NULL && bFuncF == true && meth == 3)
 296             {
 297                 pDblNg = in[iPos]->getAs<types::Double>();
 298             }
 299             else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
 300             {
 301                 if ((int)in.size() == iPos + 2)
 302                 {
 303                     if (in[iPos + 1]->isDouble() == false)
 304                     {
 305                         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 2);
 306                         DifferentialEquation::removeDifferentialEquationFunctions();
 307                         FREE(pdYData);
 308                         FREE(YSize);
 309                         return types::Function::Error;
 310                     }
 311
 312                     pDblW = in[iPos]->getAs<types::Double>();
 313                     pDblIw = in[iPos + 1]->getAs<types::Double>();
 314                     iPos++;
 315                 }
 316                 else
 317                 {
 318                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", iPos + 2);
 319                     DifferentialEquation::removeDifferentialEquationFunctions();
 320                     FREE(pdYData);
 321                     FREE(YSize);
 322                     return types::Function::Error;
 323                 }
 324             }
 325             else
 326             {
 327                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", iPos + 1);
 328                 DifferentialEquation::removeDifferentialEquationFunctions();
 329                 FREE(pdYData);
 330                 FREE(YSize);
 331                 return types::Function::Error;
 332             }
 333         }
 334         else if (in[iPos]->isCallable())
 335         {
 336             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
 337             if (bFuncF == false)
 338             {
 339                 deFunctionsManager->setFFunction(pCall);
 340                 bFuncF = true;
 341             }
 342             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
 343             {
 344                 deFunctionsManager->setJacFunction(pCall);
 345                 bFuncJac = true;
 346             }
 347             else if (bFuncG == false && meth == 3)
 348             {
 349                 deFunctionsManager->setGFunction(pCall);
 350                 bFuncG = true;
 351             }
 352             else
 353             {
 354                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
 355                 DifferentialEquation::removeDifferentialEquationFunctions();
 356                 FREE(pdYData);
 357                 FREE(YSize);
 358                 return types::Function::Error;
 359             }
 360         }
 361         else if (in[iPos]->isString())
 362         {
 363             types::String* pStr = in[iPos]->getAs<types::String>();
 364             bool bOK = false;
 365
 366             if (bFuncF == false)
 367             {
 368                 bOK = deFunctionsManager->setFFunction(pStr);
 369                 bFuncF = true;
 370             }
 371             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
 372             {
 373                 bOK = deFunctionsManager->setJacFunction(pStr);
 374                 bFuncJac = true;
 375             }
 376             else if (bFuncG == false && meth == 3)
 377             {
 378                 bOK = deFunctionsManager->setGFunction(pStr);
 379                 bFuncG = true;
 380             }
 381             else
 382             {
 383                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
 384                 DifferentialEquation::removeDifferentialEquationFunctions();
 385                 FREE(pdYData);
 386                 FREE(YSize);
 387                 return types::Function::Error;
 388             }
 389
 390             if (bOK == false)
 391             {
 392                 char* pst = wide_string_to_UTF8(pStr->get(0));
 393                 Scierror(50, _("%s: Subroutine not found: %s\n"), "ode", pst);
 394                 FREE(pst);
 395                 DifferentialEquation::removeDifferentialEquationFunctions();
 396                 FREE(pdYData);
 397                 FREE(YSize);
 398                 return types::Function::Error;
 399             }
 400         }
 401         else if (in[iPos]->isList())
 402         {
 403             types::List* pList = in[iPos]->getAs<types::List>();
 404
 405             if (pList->getSize() == 0)
 406             {
 407                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, "(string empty)");
 408                 DifferentialEquation::removeDifferentialEquationFunctions();
 409                 FREE(pdYData);
 410                 FREE(YSize);
 411                 return types::Function::Error;
 412             }
 413
 414             if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
 415             {
 416                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
 417                 DifferentialEquation::removeDifferentialEquationFunctions();
 418                 FREE(pdYData);
 419                 FREE(YSize);
 420                 return types::Function::Error;
 421             }
 422
 423             if (pList->get(0)->isString())
 424             {
 425                 types::String* pStr = pList->get(0)->getAs<types::String>();
 426                 bool bOK = false;
 427
 428                 if (bFuncF == false)
 429                 {
 430                     bFuncF = true;
 431                     bOK = deFunctionsManager->setFFunction(pStr);
 432                     sizeOfpdYData = *YSize;
 433                 }
 434                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
 435                 {
 436                     bFuncJac = true;
 437                     bOK = deFunctionsManager->setJacFunction(pStr);
 438                     if (sizeOfpdYData == 0)
 439                     {
 440                         sizeOfpdYData = *YSize;
 441                     }
 442                 }
 443                 else if (bFuncG == false && meth == 3)
 444                 {
 445                     bFuncG = true;
 446                     bOK = deFunctionsManager->setGFunction(pStr);
 447                     if (sizeOfpdYData == 0)
 448                     {
 449                         sizeOfpdYData = *YSize;
 450                     }
 451                 }
 452
 453                 if (bOK == false)
 454                 {
 455                     char* pst = wide_string_to_UTF8(pStr->get(0));
 456                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, pst);
 457                     FREE(pst);
 458                     DifferentialEquation::removeDifferentialEquationFunctions();
 459                     FREE(pdYData);
 460                     FREE(YSize);
 461                     return types::Function::Error;
 462                 }
 463
 464                 int* sizeTemp = YSize;
 465                 int totalSize = sizeOfpdYData;
 466
 467                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
 468                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
 469
 470                 std::vector<types::Double*> vpDbl;
 471                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 472                 {
 473                     if (pList->get(iter + 1)->isDouble() == false)
 474                     {
 475                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "ode", iPos + 1, iter + 1);
 476                         DifferentialEquation::removeDifferentialEquationFunctions();
 477                         FREE(pdYData);
 478                         FREE(YSize);
 479                         return types::Function::Error;
 480                     }
 481
 482                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
 483                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
 484                     totalSize += YSize[sizeOfYSize + iter];
 485                 }
 486
 487                 double* pdYDataTemp = pdYData;
 488                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
 489                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
 490
 491                 int position = sizeOfpdYData;
 492                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 493                 {
 494                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
 495                     position += vpDbl[iter]->getSize();
 496                 }
 497                 vpDbl.clear();
 498                 sizeOfpdYData = totalSize;
 499                 sizeOfYSize += pList->getSize() - 1;
 500                 FREE(pdYDataTemp);
 501                 FREE(sizeTemp);
 502             }
 503             else if (pList->get(0)->isCallable())
 504             {
 505                 if (bFuncF == false)
 506                 {
 507                     bFuncF = true;
 508                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
 509                     for (int iter = 1; iter < pList->getSize(); iter++)
 510                     {
 511                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
 512                     }
 513                 }
 514                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
 515                 {
 516                     bFuncJac = true;
 517                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
 518                     for (int iter = 1; iter < pList->getSize(); iter++)
 519                     {
 520                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
 521                     }
 522                 }
 523                 else if (bFuncG == false && meth == 3)
 524                 {
 525                     bFuncG = true;
 526                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
 527                     for (int iter = 1; iter < pList->getSize(); iter++)
 528                     {
 529                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
 530                     }
 531                 }
 532             }
 533             else
 534             {
 535                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "ode", iPos + 1);
 536                 DifferentialEquation::removeDifferentialEquationFunctions();
 537                 FREE(pdYData);
 538                 FREE(YSize);
 539                 return types::Function::Error;
 540             }
 541         }
 542         else
 543         {
 544             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "ode", iPos + 1);
 545             DifferentialEquation::removeDifferentialEquationFunctions();
 546             FREE(pdYData);
 547             FREE(YSize);
 548             return types::Function::Error;
 549         }
 550     }
 551
 552     if (bFuncF == false)
 553     {
 554         int val = (meth == 3) ? 3 : 1;
 555         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + val);
 556         DifferentialEquation::removeDifferentialEquationFunctions();
 557         FREE(pdYData);
 558         FREE(YSize);
 559         return types::Function::Error;
 560     }
 561     if (pDblNg == NULL && meth == 3)
 562     {
 563         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 2);
 564         DifferentialEquation::removeDifferentialEquationFunctions();
 565         FREE(pdYData);
 566         FREE(YSize);
 567         return types::Function::Error;
 568     }
 569     if (bFuncG == false && meth == 3)
 570     {
 571         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 1);
 572         DifferentialEquation::removeDifferentialEquationFunctions();
 573         FREE(pdYData);
 574         FREE(YSize);
 575         return types::Function::Error;
 576     }
 577
 578     // *** Initialization. ***
 579
 580     int itol        = 1;
 581     int itask       = 1; // itask = 1 for normal computation of output values of y at t = tout.
 582     int istate      = 1; // istate = integer flag (input and output).  set istate = 1.
 583     int iopt        = 0; // iopt   = 0 to indicate no optional inputs used.
 584     int jt          = 2; // jacobian type indicator.  set jt = 2.
 585     int ml          = -1;
 586     int mu          = -1;
 587
 588     // work tab
 589     double* rwork   = NULL;
 590     int* iwork      = NULL;
 591     int rworkSize   = 0;
 592     int iworkSize   = 0;
 593
 594     // contain ls0001, lsa001 and eh0001 structures
 595     double* dStructTab  = NULL;
 596     int* iStructTab     = NULL;
 597     int dStructTabSize  = 0;
 598     int iStructTabSize  = 0;
 599
 600     int rwSize  = 0; // rwSize = dStructTab + rworkSize
 601     int iwSize  = 0; // iwSize = iStructTab + iworkSize
 602
 603     // structures used by lsoda and lsode
 604     double* ls0001d = &(C2F(ls0001).tret);
 605     int* ls0001i    = &(C2F(ls0001).illin);
 606     double* lsa001d = &(C2F(lsa001).tsw);
 607     int* lsa001i    = &(C2F(lsa001).insufr);
 608     double* lsr001d = &(C2F(lsr001).rownr3[0]);
 609     int* lsr001i    = &(C2F(lsr001).lg0);
 610     int* eh0001i    = &(C2F(eh0001).mesflg);
 611
 612     // get %ODEOPTIONS
 613     types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
 614     if (pIT != NULL && pIT->isDouble())
 615     {
 616         pDblOdeOptions = pIT->getAs<types::Double>();
 617         if (pDblOdeOptions->getSize() == 12)
 618         {
 619             iopt    = 1;
 620             itask   = (int)pDblOdeOptions->get(0);
 621             jt      = (int)pDblOdeOptions->get(5);
 622             ml      = (int)pDblOdeOptions->get(10);
 623             mu      = (int)pDblOdeOptions->get(11);
 624         }
 625     }
 626
 627     if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
 628     {
 629         Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "ode", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
 630         DifferentialEquation::removeDifferentialEquationFunctions();
 631         FREE(pdYData);
 632         FREE(YSize);
 633         return types::Function::Error;
 634     }
 635
 636     if (jt < 0 || jt > 5)
 637     {
 638         Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "ode", 0, 5);
 639         DifferentialEquation::removeDifferentialEquationFunctions();
 640         FREE(pdYData);
 641         FREE(YSize);
 642         return types::Function::Error;
 643     }
 644
 645     if (iopt == 0 && bFuncJac)
 646     {
 647         jt = 1;
 648     }
 649
 650     if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
 651     {
 652         sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "ode");
 653     }
 654
 655     if (bFuncJac == false && (jt == 1 || jt == 4))
 656     {
 657         if (getWarningMode())
 658         {
 659             sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "ode");
 660         }
 661
 662         jt = 2;
 663     }
 664
 665     //compute itol and set the tolerances rtol and atol.
 666     double* rtol = NULL;
 667     double* atol = NULL;
 668
 669     if (pDblRtol)
 670     {
 671         if (pDblRtol->isScalar())
 672         {
 673             rtol = (double*)MALLOC(sizeof(double));
 674             *rtol = pDblRtol->get(0);
 675         }
 676         else
 677         {
 678             rtol = pDblRtol->get();
 679             itol += 2;
 680         }
 681     }
 682     else
 683     {
 684         rtol = (double*)MALLOC(sizeof(double));
 685         if (meth == 6 || meth == 7)
 686         {
 687             *rtol = 1.e-3;
 688         }
 689         else
 690         {
 691             *rtol = 1.0e-7;
 692         }
 693     }
 694
 695     if (pDblAtol)
 696     {
 697         if (pDblAtol->isScalar())
 698         {
 699             atol = (double*)MALLOC(sizeof(double));
 700             *atol = pDblAtol->get(0);
 701         }
 702         else
 703         {
 704             atol = pDblAtol->get();
 705             itol ++;
 706         }
 707     }
 708     else
 709     {
 710         atol = (double*)MALLOC(sizeof(double));
 711         if (meth == 6 || meth == 7)
 712         {
 713             *atol = 1.e-4;
 714         }
 715         else
 716         {
 717             *atol = 1.0e-9;
 718         }
 719     }
 720
 721     // Compute rwork, iwork size.
 722     // Create them.
 723     int nyh = (*YSize);
 724     if (pDblW) // structure ls0001 have been restored.
 725     {
 726         nyh = C2F(ls0001).nyh;
 727     }
 728
 729     int mxordn = 12;
 730     int mxords = 5;
 731
 732     if (iopt == 1)
 733     {
 734         mxordn = (int)pDblOdeOptions->get(7);
 735         mxords = (int)pDblOdeOptions->get(8);
 736     }
 737
 738     if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
 739     {
 740         if (getWarningMode())
 741         {
 742             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);
 743         }
 744
 745         mxordn = 12;
 746         mxords = 5;
 747     }
 748
 749     int maxord  = mxords;
 750     int lmat    = 0;
 751     int lrn     = 0;
 752     int lrs     = 0;
 753
 754     switch (meth)
 755     {
 756         case 3: // lsodar (root)
 757         {
 758             //             lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
 759             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
 760             //          where
 761             //             nyh    = the initial value of neq,
 762             //             mxordn = 12, unless a smaller value is given as an
 763             //                      optional input,
 764             //             mxords = 5, unless a smaller value is given as an
 765             //                      optional input,
 766             //             lmat   = length of matrix work space..
 767             //             lmat   = neq**2 + 2              if jt = 1 or 2,
 768             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
 769
 770             lrn = 3 * (int)pDblNg->get(0);
 771             lrs = lrn;
 772             dStructTabSize = 246 - 241;
 773             iStructTabSize = 59 - 50;
 774             jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
 775         }
 776         case 0: // lsoda
 777         {
 778             //             lrn = 20 + nyh*(mxordn+1) + 3*neq,
 779             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
 780             //          where
 781             //             nyh    = the initial value of neq,
 782             //             mxordn = 12, unless a smaller value is given as an
 783             //                      optional input,
 784             //             mxords = 5, unless a smaller value is given as an
 785             //                      optional input,
 786             //             lmat   = length of matrix work space..
 787             //             lmat   = neq**2 + 2              if jt = 1 or 2,
 788             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
 789
 790             if (jt == 1 || jt == 2)
 791             {
 792                 lmat = (*YSize) * (*YSize) + 2;       //  if jt = 1 or 2,
 793             }
 794             else if (jt == 4 || jt == 5)
 795             {
 796                 //ml and mu = -1 in all cases
 797                 lmat = (2 * ml + mu + 1) * (*YSize) + 2; //  if jt = 4 or 5.
 798             }
 799
 800             lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
 801             lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
 802
 803             rworkSize   = max(lrn, lrs);
 804             iworkSize   = 20 + *YSize;
 805
 806             dStructTabSize += 241;
 807             iStructTabSize += 50;
 808
 809             break;
 810         }
 811         case 1: // lsode (adams) (non stiff)
 812         {
 813             maxord = mxordn;
 814         }
 815         case 2: // lsode (stiff)
 816         {
 817             //          20 + nyh*(maxord + 1) + 3*neq + lmat
 818             //        where
 819             //          nyh    = the initial value of neq,
 820             //          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
 821             //                   smaller value is given as an optional input),
 822             //          lmat   = 0             if miter = 0,
 823             //          lmat   = neq**2 + 2    if miter is 1 or 2,
 824             //          lmat   = neq + 2       if miter = 3, and
 825             //          lmat   = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
 826
 827             if (jt == 1 || jt == 2)
 828             {
 829                 lmat = (*YSize) * (*YSize) + 2;
 830             }
 831             else if (jt == 3)
 832             {
 833                 lmat = (*YSize) + 2;
 834             }
 835             else if (jt == 4 || jt == 5)
 836             {
 837                 lmat = (2 * ml + mu + 1) * (*YSize) + 2;
 838             }
 839
 840             rworkSize = 20 + nyh * (maxord + 1) + 3 * (*YSize) + lmat;
 841             iworkSize = 20; // if jt = 0 or 3
 842             if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq  otherwise
 843             {
 844                 iworkSize += (*YSize);
 845             }
 846
 847             dStructTabSize = 219;
 848             iStructTabSize = 41;
 849             break;
 850         }
 851         case 4: // lsdisc (discrete)
 852         {
 853             rworkSize = *YSize;
 854             iworkSize = 1;
 855             break;
 856         }
 857         case 5: // lsrgk (rk)
 858         {
 859             rworkSize = 9 * (*YSize);
 860             iworkSize = 1;
 861             break;
 862         }
 863         case 6: // rkf45 (rkf)
 864         {
 865             rworkSize = 3 + 8 * (*YSize);
 866             iworkSize = 5;
 867             break;
 868         }
 869         case 7: // rksimp (fix)
 870         {
 871             rworkSize = 3 + 8 * (*YSize);
 872             iworkSize = 1;
 873             break;
 874         }
 875     }
 876
 877     rwSize = rworkSize + dStructTabSize;
 878     iwSize = iworkSize + iStructTabSize;
 879
 880     rwork = (double*)MALLOC(rworkSize * sizeof(double));
 881     iwork = (int*)MALLOC(iworkSize * sizeof(int));
 882
 883     if (meth < 4)
 884     {
 885         if (pDblW && pDblIw)
 886         {
 887             if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
 888             {
 889                 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "ode", rwSize, iwSize);
 890                 DifferentialEquation::removeDifferentialEquationFunctions();
 891                 FREE(pdYData);
 892                 FREE(YSize);
 893                 FREE(rwork);
 894                 FREE(iwork);
 895                 if (jroot)
 896                 {
 897                     FREE(jroot);
 898                 }
 899                 if (itol == 1 || itol == 3)
 900                 {
 901                     FREE(atol);
 902                 }
 903                 if (itol < 3)
 904                 {
 905                     FREE(rtol);
 906                 }
 907                 return types::Function::Error;
 908             }
 909
 910             istate = 2; // 1 means this is the first call | 2  means this is not the first call
 911
 912             dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
 913             iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
 914
 915             C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
 916             C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
 917
 918             for (int i = 0; i < iworkSize; i++)
 919             {
 920                 iwork[i] = (int)pDblIw->get(i);
 921             }
 922
 923             for (int i = 0; i < iStructTabSize; i++)
 924             {
 925                 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
 926             }
 927         }
 928         else
 929         {
 930             //if iopt = 0 default value are used without set it.
 931             if (iopt == 1)
 932             {
 933                 if (itask >= 4)
 934                 {
 935                     rwork[0] = pDblOdeOptions->get(1); // tcrit
 936                 }
 937
 938                 rwork[4] = pDblOdeOptions->get(2);      // h0 =0
 939                 rwork[5] = pDblOdeOptions->get(3);      // hmax = %inf
 940                 rwork[6] = pDblOdeOptions->get(4);      // hmin = 0
 941
 942                 if (jt == 4 || jt == 5)
 943                 {
 944                     iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
 945                     iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
 946                 }
 947
 948                 if (meth == 0 || meth == 3)
 949                 {
 950                     // lsoda/lsodar
 951                     iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
 952                     iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
 953                     iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
 954                 }
 955                 else // 1 or 2
 956                 {
 957                     // lsode
 958                     if (meth == 1)
 959                     {
 960                         iwork[4] = (int)pDblOdeOptions->get(7); // mxordn = 12
 961                     }
 962                     else // meth == 2
 963                     {
 964                         iwork[4] = (int)pDblOdeOptions->get(8); // mxords = 5
 965                     }
 966                 }
 967
 968                 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
 969                 iwork[6] = 10;  // mxhnil = 10 maximum number of messages printed per problem
 970             }
 971         }
 972     }
 973
 974     if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
 975     {
 976         sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\tiopt = %d\n"), itask, meth, jt, ml, mu, iopt);
 977         sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
 978     }
 979
 980     if (meth < 4)
 981     {
 982         if (pDblW && pDblIw)
 983         {
 984             int dPos    = 0;
 985             int iPos    = 0;
 986             int dSize   = 219;
 987
 988             //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
 989             C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
 990             memcpy(ls0001i, iStructTab, 39 * sizeof(int));
 991
 992             dPos = dSize;
 993             iPos = 39;
 994
 995             //restore lsa001 from w (rwork) and iw (iwork).
 996             if (meth == 0 || meth == 3)
 997             {
 998                 dSize = 22;
 999                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
1000                 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
1001
1002                 dPos += dSize;
1003                 iPos += 9;
1004             }
1005
1006             //restore lsr001
1007             if (meth == 3)
1008             {
1009                 dSize = 5;
1010                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
1011                 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
1012             }
1013         }
1014     }
1015
1016     // *** Perform operation. ***
1017     int err = 0;
1018     double t0 = pDblT0->get(0);
1019     bool bOneStep = false;
1020
1021     if (itask == 5 || itask == 3 || itask == 2)
1022     {
1023         bOneStep = true;
1024         if (getWarningMode() && pDblT->isScalar() == false)
1025         {
1026             sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1027         }
1028     }
1029
1030     if (bOneStep)
1031     {
1032         std::list<double*> pDblYOutList = std::list<double*>();
1033         std::list<double> pDblTOutList = std::list<double>();
1034         int iLastT = pDblT->getSize() - 1;
1035         double t = pDblT->get(iLastT);
1036         int iDir = t - t0 < 0 ? -1 : 1;
1037
1038         do
1039         {
1040             std::string strMeth;
1041             try
1042             {
1043                 switch (meth)
1044                 {
1045                     case 1: // lsode (adams)
1046                     case 2: // lsode (stiff)
1047                     {
1048                         strMeth = "lsode";
1049                         int jacType = 10 * meth + jt;
1050                         C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1051                         break;
1052                     }
1053                     case 3: // lsodar
1054                     {
1055                         strMeth = "lsodar";
1056                         int ng = (int)pDblNg->get(0);
1057                         C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1058                         break;
1059                     }
1060                     case 4: // lsdisc (discrete)
1061                     {
1062                         strMeth = "lsdisc";
1063                         C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1064                         break;
1065                     }
1066                     case 5: // lsrgk (rk)
1067                     {
1068                         strMeth = "lsrgk";
1069                         C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1070                         break;
1071                     }
1072                     case 6: // rkf45 (rkf)
1073                     {
1074                         strMeth = "rkf45";
1075                         C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1076                         break;
1077                     }
1078                     case 7: // rksimp (fix)
1079                     {
1080                         strMeth = "rksimp";
1081                         C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1082                         break;
1083                     }
1084                     default: // case 0: // lsoda
1085                     {
1086                         strMeth = "lsoda";
1087                         C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1088                     }
1089                 }
1090
1091                 // check error
1092                 err = checkOdeError(meth, istate);
1093                 if (err == 1)
1094                 {
1095                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1096                 }
1097             }
1098             catch (ast::ScilabError &e)
1099             {
1100                 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1101                 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1102                 Scierror(999, pstrMsg);
1103                 FREE(pstrMsg);
1104                 err = 1;
1105             }
1106
1107             // FREE allocated data
1108             if (err == 1) // error case
1109             {
1110                 DifferentialEquation::removeDifferentialEquationFunctions();
1111                 FREE(pdYData);
1112                 FREE(YSize);
1113                 FREE(rwork);
1114                 FREE(iwork);
1115                 if (jroot)
1116                 {
1117                     FREE(jroot);
1118                 }
1119                 if (dStructTab)
1120                 {
1121                     FREE(dStructTab);
1122                 }
1123                 if (iStructTab)
1124                 {
1125                     FREE(iStructTab);
1126                 }
1127                 if (itol == 1 || itol == 3)
1128                 {
1129                     FREE(atol);
1130                 }
1131                 if (itol < 3)
1132                 {
1133                     FREE(rtol);
1134                 }
1135                 return types::Function::Error;
1136             }
1137
1138             pDblYOutList.push_back(pdYData);
1139             pDblTOutList.push_back(t0);
1140
1141             if (err == 2) // warning case
1142             {
1143                 if (getWarningMode())
1144                 {
1145                     sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1146                 }
1147                 break;
1148             }
1149
1150             if (meth == 3 && istate == 3)
1151             {
1152                 // istate == 3  means the integration was successful, and one or more
1153                 //              roots were found before satisfying the stop condition
1154                 //              specified by itask.  see jroot.
1155                 if (getWarningMode())
1156                 {
1157                     sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1158                     for (int k = 0; k < pDblNg->get(0); k++)
1159                     {
1160                         sciprint("\t%d", jroot[k]);
1161                     }
1162                     sciprint("\n");
1163                 }
1164
1165                 break;
1166             }
1167         }
1168         while ((t0 - t) * iDir < 0);
1169
1170         int iSizeList = (int)pDblYOutList.size();
1171         /*
1172                 if(pPolyY0) // pPolyY0 is scalar
1173                 {
1174                     // result format = [2 x iSizeList]
1175                     // first lime is the time of result stored in second line.
1176                     int size = 2 * iSizeList;
1177                     int* iRanks = (int*)MALLOC(size * sizeof(int));
1178                     int iMaxRank = pPolyY0->getMaxRank();
1179                     for(int i = 0; i < iSizeList; i++)
1180                     {
1181                         iRanks[i * 2] = 1; // time rank
1182                         iRanks[i * 2 + 1] = iMaxRank; // result rank
1183                     }
1184
1185                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
1186                     pDblYOut = new types::Double(pDblY0->getRows(), 1);
1187                     pDblTOut = new types::Double(1, 1);
1188
1189                     for(int i = 0; i < iSizeList; i++)
1190                     {
1191                         // pDblY0 contain pPolyY0 (ie: pPolyY0 = 2 + x => pDblY0 = [2 1])
1192                         for(int j = 0; j < pDblY0->getRows(); j++)
1193                         {
1194                             pDblYOut->set(j, pDblYOutList.front()[j]);
1195                         }
1196                         pDblTOut->set(0, pDblTOutList.front());
1197
1198                         pPolyYOut->setCoef(i * 2, pDblTOut);
1199                         pPolyYOut->setCoef(i * 2 + 1, pDblYOut);
1200
1201                         pDblYOutList.pop_front();
1202                         pDblTOutList.pop_front();
1203                     }
1204                 }
1205                 else
1206                 {
1207         */
1208         pDblYOut = new types::Double(pDblY0->getRows() + 1, iSizeList);
1209         for (int i = 0; i < iSizeList; i++)
1210         {
1211             pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1212             for (int j = 0; j < *YSize; j++)
1213             {
1214                 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1215             }
1216             pDblYOutList.pop_front();
1217             pDblTOutList.pop_front();
1218         }
1219         //        }
1220     }
1221     else
1222     {
1223         // Create result
1224         /*        if(pPolyY0)
1225                 {
1226                     int size = pDblT->getSize();
1227                     int* iRanks = (int*)MALLOC(size * sizeof(int));
1228                     int iMaxRank = pPolyY0->getMaxRank();
1229                     for(int i = 0; i < size; i++)
1230                     {
1231                         iRanks[i] = iMaxRank;
1232                     }
1233
1234                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
1235
1236                     pDblYOut = new types::Double(pDblY0->getRows(), 1);
1237                 }
1238                 else
1239                 {
1240         */
1241         pDblYOut = new types::Double(pDblY0->getRows(), pDblY0->getCols() * pDblT->getSize());
1242         //        }
1243         bool bBreak = false;
1244         for (int i = 0; i < pDblT->getSize(); i++)
1245         {
1246             double t = pDblT->get(i);
1247             std::string strMeth;
1248
1249             if (itask >= 4 && t > rwork[0]) // rwork[0] => tcrit
1250             {
1251                 t = rwork[0];
1252                 bBreak = true;
1253             }
1254
1255             try
1256             {
1257                 switch (meth)
1258                 {
1259                     case 1: // lsode (adams)
1260                     case 2: // lsode (stiff)
1261                     {
1262                         strMeth = "lsode";
1263                         int jacType = 10 * meth + jt;
1264                         C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1265                         break;
1266                     }
1267                     case 3: // lsodar
1268                     {
1269                         strMeth = "lsodar";
1270                         int ng = (int)pDblNg->get(0);
1271                         C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1272                         break;
1273                     }
1274                     case 4: // lsdisc (discrete)
1275                     {
1276                         strMeth = "lsdisc";
1277                         C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1278                         break;
1279                     }
1280                     case 5: // lsrgk (rk)
1281                     {
1282                         strMeth = "lsrgk";
1283                         C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1284                         break;
1285                     }
1286                     case 6: // rkf45 (rkf)
1287                     {
1288                         strMeth = "rkf45";
1289                         C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1290                         break;
1291                     }
1292                     case 7: // rksimp (fix)
1293                     {
1294                         strMeth = "rksimp";
1295                         C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1296                         break;
1297                     }
1298                     default: // case 0: // lsoda
1299                     {
1300                         strMeth = "lsoda";
1301                         C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1302                     }
1303                 }
1304
1305                 // check error
1306                 err = checkOdeError(meth, istate);
1307                 if (err == 1)
1308                 {
1309                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1310                 }
1311             }
1312             catch (ast::ScilabError &e)
1313             {
1314                 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1315                 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1316                 Scierror(999, pstrMsg);
1317                 FREE(pstrMsg);
1318                 err = 1;
1319             }
1320
1321             // FREE allocated data
1322             if (err == 1) // error case
1323             {
1324                 DifferentialEquation::removeDifferentialEquationFunctions();
1325                 FREE(pdYData);
1326                 FREE(YSize);
1327                 FREE(rwork);
1328                 FREE(iwork);
1329                 if (jroot)
1330                 {
1331                     FREE(jroot);
1332                 }
1333                 if (dStructTab)
1334                 {
1335                     FREE(dStructTab);
1336                 }
1337                 if (iStructTab)
1338                 {
1339                     FREE(iStructTab);
1340                 }
1341                 if (itol == 1 || itol == 3)
1342                 {
1343                     FREE(atol);
1344                 }
1345                 if (itol < 3)
1346                 {
1347                     FREE(rtol);
1348                 }
1349                 return types::Function::Error;
1350             }
1351
1352             /*            if(pPolyYOut)
1353                         {
1354                             for(int j = 0; j < pDblY0->getRows(); j++)
1355                             {
1356                                 pDblYOut->set(j, pdYData[j]);
1357                             }
1358                             pPolyYOut->setCoef(i, pDblYOut);
1359                         }
1360                         else
1361                         {
1362             */
1363             for (int j = 0; j < *YSize; j++)
1364             {
1365                 pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1366             }
1367             //            }
1368
1369             if (err == 2) // warning case
1370             {
1371                 if (getWarningMode())
1372                 {
1373                     sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1374                 }
1375
1376                 types::Double* pDblYOutTemp = pDblYOut;
1377                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1378                 for (int k = 0; k < pDblYOut->getSize(); k++)
1379                 {
1380                     pDblYOut->set(k, pDblYOutTemp->get(k));
1381                 }
1382                 delete pDblYOutTemp;
1383
1384                 break;
1385             }
1386
1387             if (meth == 3 && istate == 3)
1388             {
1389                 // istate == 3  means the integration was successful, and one or more
1390                 //              roots were found before satisfying the stop condition
1391                 //              specified by itask.  see jroot.
1392
1393                 if (getWarningMode())
1394                 {
1395                     sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1396                     for (int k = 0; k < pDblNg->get(0); k++)
1397                     {
1398                         sciprint("\t%d", jroot[k]);
1399                     }
1400                     sciprint("\n");
1401                 }
1402
1403                 types::Double* pDblYOutTemp = pDblYOut;
1404                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1405                 for (int k = 0; k < pDblYOut->getSize(); k++)
1406                 {
1407                     pDblYOut->set(k, pDblYOutTemp->get(k));
1408                 }
1409                 delete pDblYOutTemp;
1410                 break;
1411             }
1412
1413             if (itask >= 4 && bBreak)
1414             {
1415                 types::Double* pDblYOutTemp = pDblYOut;
1416                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1417                 for (int k = 0; k < pDblYOut->getSize(); k++)
1418                 {
1419                     pDblYOut->set(k, pDblYOutTemp->get(k));
1420                 }
1421                 delete pDblYOutTemp;
1422                 break;
1423             }
1424         }
1425     }
1426
1427     if (meth < 4)
1428     {
1429         if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1430         {
1431             int dPos    = 0;
1432             int iPos    = 0;
1433             int dSize   = 219;
1434
1435             if (dStructTab == NULL)
1436             {
1437                 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1438                 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1439             }
1440
1441             C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1442             memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1443
1444             dPos = dSize;
1445             iPos = 39;
1446
1447             //save lsa001
1448             if (meth == 0 || meth == 3)
1449             {
1450                 dSize = 22;
1451                 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1452                 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1453
1454                 dPos += dSize;
1455                 iPos += 9;
1456             }
1457
1458             //save lsr001
1459             if (meth == 3)
1460             {
1461                 dSize = 5;
1462                 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1463                 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1464
1465                 iPos += 9;
1466             }
1467
1468             memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1469         }
1470     }
1471
1472     // *** Return result in Scilab. ***
1473     /*
1474         if(pPolyYOut)
1475         {
1476             out.push_back(pPolyYOut); // y
1477         }
1478         else
1479         {
1480     */
1481     out.push_back(pDblYOut); // y
1482     //    }
1483
1484     if (meth == 3 && _iRetCount >= 2)
1485     {
1486         int sizeOfRd = 1;
1487         int k = 0;
1488
1489         for (int i = 0; i < pDblNg->get(0); i++)
1490         {
1491             if (jroot[i])
1492             {
1493                 sizeOfRd++;
1494             }
1495         }
1496
1497         types::Double* pDblRd = new types::Double(1, sizeOfRd);
1498         //rd: The first entry contains the stopping time.
1499         pDblRd->set(0, C2F(lsr001).tlast);
1500         for (int i = 0; i < pDblNg->get(0); i++)
1501         {
1502             if (jroot[i])
1503             {
1504                 k++;
1505                 pDblRd->set(k, (double)i + 1);
1506             }
1507         }
1508         out.push_back(pDblRd); // rd
1509     }
1510
1511     if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1512     {
1513         types::Double* pDblWOut = new types::Double(1, rwSize);
1514         C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1515         C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1516
1517         types::Double* pDblIwOut = new types::Double(1, iwSize);
1518         for (int i = 0; i < iworkSize; i++)
1519         {
1520             pDblIwOut->set(i, (double)iwork[i]);
1521         }
1522
1523         for (int i = 0; i < iStructTabSize; i++)
1524         {
1525             pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1526         }
1527
1528         out.push_back(pDblWOut);  // w
1529         out.push_back(pDblIwOut); // iw
1530     }
1531
1532     // *** FREE. ***
1533     if (itol == 1 || itol == 3) // atol is scalar
1534     {
1535         FREE(atol);
1536     }
1537
1538     if (itol < 3) // rtol is scalar
1539     {
1540         FREE(rtol);
1541     }
1542
1543     FREE(pdYData);
1544     FREE(YSize);
1545     FREE(rwork);
1546     FREE(iwork);
1547
1548     if (jroot)
1549     {
1550         FREE(jroot);
1551     }
1552
1553     if (dStructTab)
1554     {
1555         FREE(dStructTab);
1556         FREE(iStructTab);
1557     }
1558
1559     DifferentialEquation::removeDifferentialEquationFunctions();
1560
1561     return types::Function::OK;
1562 }
1563