scilab/modules/differential_equations/sci_gateway/cpp/sci_dassl.cpp

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
   4 *
   5 * This file must be used under the terms of the CeCILL.
   6 * This source file is licensed as described in the file COPYING, which
   7 * you should have received as part of this distribution.  The terms
   8 * are also available at
   9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  10 *
  11 */
  12 /*--------------------------------------------------------------------------*/
  13
  14 #include "differential_equations_gw.hxx"
  15 #include "function.hxx"
  16 #include "double.hxx"
  17 #include "string.hxx"
  18 #include "list.hxx"
  19 #include "callable.hxx"
  20 #include "differentialequationfunctions.hxx"
  21 #include "runvisitor.hxx"
  22 #include "checkodeerror.hxx"
  23
  24 extern "C"
  25 {
  26 #include "sci_malloc.h"
  27 #include "localization.h"
  28 #include "Scierror.h"
  29 #include "sciprint.h"
  30 #include "scifunctions.h"
  31 #include "elem_common.h"
  32 #include "xerhlt.h"
  33 }
  34
  35 /*--------------------------------------------------------------------------*/
  36 types::Function::ReturnValue sci_dassl(types::typed_list &in, int _iRetCount, types::typed_list &out)
  37 {
  38     // input args
  39     types::Double* pDblX0   = NULL;
  40     types::Double* pDblT0   = NULL;
  41     types::Double* pDblT    = NULL;
  42     types::Double* pDblRtol = NULL;
  43     types::Double* pDblAtol = NULL;
  44     types::Double* pDblHd   = NULL;
  45
  46     // x0 = [y0, ydot0]
  47     double* pdYData         = NULL; // contain y0 following by all args data in list case.
  48     double* pdYdotData      = NULL;
  49     int sizeOfpdYData       = 0;
  50
  51     int iPos    = 0;
  52     int one     = 1;
  53
  54     int info[15]    = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  55
  56     double tstop    = 0;
  57     double maxstep  = 0;
  58     double stepin   = 0;
  59     int mu          = 0;
  60     int ml          = 0;
  61
  62     int sizeOfYSize = 1;
  63     int* YSize      = NULL;    // YSize(1) = size of y0,
  64     // YSize(n) = size of Args(n) in list case.
  65
  66     // Indicate if the function is given.
  67     bool bFuncF     = false;
  68     bool bFuncJac   = false;
  69
  70     // Indicate if info list is given.
  71     bool bListInfo  = false;
  72
  73     // *** check the minimal number of input args. ***
  74     if (in.size() < 4 || in.size() > 9)
  75     {
  76         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "dassl", 4, 9);
  77         return types::Function::Error;
  78     }
  79
  80     // *** check number of output args ***
  81     if (_iRetCount > 2)
  82     {
  83         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "dassl", 1, 2);
  84         return types::Function::Error;
  85     }
  86
  87     // *** check type of input args and get it. ***
  88     // x0 = [y0, yd0]
  89     if (in[iPos]->isDouble() == false)
  90     {
  91         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
  92         return types::Function::Error;
  93     }
  94
  95     pDblX0 = in[iPos]->getAs<types::Double>();
  96
  97     if (pDblX0->isComplex())
  98     {
  99         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "dassl", iPos + 1);
 100         return types::Function::Error;
 101     }
 102
 103     if (pDblX0->getCols() > 2)
 104     {
 105         Scierror(999, _("%s: Wrong size for input argument #%d: A real matrix with %d to %d colomn(s) expected.\n"), "dassl", iPos + 1, 1, 2);
 106         return types::Function::Error;
 107     }
 108
 109     if (pDblX0->getCols() == 1)
 110     {
 111         info[10] = 1;
 112     }
 113
 114     // t0
 115     iPos++;
 116     if (in[iPos]->isDouble() == false)
 117     {
 118         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "dassl", iPos + 1);
 119         return types::Function::Error;
 120     }
 121
 122     pDblT0 = in[iPos]->getAs<types::Double>();
 123
 124     if (pDblT0->isScalar() == false)
 125     {
 126         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "dassl", iPos + 1);
 127         return types::Function::Error;
 128     }
 129
 130     // t
 131     iPos++;
 132     if (in[iPos]->isDouble() == false)
 133     {
 134         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
 135         return types::Function::Error;
 136     }
 137
 138     pDblT = in[iPos]->getAs<types::Double>();
 139
 140     if (pDblT->isComplex())
 141     {
 142         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "dassl", iPos + 1);
 143         return types::Function::Error;
 144     }
 145
 146     // get next inputs
 147     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"dassl");
 148     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
 149
 150     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
 151     *YSize = pDblX0->getRows();
 152
 153     pdYData = (double*)MALLOC(*YSize * sizeof(double));
 154     pdYdotData = (double*)MALLOC(*YSize * sizeof(double));
 155
 156     C2F(dcopy)(YSize, pDblX0->get(), &one, pdYData, &one);
 157     if (pDblX0->getCols() == 2)
 158     {
 159         C2F(dcopy)(YSize, pDblX0->get() + *YSize, &one, pdYdotData, &one);
 160     }
 161     else
 162     {
 163         memset(pdYdotData, 0x00, *YSize);
 164     }
 165
 166     deFunctionsManager->setOdeYRows(pDblX0->getRows());
 167
 168     for (iPos++; iPos < in.size(); iPos++)
 169     {
 170         if (in[iPos]->isDouble())
 171         {
 172             if (pDblAtol == NULL && bFuncF == false)
 173             {
 174                 pDblAtol = in[iPos]->getAs<types::Double>();
 175                 if (pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
 176                 {
 177                     Scierror(267, _("%s: Wrong size for input argument #%d: A scalar or a matrix of size %d expected.\n"), "dassl", iPos + 1, pDblX0->getRows());
 178                     DifferentialEquation::removeDifferentialEquationFunctions();
 179                     FREE(pdYdotData);
 180                     FREE(pdYData);
 181                     FREE(YSize);
 182                     return types::Function::Error;
 183                 }
 184             }
 185             else if (pDblRtol == NULL && bFuncF == false)
 186             {
 187                 pDblRtol = in[iPos]->getAs<types::Double>();
 188                 if (pDblAtol->getSize() != pDblRtol->getSize())
 189                 {
 190                     Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "dassl", iPos + 1, pDblX0->getRows());
 191                     DifferentialEquation::removeDifferentialEquationFunctions();
 192                     FREE(pdYdotData);
 193                     FREE(pdYData);
 194                     FREE(YSize);
 195                     return types::Function::Error;
 196                 }
 197             }
 198             else if (pDblHd == NULL && bFuncF == true)
 199             {
 200                 pDblHd = in[iPos]->getAs<types::Double>();
 201                 if (in.size() != iPos + 1)
 202                 {
 203                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", iPos + 1);
 204                     DifferentialEquation::removeDifferentialEquationFunctions();
 205                     FREE(pdYdotData);
 206                     FREE(pdYData);
 207                     FREE(YSize);
 208                     return types::Function::Error;
 209                 }
 210             }
 211             else
 212             {
 213                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "dassl", iPos + 1);
 214                 DifferentialEquation::removeDifferentialEquationFunctions();
 215                 FREE(pdYdotData);
 216                 FREE(pdYData);
 217                 FREE(YSize);
 218                 return types::Function::Error;
 219             }
 220         }
 221         else if (in[iPos]->isCallable())
 222         {
 223             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
 224             if (bFuncF == false)
 225             {
 226                 deFunctionsManager->setFFunction(pCall);
 227                 bFuncF = true;
 228             }
 229             else if (bFuncJac == false)
 230             {
 231                 deFunctionsManager->setJacFunction(pCall);
 232                 bFuncJac = true;
 233             }
 234             else
 235             {
 236                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
 237                 DifferentialEquation::removeDifferentialEquationFunctions();
 238                 FREE(pdYdotData);
 239                 FREE(pdYData);
 240                 FREE(YSize);
 241                 return types::Function::Error;
 242             }
 243         }
 244         else if (in[iPos]->isString())
 245         {
 246             types::String* pStr = in[iPos]->getAs<types::String>();
 247             bool bOK = false;
 248
 249             if (bFuncF == false)
 250             {
 251                 bOK = deFunctionsManager->setFFunction(pStr);
 252                 bFuncF = true;
 253             }
 254             else if (bFuncJac == false)
 255             {
 256                 bOK = deFunctionsManager->setJacFunction(pStr);
 257                 bFuncJac = true;
 258             }
 259             else
 260             {
 261                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
 262                 DifferentialEquation::removeDifferentialEquationFunctions();
 263                 FREE(pdYdotData);
 264                 FREE(pdYData);
 265                 FREE(YSize);
 266                 return types::Function::Error;
 267             }
 268
 269             if (bOK == false)
 270             {
 271                 char* pst = wide_string_to_UTF8(pStr->get(0));
 272                 Scierror(50, _("%s: Subroutine not found: %s\n"), "dassl", pst);
 273                 FREE(pst);
 274                 DifferentialEquation::removeDifferentialEquationFunctions();
 275                 FREE(pdYdotData);
 276                 FREE(pdYData);
 277                 FREE(YSize);
 278                 return types::Function::Error;
 279             }
 280         }
 281         else if (in[iPos]->isList())
 282         {
 283             types::List* pList = in[iPos]->getAs<types::List>();
 284
 285             if (pList->getSize() == 0)
 286             {
 287                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, "(string empty)");
 288                 DifferentialEquation::removeDifferentialEquationFunctions();
 289                 FREE(pdYdotData);
 290                 FREE(pdYData);
 291                 FREE(YSize);
 292                 return types::Function::Error;
 293             }
 294
 295             if (bFuncF && bListInfo)
 296             {
 297                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
 298                 DifferentialEquation::removeDifferentialEquationFunctions();
 299                 FREE(pdYdotData);
 300                 FREE(pdYData);
 301                 FREE(YSize);
 302                 return types::Function::Error;
 303             }
 304
 305             if (pList->get(0)->isString())
 306             {
 307                 types::String* pStr = pList->get(0)->getAs<types::String>();
 308                 bool bOK = false;
 309
 310                 if (bFuncF == false)
 311                 {
 312                     bFuncF = true;
 313                     bOK = deFunctionsManager->setFFunction(pStr);
 314                     sizeOfpdYData = *YSize;
 315                 }
 316                 else if (bFuncJac == false)
 317                 {
 318                     bFuncJac = true;
 319                     bOK = deFunctionsManager->setJacFunction(pStr);
 320                     if (sizeOfpdYData == 0)
 321                     {
 322                         sizeOfpdYData = *YSize;
 323                     }
 324                 }
 325
 326                 if (bOK == false)
 327                 {
 328                     char* pst = wide_string_to_UTF8(pStr->get(0));
 329                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, pst);
 330                     FREE(pst);
 331                     DifferentialEquation::removeDifferentialEquationFunctions();
 332                     FREE(pdYdotData);
 333                     FREE(pdYData);
 334                     FREE(YSize);
 335                     return types::Function::Error;
 336                 }
 337
 338                 int* sizeTemp = YSize;
 339                 int totalSize = sizeOfpdYData;
 340
 341                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
 342                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
 343
 344                 std::vector<types::Double*> vpDbl;
 345                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 346                 {
 347                     if (pList->get(iter + 1)->isDouble() == false)
 348                     {
 349                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "dassl", iPos + 1, iter + 1);
 350                         DifferentialEquation::removeDifferentialEquationFunctions();
 351                         FREE(pdYdotData);
 352                         FREE(pdYData);
 353                         FREE(YSize);
 354                         return types::Function::Error;
 355                     }
 356
 357                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
 358                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
 359                     totalSize += YSize[sizeOfYSize + iter];
 360                 }
 361
 362                 double* pdYDataTemp = pdYData;
 363                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
 364                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
 365
 366                 int position = sizeOfpdYData;
 367                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 368                 {
 369                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
 370                     position += vpDbl[iter]->getSize();
 371                 }
 372                 vpDbl.clear();
 373                 sizeOfpdYData = totalSize;
 374                 sizeOfYSize += pList->getSize() - 1;
 375                 FREE(pdYDataTemp);
 376                 FREE(sizeTemp);
 377             }
 378             else if (pList->get(0)->isCallable())
 379             {
 380                 if (bFuncF == false)
 381                 {
 382                     bFuncF = true;
 383                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
 384                     for (int iter = 1; iter < pList->getSize(); iter++)
 385                     {
 386                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
 387                     }
 388                 }
 389                 else if (bFuncJac == false)
 390                 {
 391                     bFuncJac = true;
 392                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
 393                     for (int iter = 1; iter < pList->getSize(); iter++)
 394                     {
 395                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
 396                     }
 397                 }
 398             }
 399             else if (pList->get(0)->isDouble() && bFuncF == true)
 400             {
 401                 if (pList->getSize() != 7)
 402                 {
 403                     Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "dassl", iPos + 1, 7);
 404                     DifferentialEquation::removeDifferentialEquationFunctions();
 405                     FREE(pdYdotData);
 406                     FREE(pdYData);
 407                     FREE(YSize);
 408                     return types::Function::Error;
 409                 }
 410
 411                 for (int i = 0; i < 7; i++) // info = list([],0,[],[],[],0,0)
 412                 {
 413                     if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 6)))
 414                     {
 415                         if (i == 1 || i == 5 || i == 6)
 416                         {
 417                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "dassl", iPos + 1, i);
 418                         }
 419                         else
 420                         {
 421                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "dassl", iPos + 1, i);
 422                         }
 423                         DifferentialEquation::removeDifferentialEquationFunctions();
 424                         FREE(pdYdotData);
 425                         FREE(pdYData);
 426                         FREE(YSize);
 427                         return types::Function::Error;
 428                     }
 429                 }
 430
 431                 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
 432                 if (pDblTemp->getSize() != 0)
 433                 {
 434                     info[3] = 1;
 435                     tstop = pDblTemp->get(0);
 436                 }
 437
 438                 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
 439
 440                 pDblTemp = pList->get(2)->getAs<types::Double>();
 441                 if (pDblTemp->getSize() == 2)
 442                 {
 443                     info[5] = 1;
 444                     ml = (int)pDblTemp->get(0);
 445                     mu = (int)pDblTemp->get(1);
 446                     deFunctionsManager->setMl(ml);
 447                     deFunctionsManager->setMu(mu);
 448                 }
 449                 else if (pDblTemp->getSize() != 0)
 450                 {
 451                     Scierror(267, _("%s: Wrong size for input argument #%d: Argument %d in te list must be of size %d.\n"), "dassl", iPos + 1, 3, 2);
 452                     DifferentialEquation::removeDifferentialEquationFunctions();
 453                     FREE(pdYdotData);
 454                     FREE(pdYData);
 455                     FREE(YSize);
 456                     return types::Function::Error;
 457                 }
 458
 459                 pDblTemp = pList->get(3)->getAs<types::Double>();
 460                 if (pDblTemp->getSize() != 0)
 461                 {
 462                     info[6] = 1;
 463                     maxstep = pDblTemp->get(0);
 464                 }
 465
 466                 pDblTemp = pList->get(4)->getAs<types::Double>();
 467                 if (pDblTemp->getSize() != 0)
 468                 {
 469                     info[7] = 1;
 470                     stepin = pDblTemp->get(0);
 471                 }
 472
 473                 info[9]  = (int)pList->get(5)->getAs<types::Double>()->get(0);
 474                 if (pList->get(6)->getAs<types::Double>()->get(0) == 1)
 475                 {
 476                     info[10] = 1;
 477                 }
 478
 479                 bListInfo = true;
 480             }
 481             else
 482             {
 483                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string, a function or a matrix in case of argument info.\n"), "dassl", iPos + 1);
 484                 DifferentialEquation::removeDifferentialEquationFunctions();
 485                 FREE(pdYdotData);
 486                 FREE(pdYData);
 487                 FREE(YSize);
 488                 return types::Function::Error;
 489             }
 490         }
 491         else
 492         {
 493             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "dassl", iPos + 1);
 494             DifferentialEquation::removeDifferentialEquationFunctions();
 495             FREE(pdYdotData);
 496             FREE(pdYData);
 497             FREE(YSize);
 498             return types::Function::Error;
 499         }
 500     }
 501
 502     if (bFuncF == false)
 503     {
 504         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", in.size() + 1);
 505         DifferentialEquation::removeDifferentialEquationFunctions();
 506         FREE(pdYdotData);
 507         FREE(pdYData);
 508         FREE(YSize);
 509         return types::Function::Error;
 510     }
 511
 512     if (bFuncJac == true)
 513     {
 514         info[4] = 1;
 515     }
 516
 517     // *** Initialization. ***
 518     double t0   = pDblT0->get(0);
 519     double rpar = 0;
 520     int ipar    = 0;
 521     int idid    = 0;
 522     int maxord  = 5;
 523
 524     //compute itol and set the tolerances rtol and atol.
 525     double* rtol = NULL;
 526     double* atol = NULL;
 527
 528     if (pDblAtol)
 529     {
 530         if (pDblAtol->isScalar())
 531         {
 532             atol  = (double*)MALLOC(sizeof(double));
 533             *atol = pDblAtol->get(0);
 534         }
 535         else
 536         {
 537             atol    = pDblAtol->get();
 538             info[1] = 1;
 539         }
 540     }
 541     else
 542     {
 543         atol  = (double*)MALLOC(sizeof(double));
 544         *atol = 1.e-7;
 545     }
 546
 547     if (pDblRtol)
 548     {
 549         if (pDblRtol->isScalar())
 550         {
 551             rtol  = (double*)MALLOC(sizeof(double));
 552             *rtol = pDblRtol->get(0);
 553         }
 554         else
 555         {
 556             rtol = pDblRtol->get();
 557         }
 558     }
 559     else // if rtol is not given atol will be used as a scalar.
 560     {
 561         if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
 562         {
 563             double dblSrc = 1.e-9;
 564             int iSize = pDblAtol->getSize();
 565             int iOne = 1;
 566             int iZero = 0;
 567
 568             rtol = (double*)MALLOC(iSize * sizeof(double));
 569             C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
 570         }
 571         else
 572         {
 573             rtol    = (double*)MALLOC(sizeof(double));
 574             *rtol   = 1.e-9;
 575         }
 576     }
 577
 578     // Compute rwork, iwork size.
 579     // Create them.
 580     int iworksize   = 20 + pDblX0->getRows();
 581     int rworksize   = 0;
 582     int* iwork      = NULL;
 583     double* rwork   = NULL;
 584
 585     if (info[5] == 0)
 586     {
 587         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + pDblX0->getRows() * pDblX0->getRows();
 588     }
 589     else if (info[4] == 1)
 590     {
 591         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows();
 592     }
 593     else if (info[4] == 0)
 594     {
 595         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
 596     }
 597
 598     iwork = (int*)MALLOC(iworksize * sizeof(int));
 599     rwork = (double*)MALLOC(rworksize * sizeof(double));
 600
 601     if (pDblHd != NULL)
 602     {
 603         if (iworksize + rworksize != pDblHd->getSize())
 604         {
 605             Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "dassl", in.size(), iworksize + rworksize);
 606             DifferentialEquation::removeDifferentialEquationFunctions();
 607             FREE(pdYdotData);
 608             FREE(pdYData);
 609             FREE(YSize);
 610             FREE(iwork);
 611             FREE(rwork);
 612             if (pDblAtol == NULL || pDblAtol->isScalar())
 613             {
 614                 FREE(atol);
 615             }
 616             if (pDblRtol == NULL || pDblRtol->isScalar())
 617             {
 618                 FREE(rtol);
 619             }
 620             return types::Function::Error;
 621         }
 622
 623         C2F(dcopy)(&rworksize, pDblHd->get(), &one, rwork, &one);
 624
 625         for (int i = 0; i < iworksize; i++)
 626         {
 627             iwork[i] = (int)pDblHd->get(rworksize + i);
 628         }
 629
 630         info[0] = 1;
 631     }
 632
 633     if (info[3] == 1)
 634     {
 635         rwork[0] = tstop;
 636     }
 637
 638     if (info[6] == 1)
 639     {
 640         rwork[1] = maxstep;
 641     }
 642
 643     if (info[7] == 1)
 644     {
 645         rwork[2] = stepin;
 646     }
 647
 648     if (info[5] == 1)
 649     {
 650         iwork[0] = ml;
 651         iwork[1] = mu;
 652     }
 653
 654     // *** Perform operation. ***
 655     std::list<types::Double*> lpDblOut;
 656     int size = pDblX0->getRows();
 657     int rowsOut = 1 + pDblX0->getRows() * 2;
 658     int iret = 0;
 659
 660     for (int i = 0; i < pDblT->getSize(); i++)
 661     {
 662         types::Double* pDblOut = new types::Double(rowsOut, 1);
 663         lpDblOut.push_back(pDblOut);
 664
 665         double t = pDblT->get(i);
 666         int pos  = 0;
 667         pDblOut->set(pos, t);
 668
 669         if (t == t0)
 670         {
 671             pos++;
 672             C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
 673             pos += pDblX0->getRows();
 674             C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
 675
 676             continue;
 677         }
 678
 679         try
 680         {
 681             C2F(dassl)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t, info, rtol, atol, &idid, rwork, &rworksize, iwork, &iworksize, &rpar, &ipar, dassl_jac);
 682
 683             iret = checkError(idid, "dassl");
 684             if (iret == 1) // error
 685             {
 686                 Scierror(999, _("%s: dassl return with state %d.\n"), "dassl", idid);
 687             }
 688         }
 689         catch (ast::ScilabError &e)
 690         {
 691             char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
 692             sciprint(_("%s: exception caught in '%s' subroutine.\n"), "dassl", "dassl");
 693             Scierror(999, pstrMsg);
 694             // set iret to 1 for FREE allocated memory
 695             iret = 1;
 696         }
 697
 698         if (iret == 1) // error
 699         {
 700             lpDblOut.clear();
 701             DifferentialEquation::removeDifferentialEquationFunctions();
 702             FREE(pdYdotData);
 703             FREE(pdYData);
 704             FREE(YSize);
 705             FREE(iwork);
 706             FREE(rwork);
 707             if (pDblAtol == NULL || pDblAtol->isScalar())
 708             {
 709                 FREE(atol);
 710             }
 711             if (pDblRtol == NULL || pDblRtol->isScalar())
 712             {
 713                 FREE(rtol);
 714             }
 715             return types::Function::Error;
 716         }
 717
 718         pos++;
 719         C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
 720         pos += size;
 721         C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
 722
 723         if (iret == 2) // warning
 724         {
 725             break;
 726         }
 727
 728         // iret == 0
 729         if (idid == 1)
 730         {
 731             pDblOut->set(0, t0);
 732             i--;
 733         }
 734         else if (idid == -2)
 735         {
 736             t0 = t;
 737             i--;
 738         }
 739         else
 740         {
 741             t0 = t;
 742         }
 743
 744         info[0] = 1;
 745     }
 746
 747     // *** Return result in Scilab. ***
 748     types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
 749
 750     int sizeOfList = (int)lpDblOut.size();
 751     for (int i = 0; i < sizeOfList; i++)
 752     {
 753         int pos = i * rowsOut;
 754         C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &one, pDblOut->get() + pos, &one);
 755         lpDblOut.pop_front();
 756     }
 757     out.push_back(pDblOut);
 758
 759     if (_iRetCount == 2)
 760     {
 761         types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
 762         C2F(dcopy)(&rworksize, rwork, &one, pDblHdOut->get(), &one);
 763
 764         for (int i = 0; i < iworksize; i++)
 765         {
 766             pDblHdOut->set(rworksize + i, (double)iwork[i]);
 767         }
 768
 769         out.push_back(pDblHdOut);
 770     }
 771
 772     // *** FREE. ***
 773     if (pDblAtol == NULL || pDblAtol->isScalar())
 774     {
 775         FREE(atol);
 776     }
 777
 778     if (pDblRtol == NULL || pDblRtol->isScalar())
 779     {
 780         FREE(rtol);
 781     }
 782
 783     FREE(pdYData);
 784     FREE(pdYdotData);
 785     FREE(YSize);
 786     FREE(rwork);
 787     FREE(iwork);
 788
 789     DifferentialEquation::removeDifferentialEquationFunctions();
 790
 791     return types::Function::OK;
 792 }
 793