scilab/modules/differential_equations/sci_gateway/cpp/sci_impl.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
  23 extern "C"
  24 {
  25 #include "sci_malloc.h"
  26 #include "localization.h"
  27 #include "Scierror.h"
  28 #include "scifunctions.h"
  29 #include "elem_common.h"
  30 #include "checkodeerror.h"
  31 #include "sciprint.h"
  32 #include "common_structure.h"
  33 }
  34 /*--------------------------------------------------------------------------*/
  35
  36 types::Function::ReturnValue sci_impl(types::typed_list &in, int _iRetCount, types::typed_list &out)
  37 {
  38     // Methode
  39     types::String* pStrType     = NULL;
  40     const wchar_t * wcsType     = L"lsoda";
  41     int meth                    = 2;// default methode is stiff
  42
  43     // y0
  44     types::Double* pDblY0       = NULL;
  45     double* pdYData             = NULL; // contain y0 following by all args data in list case.
  46     int sizeOfpdYData           = 0;
  47
  48     // Other input args
  49     types::Double* pDblYdot0    = NULL;
  50     types::Double* pDblT0       = NULL;
  51     types::Double* pDblT        = NULL;
  52     types::Double* pDblRtol     = NULL;
  53     types::Double* pDblAtol     = NULL;
  54     types::Double* pDblW        = NULL;
  55     types::Double* pDblIw       = NULL;
  56
  57     // Result
  58     types::Double* pDblYOut     = NULL;
  59
  60     // Indicate if the function is given.
  61     bool bFuncF     = false; // res
  62     bool bFuncJac   = false; // jac
  63     bool bFuncG     = false; // adda
  64
  65     int iPos        = 0; // Position in types::typed_list in
  66     int maxord      = 5; // maxord = 12 (if meth = 1) or 5 (if meth = 2)
  67
  68     int sizeOfYSize = 1;
  69     int* YSize      = NULL;    // YSize(1) = size of y0,
  70     // YSize(n) = size of Args(n) in list case.
  71
  72     C2F(eh0001).mesflg  = 1; // flag to control printing of error messages in lapack routine.
  73     // 1 means print, 0 means no printing.
  74     C2F(eh0001).lunit   = 6; // 6 = stdout
  75
  76     int one = 1; // use in dcopy
  77
  78     // *** check the minimal number of input args. ***
  79     if (in.size() < 6 || in.size() > 12)
  80     {
  81         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "impl", 6, 12);
  82         return types::Function::Error;
  83     }
  84
  85     // *** check number of output args ***
  86     if (_iRetCount > 3 || _iRetCount == 2)
  87     {
  88         Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "impl", 1, 3);
  89         return types::Function::Error;
  90     }
  91
  92     // *** Get the methode. ***
  93     if (in[0]->isString())
  94     {
  95         pStrType = in[0]->getAs<types::String>();
  96         wcsType = pStrType->get(0);
  97         iPos++;
  98     }
  99
 100     if (iPos)
 101     {
 102         if (wcscmp(wcsType, L"adams") == 0)
 103         {
 104             meth = 1;
 105             maxord = 12;
 106         }
 107         else if (wcscmp(wcsType, L"stiff") == 0)
 108         {
 109             meth = 2;
 110         }
 111         else
 112         {
 113             Scierror(999, _("%s: Wrong value for input argument #%d: It must be one of the following strings: adams or stiff.\n"), "impl", 1);
 114             return types::Function::Error;
 115         }
 116     }
 117
 118     // *** check type of input args and get it. ***
 119     // y0
 120     if (in[iPos]->isDouble() == false)
 121     {
 122         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 123         return types::Function::Error;
 124     }
 125
 126     pDblY0 = in[iPos]->getAs<types::Double>();
 127
 128     if (pDblY0->isComplex())
 129     {
 130         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "impl", iPos + 1);
 131         return types::Function::Error;
 132     }
 133
 134     if (pDblY0->getCols() != 1 && pDblY0->getRows() != 1)
 135     {
 136         Scierror(999, _("%s: Wrong type for input argument #%d: A real vector expected.\n"), "impl", iPos + 1);
 137         return types::Function::Error;
 138     }
 139
 140     // ydot0
 141     iPos++;
 142     if (in[iPos]->isDouble() == false)
 143     {
 144         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 145         return types::Function::Error;
 146     }
 147
 148     pDblYdot0 = in[iPos]->getAs<types::Double>();
 149
 150     if (pDblYdot0->isComplex())
 151     {
 152         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "impl", iPos + 1);
 153         return types::Function::Error;
 154     }
 155
 156     if (pDblYdot0->getCols() != 1 && pDblYdot0->getRows() != 1)
 157     {
 158         Scierror(999, _("%s: Wrong type for input argument #%d: A real vector expected.\n"), "impl", iPos + 1);
 159         return types::Function::Error;
 160     }
 161
 162     // t0
 163     iPos++;
 164     if (in[iPos]->isDouble() == false)
 165     {
 166         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "impl", iPos + 1);
 167         return types::Function::Error;
 168     }
 169
 170     pDblT0 = in[iPos]->getAs<types::Double>();
 171
 172     if (pDblT0->isScalar() == false)
 173     {
 174         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "impl", iPos + 1);
 175         return types::Function::Error;
 176     }
 177
 178     // t
 179     iPos++;
 180     if (in[iPos]->isDouble() == false)
 181     {
 182         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 183         return types::Function::Error;
 184     }
 185
 186     pDblT = in[iPos]->getAs<types::Double>();
 187
 188     // get next inputs
 189     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"impl");
 190     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
 191
 192     YSize = (int*)malloc(sizeOfYSize * sizeof(int));
 193     *YSize = pDblY0->getSize();
 194     pdYData = (double*)malloc(pDblY0->getSize() * sizeof(double));
 195     C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
 196
 197     for (iPos++; iPos < in.size(); iPos++)
 198     {
 199         if (in[iPos]->isDouble())
 200         {
 201             if (pDblAtol == NULL && bFuncF == false)
 202             {
 203                 pDblAtol = in[iPos]->getAs<types::Double>();
 204                 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
 205                 {
 206                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "impl", pStrType ? 2 : 1, iPos + 1);
 207                     DifferentialEquation::removeDifferentialEquationFunctions();
 208                     free(pdYData);
 209                     free(YSize);
 210                     return types::Function::Error;
 211                 }
 212             }
 213             else if (pDblRtol == NULL && bFuncF == false)
 214             {
 215                 pDblRtol = in[iPos]->getAs<types::Double>();
 216                 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
 217                 {
 218                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "impl", pStrType ? 2 : 1, iPos + 1);
 219                     DifferentialEquation::removeDifferentialEquationFunctions();
 220                     free(pdYData);
 221                     free(YSize);
 222                     return types::Function::Error;
 223                 }
 224             }
 225             else if (pDblW == NULL && bFuncG == true)
 226             {
 227                 if (in.size() == iPos + 2)
 228                 {
 229                     if (in[iPos + 1]->isDouble() == false)
 230                     {
 231                         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 2);
 232                         DifferentialEquation::removeDifferentialEquationFunctions();
 233                         free(pdYData);
 234                         free(YSize);
 235                         return types::Function::Error;
 236                     }
 237
 238                     pDblW = in[iPos]->getAs<types::Double>();
 239                     pDblIw = in[iPos + 1]->getAs<types::Double>();
 240                     iPos++;
 241                 }
 242                 else
 243                 {
 244                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "impl", iPos + 2);
 245                     DifferentialEquation::removeDifferentialEquationFunctions();
 246                     free(pdYData);
 247                     free(YSize);
 248                     return types::Function::Error;
 249                 }
 250             }
 251             else
 252             {
 253                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "impl", iPos + 1);
 254                 DifferentialEquation::removeDifferentialEquationFunctions();
 255                 free(pdYData);
 256                 free(YSize);
 257                 return types::Function::Error;
 258             }
 259         }
 260         else if (in[iPos]->isCallable())
 261         {
 262             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
 263             if (bFuncF == false)
 264             {
 265                 deFunctionsManager->setFFunction(pCall);
 266                 bFuncF = true;
 267             }
 268             else if (bFuncG == false)
 269             {
 270                 deFunctionsManager->setGFunction(pCall);
 271                 bFuncG = true;
 272             }
 273             else if (bFuncJac == false)
 274             {
 275                 deFunctionsManager->setJacFunction(pCall);
 276                 bFuncJac = true;
 277             }
 278             else
 279             {
 280                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 281                 DifferentialEquation::removeDifferentialEquationFunctions();
 282                 free(pdYData);
 283                 free(YSize);
 284                 return types::Function::Error;
 285             }
 286         }
 287         else if (in[iPos]->isString())
 288         {
 289             types::String* pStr = in[iPos]->getAs<types::String>();
 290             bool bOK = false;
 291
 292             if (bFuncF == false)
 293             {
 294                 bOK = deFunctionsManager->setFFunction(pStr);
 295                 bFuncF = true;
 296             }
 297             else if (bFuncG == false)
 298             {
 299                 bOK = deFunctionsManager->setGFunction(pStr);
 300                 bFuncG = true;
 301             }
 302             else if (bFuncJac == false)
 303             {
 304                 bOK = deFunctionsManager->setJacFunction(pStr);
 305                 bFuncJac = true;
 306             }
 307             else
 308             {
 309                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 310                 DifferentialEquation::removeDifferentialEquationFunctions();
 311                 free(pdYData);
 312                 free(YSize);
 313                 return types::Function::Error;
 314             }
 315
 316             if (bOK == false)
 317             {
 318                 char* pst = wide_string_to_UTF8(pStr->get(0));
 319                 Scierror(50, _("%s: Subroutine not found: %s\n"), "impl", pst);
 320                 FREE(pst);
 321                 DifferentialEquation::removeDifferentialEquationFunctions();
 322                 free(pdYData);
 323                 free(YSize);
 324                 return types::Function::Error;
 325             }
 326         }
 327         else if (in[iPos]->isList())
 328         {
 329             types::List* pList = in[iPos]->getAs<types::List>();
 330
 331             if (pList->getSize() == 0)
 332             {
 333                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "impl", iPos + 1, "(string empty)");
 334                 DifferentialEquation::removeDifferentialEquationFunctions();
 335                 free(pdYData);
 336                 free(YSize);
 337                 return types::Function::Error;
 338             }
 339
 340             if (bFuncF && bFuncG && bFuncJac)
 341             {
 342                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "impl", iPos + 1);
 343                 DifferentialEquation::removeDifferentialEquationFunctions();
 344                 free(pdYData);
 345                 free(YSize);
 346                 return types::Function::Error;
 347             }
 348
 349             if (pList->get(0)->isString())
 350             {
 351                 types::String* pStr = pList->get(0)->getAs<types::String>();
 352                 bool bOK = false;
 353
 354                 if (bFuncF == false)
 355                 {
 356                     bFuncF = true;
 357                     bOK = deFunctionsManager->setFFunction(pStr);
 358                     sizeOfpdYData = *YSize;
 359                 }
 360                 else if (bFuncG == false)
 361                 {
 362                     bFuncG = true;
 363                     bOK = deFunctionsManager->setGFunction(pStr);
 364                     if (sizeOfpdYData == 0)
 365                     {
 366                         sizeOfpdYData = *YSize;
 367                     }
 368                 }
 369                 else if (bFuncJac == false)
 370                 {
 371                     bFuncJac = true;
 372                     bOK = deFunctionsManager->setJacFunction(pStr);
 373                     if (sizeOfpdYData == 0)
 374                     {
 375                         sizeOfpdYData = *YSize;
 376                     }
 377                 }
 378
 379                 if (bOK == false)
 380                 {
 381                     char* pst = wide_string_to_UTF8(pStr->get(0));
 382                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "impl", iPos + 1, pst);
 383                     FREE(pst);
 384                     DifferentialEquation::removeDifferentialEquationFunctions();
 385                     free(pdYData);
 386                     free(YSize);
 387                     return types::Function::Error;
 388                 }
 389
 390                 int* sizeTemp = YSize;
 391                 int totalSize = sizeOfpdYData;
 392
 393                 YSize = (int*)malloc((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
 394                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
 395
 396                 std::vector<types::Double*> vpDbl;
 397                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 398                 {
 399                     if (pList->get(iter + 1)->isDouble() == false)
 400                     {
 401                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "impl", iPos + 1, iter + 1);
 402                         DifferentialEquation::removeDifferentialEquationFunctions();
 403                         free(pdYData);
 404                         free(YSize);
 405                         return types::Function::Error;
 406                     }
 407
 408                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
 409                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
 410                     totalSize += YSize[sizeOfYSize + iter];
 411                 }
 412
 413                 double* pdYDataTemp = pdYData;
 414                 pdYData = (double*)malloc(totalSize * sizeof(double));
 415                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
 416
 417                 int position = sizeOfpdYData;
 418                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 419                 {
 420                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
 421                     position += vpDbl[iter]->getSize();
 422                 }
 423                 vpDbl.clear();
 424                 sizeOfpdYData = totalSize;
 425                 sizeOfYSize += pList->getSize() - 1;
 426                 free(pdYDataTemp);
 427                 free(sizeTemp);
 428             }
 429             else if (pList->get(0)->isCallable())
 430             {
 431                 if (bFuncF == false)
 432                 {
 433                     bFuncF = true;
 434                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
 435                     for (int iter = 1; iter < pList->getSize(); iter++)
 436                     {
 437                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
 438                     }
 439                 }
 440                 else if (bFuncG == false)
 441                 {
 442                     bFuncG = true;
 443                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
 444                     for (int iter = 1; iter < pList->getSize(); iter++)
 445                     {
 446                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
 447                     }
 448                 }
 449                 else if (bFuncJac == false)
 450                 {
 451                     bFuncJac = true;
 452                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
 453                     for (int iter = 1; iter < pList->getSize(); iter++)
 454                     {
 455                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
 456                     }
 457                 }
 458             }
 459             else
 460             {
 461                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "impl", iPos + 1);
 462                 DifferentialEquation::removeDifferentialEquationFunctions();
 463                 free(pdYData);
 464                 free(YSize);
 465                 return types::Function::Error;
 466             }
 467         }
 468         else
 469         {
 470             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "impl", iPos + 1);
 471             DifferentialEquation::removeDifferentialEquationFunctions();
 472             free(pdYData);
 473             free(YSize);
 474             return types::Function::Error;
 475         }
 476     }
 477
 478     if (bFuncF == false)
 479     {
 480         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "impl", in.size() + 1);
 481         DifferentialEquation::removeDifferentialEquationFunctions();
 482         free(pdYData);
 483         free(YSize);
 484         return types::Function::Error;
 485     }
 486
 487     if (bFuncG == false)
 488     {
 489         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "impl", in.size() + 1);
 490         DifferentialEquation::removeDifferentialEquationFunctions();
 491         free(pdYData);
 492         free(YSize);
 493         return types::Function::Error;
 494     }
 495
 496     // *** Initialization. ***
 497     double t0   = pDblT0->get(0);
 498     int itol    = 1;
 499     int iopt    = 0;
 500     int istate  = 1;
 501     int itask   = 1;
 502     int jt = bFuncJac ? 1 : 2;
 503     int jacType = 10 * meth + jt;
 504
 505     pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize());
 506
 507     // work tab
 508     double* rwork   = NULL;
 509     int* iwork      = NULL;
 510     int rworkSize   = 0;
 511     int iworkSize   = 0;
 512
 513     // contain ls0001, lsa001 and eh0001 structures
 514     double* dStructTab  = NULL;
 515     int* iStructTab     = NULL;
 516     int dStructTabSize  = 219;  // number of double in ls0001
 517     int iStructTabSize  = 41;   // number of int in ls0001 (39) + eh0001 (2)
 518
 519     int rwSize  = 0; // rwSize = dStructTab + rworkSize
 520     int iwSize  = 0; // iwSize = iStructTab + iworkSize
 521
 522     // structures used by lsoda and lsode
 523     double* ls0001d = &(C2F(ls0001).tret);
 524     int* ls0001i    = &(C2F(ls0001).illin);
 525     int* eh0001i    = &(C2F(eh0001).mesflg);
 526
 527     //compute itol and set the tolerances rtol and atol.
 528     double* rtol = NULL;
 529     double* atol = NULL;
 530
 531     if (pDblRtol)
 532     {
 533         if (pDblRtol->isScalar())
 534         {
 535             rtol = (double*)malloc(sizeof(double));
 536             *rtol = pDblRtol->get(0);
 537         }
 538         else
 539         {
 540             rtol = pDblRtol->get();
 541             itol += 2;
 542         }
 543     }
 544     else
 545     {
 546         rtol = (double*)malloc(sizeof(double));
 547         *rtol = 1.e-9;
 548     }
 549
 550     if (pDblAtol)
 551     {
 552         if (pDblAtol->isScalar())
 553         {
 554             atol = (double*)malloc(sizeof(double));
 555             *atol = pDblAtol->get(0);
 556         }
 557         else
 558         {
 559             atol = pDblAtol->get();
 560             itol ++;
 561         }
 562     }
 563     else
 564     {
 565         atol = (double*)malloc(sizeof(double));
 566         *atol = 1.e-7;
 567     }
 568
 569     // Compute rwork, iwork size.
 570     // Create them.
 571
 572     int nyh = (*YSize);
 573     if (pDblW) // structure ls0001 have been restored.
 574     {
 575         nyh = C2F(ls0001).nyh;
 576     }
 577
 578     rworkSize = 20 + nyh * (maxord + 1) + 3 * *YSize + *YSize * *YSize + 2;
 579     iworkSize = 20 + *YSize;
 580
 581     rwSize = rworkSize + dStructTabSize;
 582     iwSize = iworkSize + iStructTabSize;
 583
 584     rwork = (double*)malloc(rworkSize * sizeof(double));
 585     iwork = (int*)malloc(iworkSize * sizeof(int));
 586
 587     if (pDblW && pDblIw)
 588     {
 589         if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
 590         {
 591             Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "impl", rwSize, iwSize);
 592             DifferentialEquation::removeDifferentialEquationFunctions();
 593             free(pdYData);
 594             free(YSize);
 595             free(rwork);
 596             free(iwork);
 597             if (itol == 1 || itol == 3)
 598             {
 599                 free(atol);
 600             }
 601             if (itol < 3)
 602             {
 603                 free(rtol);
 604             }
 605             return types::Function::Error;
 606         }
 607
 608         istate = 2; // 1 means this is the first call | 2  means this is not the first call
 609
 610         // restore rwork from pDblW
 611         C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
 612
 613         // restore iwork from pDblIw
 614         iStructTab = (int*)malloc(iStructTabSize * sizeof(int));
 615         for (int i = 0; i < iworkSize; i++)
 616         {
 617             iwork[i] = (int)pDblIw->get(i);
 618         }
 619
 620         //restore ls0001d from pDblW
 621         C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, ls0001d, &one);
 622
 623         //restore ls0001i from pDblIw
 624         for (int i = 0; i < iStructTabSize; i++)
 625         {
 626             iStructTab[i] = (int)pDblIw->get(i + iworkSize);
 627         }
 628         memcpy(ls0001i, iStructTab, 39 * sizeof(int));
 629     }
 630
 631     // *** Perform operation. ***
 632     int err = 0;
 633     for (int i = 0; i < pDblT->getSize(); i++)
 634     {
 635         double t = pDblT->get(i);
 636         try
 637         {
 638             C2F(lsodi)(impl_f, impl_g, impl_jac, YSize, pdYData, pDblYdot0->get(), &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, &jacType);
 639
 640             // check error
 641             if (istate == 3)
 642             {
 643                 sciprint(_("The user-supplied subroutine res signalled lsodi to halt the integration and return (ires=2). Execution of the external function has failed.\n"));
 644                 err = 1;
 645                 Scierror(999, _("%s: %s exit with state %d.\n"), "impl", "lsodi", istate);
 646             }
 647             else
 648             {
 649                 err = checkOdeError(meth, istate);
 650                 if (err == 1)
 651                 {
 652                     Scierror(999, _("%s: %s exit with state %d.\n"), "impl", "lsodi", istate);
 653                 }
 654             }
 655         }
 656         catch (ast::ScilabError &e)
 657         {
 658             char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
 659             sciprint(_("%s: exception caught in '%s' subroutine.\n"), "impl", "lsodi");
 660             Scierror(999, pstrMsg);
 661             FREE(pstrMsg);
 662             err = 1;
 663         }
 664
 665         if (err == 1)
 666         {
 667             DifferentialEquation::removeDifferentialEquationFunctions();
 668             free(pdYData);
 669             free(YSize);
 670             free(rwork);
 671             free(iwork);
 672             if (iStructTab)
 673             {
 674                 free(iStructTab);
 675             }
 676             if (itol == 1 || itol == 3)
 677             {
 678                 free(atol);
 679             }
 680             if (itol < 3)
 681             {
 682                 free(rtol);
 683             }
 684             return types::Function::Error;
 685         }
 686
 687         for (int j = 0; j < *YSize; j++)
 688         {
 689             pDblYOut->set(i * (*YSize) + j, pdYData[j]);
 690         }
 691     }
 692
 693     if (_iRetCount > 2) //save ls0001 and eh0001 following pDblW and pDblIw.
 694     {
 695         int dSize   = 219;
 696
 697         if (iStructTab == NULL)
 698         {
 699             iStructTab = (int*)malloc(iStructTabSize * sizeof(int));
 700         }
 701
 702         if (dStructTab == NULL)
 703         {
 704             dStructTab = (double*)malloc(dStructTabSize * sizeof(double));
 705         }
 706
 707         // save ls0001
 708         C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
 709         memcpy(iStructTab, ls0001i, 39 * sizeof(int));
 710
 711         // save eh0001
 712         memcpy(&iStructTab[39], eh0001i, 2 * sizeof(int));
 713     }
 714
 715     // *** Return result in Scilab. ***
 716     out.push_back(pDblYOut);
 717
 718     if (_iRetCount > 2)
 719     {
 720         types::Double* pDblWOut = new types::Double(1, rwSize);
 721         C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
 722         C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
 723
 724         types::Double* pDblIwOut = new types::Double(1, iwSize);
 725         for (int i = 0; i < iworkSize; i++)
 726         {
 727             pDblIwOut->set(i, (double)iwork[i]);
 728         }
 729
 730         for (int i = 0; i < iStructTabSize; i++)
 731         {
 732             pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
 733         }
 734
 735         out.push_back(pDblWOut);
 736         out.push_back(pDblIwOut);
 737     }
 738
 739     // *** free. ***
 740     if (itol == 1 || itol == 3) // atol is scalar
 741     {
 742         free(atol);
 743     }
 744
 745     if (itol < 3) // rtol is scalar
 746     {
 747         free(rtol);
 748     }
 749
 750     free(pdYData);
 751     free(YSize);
 752     free(rwork);
 753     free(iwork);
 754
 755     if (dStructTab)
 756     {
 757         free(dStructTab);
 758     }
 759
 760     if (iStructTab)
 761     {
 762         free(iStructTab);
 763     }
 764
 765     DifferentialEquation::removeDifferentialEquationFunctions();
 766
 767     return types::Function::OK;
 768 }
 769