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