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