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