scilab/modules/differential_equations/sci_gateway/cpp/sci_daskr.cpp

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
   4 * Copyright (C) 2013 - Scilab Enterprises - Cedric DELAMARRE
   5 *
   6 * This file must be used under the terms of the CeCILL.
   7 * This source file is licensed as described in the file COPYING, which
   8 * you should have received as part of this distribution.  The terms
   9 * are also available at
  10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  11 *
  12 */
  13 /*--------------------------------------------------------------------------*/
  14
  15 #include "differential_equations_gw.hxx"
  16 #include "function.hxx"
  17 #include "double.hxx"
  18 #include "string.hxx"
  19 #include "list.hxx"
  20 #include "callable.hxx"
  21 #include "differentialequationfunctions.hxx"
  22 #include "runvisitor.hxx"
  23 #include "checkodeerror.hxx"
  24
  25 extern "C"
  26 {
  27 #include "sci_malloc.h"
  28 #include "localization.h"
  29 #include "Scierror.h"
  30 #include "sciprint.h"
  31 #include "scifunctions.h"
  32 #include "elem_common.h"
  33 }
  34
  35 /*--------------------------------------------------------------------------*/
  36 //  This code was inspired by sci_f_daskr.f from the same folder.
  37 //  daskr() has more optional parameters than daskr(), so the gateway
  38 //  has to adapt to the new options
  39 /*--------------------------------------------------------------------------*/
  40 //[y0,nvs[,hotdata]]=daskr(y0,t0,t1[,atol[,rtol]],res[,jac],nh,h[,info[,psol][,pjac]][,hotdata])
  41 types::Function::ReturnValue sci_daskr(types::typed_list &in, int _iRetCount, types::typed_list &out)
  42 {
  43     // input args
  44     types::Double* pDblX0   = NULL;
  45     types::Double* pDblT0   = NULL;
  46     types::Double* pDblT    = NULL;
  47     types::Double* pDblRtol = NULL;
  48     types::Double* pDblAtol = NULL;
  49     types::Double* pDblHd   = NULL;
  50     types::Double* pDblNg   = NULL;
  51     types::Double* pDblE7   = NULL; // 7th element of list info
  52     types::Double* pDblE12  = NULL; // 12th element of list info
  53
  54
  55     // x0 = [y0, ydot0]
  56     double* pdYData         = NULL; // contain y0 following by all args data in list case.
  57     double* pdYdotData      = NULL;
  58     int sizeOfpdYData       = 0;
  59
  60     int sizeOfYSize = 1;
  61     int* YSize      = NULL;    // YSize(1) = size of y0,
  62     // YSize(n) = size of Args(n) in list case.
  63     int iPos    = 0;
  64     int iOne    = 1;
  65
  66     int info[20]    = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  67
  68     double tstop    = 0;
  69     double maxstep  = 0;
  70     double stepin   = 0;
  71     double epli     = 0;
  72     double dTemp    = 0;
  73     double stptol   = 0;
  74     double epinit   = 0;
  75
  76     int ng      = 0;
  77     int mu      = 0;
  78     int ml      = 0;
  79     int maxl    = 0;
  80     int kmp     = 0;
  81     int nrmax   = 0;
  82     int mxnit   = 0;
  83     int mxnj    = 0;
  84     int mxnh    = 0;
  85     int lsoff   = 0;
  86     int LID     = 0;
  87     int LENIWP  = 0;
  88     int LENWP   = 0;
  89
  90     // Indicate if the function is given.
  91     bool bFuncF     = false;
  92     bool bFuncJac   = false;
  93     bool bFuncG     = false;
  94     bool bFuncPsol  = false;
  95     bool bFuncPjac  = false;
  96
  97     // Indicate if info list is given.
  98     bool bListInfo  = false;
  99
 100     // *** check the minimal number of input args. ***
 101     if ((int)in.size() < 6 || (int)in.size() > 13)
 102     {
 103         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "daskr", 6, 13);
 104         return types::Function::Error;
 105     }
 106
 107     // *** check number of output args ***
 108     if (_iRetCount != 3 && _iRetCount != 2)
 109     {
 110         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "daskr", 2, 3);
 111         return types::Function::Error;
 112     }
 113
 114     // *** check type of input args and get it. ***
 115     // x0 = [y0, yd0]
 116     if (in[iPos]->isDouble() == false)
 117     {
 118         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
 119         return types::Function::Error;
 120     }
 121
 122     pDblX0 = in[iPos]->getAs<types::Double>();
 123
 124     if (pDblX0->isComplex())
 125     {
 126         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
 127         return types::Function::Error;
 128     }
 129
 130     if (pDblX0->getCols() > 2)
 131     {
 132         Scierror(999, _("%s: Wrong size for input argument #%d: A real matrix with %d to %d colomn(s) expected.\n"), "daskr", iPos + 1, 1, 2);
 133         return types::Function::Error;
 134     }
 135
 136     if (pDblX0->getCols() == 1)
 137     {
 138         info[10] = 1;
 139     }
 140
 141     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
 142     *YSize = pDblX0->getRows();
 143
 144     pdYData = (double*)MALLOC(*YSize * sizeof(double));
 145     pdYdotData = (double*)MALLOC(*YSize * sizeof(double));
 146
 147     C2F(dcopy)(YSize, pDblX0->get(), &iOne, pdYData, &iOne);
 148     if (pDblX0->getCols() == 2)
 149     {
 150         C2F(dcopy)(YSize, pDblX0->get() + *YSize, &iOne, pdYdotData, &iOne);
 151     }
 152     else
 153     {
 154         memset(pdYdotData, 0x00, *YSize);
 155     }
 156
 157     // t0
 158     iPos++;
 159     if (in[iPos]->isDouble() == false)
 160     {
 161         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
 162         return types::Function::Error;
 163     }
 164
 165     pDblT0 = in[iPos]->getAs<types::Double>();
 166
 167     if (pDblT0->isScalar() == false)
 168     {
 169         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
 170         return types::Function::Error;
 171     }
 172
 173     // t
 174     iPos++;
 175     if (in[iPos]->isDouble() == false)
 176     {
 177         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
 178         return types::Function::Error;
 179     }
 180
 181     pDblT = in[iPos]->getAs<types::Double>();
 182
 183     if (pDblT->isComplex())
 184     {
 185         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
 186         return types::Function::Error;
 187     }
 188
 189     // get next inputs
 190     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"daskr");
 191     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
 192
 193     deFunctionsManager->setOdeYRows(pDblX0->getRows());
 194
 195     for (iPos++; iPos < in.size(); iPos++)
 196     {
 197         if (in[iPos]->isDouble())
 198         {
 199             if (pDblAtol == NULL && bFuncF == false)
 200             {
 201                 pDblAtol = in[iPos]->getAs<types::Double>();
 202                 if (pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
 203                 {
 204                     Scierror(267, _("%s: Wrong size for input argument #%d: A scalar or a matrix of size %d expected.\n"), "daskr", iPos + 1, pDblX0->getRows());
 205                     DifferentialEquation::removeDifferentialEquationFunctions();
 206                     FREE(pdYdotData);
 207                     FREE(pdYData);
 208                     FREE(YSize);
 209                     return types::Function::Error;
 210                 }
 211             }
 212             else if (pDblRtol == NULL && bFuncF == false)
 213             {
 214                 pDblRtol = in[iPos]->getAs<types::Double>();
 215                 if (pDblAtol->getSize() != pDblRtol->getSize())
 216                 {
 217                     Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "daskr", iPos + 1, pDblX0->getRows());
 218                     DifferentialEquation::removeDifferentialEquationFunctions();
 219                     FREE(pdYdotData);
 220                     FREE(pdYData);
 221                     FREE(YSize);
 222                     return types::Function::Error;
 223                 }
 224             }
 225             else if (pDblNg == NULL && bFuncF == true)
 226             {
 227                 pDblNg = in[iPos]->getAs<types::Double>();
 228                 if (pDblNg->isScalar() == false)
 229                 {
 230                     Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
 231                     DifferentialEquation::removeDifferentialEquationFunctions();
 232                     FREE(pdYdotData);
 233                     FREE(pdYData);
 234                     FREE(YSize);
 235                     return types::Function::Error;
 236                 }
 237                 ng = (int)pDblNg->get(0);
 238             }
 239             else if (pDblHd == NULL && bFuncF == true)
 240             {
 241                 pDblHd = in[iPos]->getAs<types::Double>();
 242                 if ((int)in.size() != iPos + 1)
 243                 {
 244                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", iPos + 1);
 245                     DifferentialEquation::removeDifferentialEquationFunctions();
 246                     FREE(pdYdotData);
 247                     FREE(pdYData);
 248                     FREE(YSize);
 249                     return types::Function::Error;
 250                 }
 251             }
 252             else
 253             {
 254                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "daskr", iPos + 1);
 255                 DifferentialEquation::removeDifferentialEquationFunctions();
 256                 FREE(pdYdotData);
 257                 FREE(pdYData);
 258                 FREE(YSize);
 259                 return types::Function::Error;
 260             }
 261         }
 262         else if (in[iPos]->isCallable())
 263         {
 264             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
 265             if (bFuncF == false)
 266             {
 267                 deFunctionsManager->setFFunction(pCall);
 268                 bFuncF = true;
 269             }
 270             else if (bFuncJac == false && pDblNg == NULL)
 271             {
 272                 deFunctionsManager->setJacFunction(pCall);
 273                 bFuncJac = true;
 274             }
 275             else if (bFuncG == false && pDblNg)
 276             {
 277                 deFunctionsManager->setGFunction(pCall);
 278                 bFuncG = true;
 279             }
 280             else if (bFuncG && bFuncPsol == false)
 281             {
 282                 deFunctionsManager->setPsolFunction(pCall);
 283                 bFuncPsol = true;
 284             }
 285             else if (bFuncPsol && bFuncPjac == false)
 286             {
 287                 deFunctionsManager->setPjacFunction(pCall);
 288                 bFuncPjac = true;
 289             }
 290             else
 291             {
 292                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
 293                 DifferentialEquation::removeDifferentialEquationFunctions();
 294                 FREE(pdYdotData);
 295                 FREE(pdYData);
 296                 FREE(YSize);
 297                 return types::Function::Error;
 298             }
 299         }
 300         else if (in[iPos]->isString())
 301         {
 302             types::String* pStr = in[iPos]->getAs<types::String>();
 303             bool bOK = false;
 304
 305             if (bFuncF == false)
 306             {
 307                 bOK = deFunctionsManager->setFFunction(pStr);
 308                 bFuncF = true;
 309             }
 310             else if (bFuncJac == false && pDblNg == NULL)
 311             {
 312                 bOK = deFunctionsManager->setJacFunction(pStr);
 313                 bFuncJac = true;
 314             }
 315             else if (bFuncG == false && pDblNg)
 316             {
 317                 bOK = deFunctionsManager->setGFunction(pStr);
 318                 bFuncG = true;
 319             }
 320             else if (bFuncG && bFuncPsol == false)
 321             {
 322                 bOK = deFunctionsManager->setPsolFunction(pStr);
 323                 bFuncPsol = true;
 324             }
 325             else if (bFuncPsol && bFuncPjac == false)
 326             {
 327                 bOK = deFunctionsManager->setPjacFunction(pStr);
 328                 bFuncPjac = true;
 329             }
 330             else
 331             {
 332                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
 333                 DifferentialEquation::removeDifferentialEquationFunctions();
 334                 FREE(pdYdotData);
 335                 FREE(pdYData);
 336                 FREE(YSize);
 337                 return types::Function::Error;
 338             }
 339
 340             if (bOK == false)
 341             {
 342                 char* pst = wide_string_to_UTF8(pStr->get(0));
 343                 Scierror(50, _("%s: Subroutine not found: %s\n"), "daskr", pst);
 344                 FREE(pst);
 345                 DifferentialEquation::removeDifferentialEquationFunctions();
 346                 FREE(pdYdotData);
 347                 FREE(pdYData);
 348                 FREE(YSize);
 349                 return types::Function::Error;
 350             }
 351         }
 352         else if (in[iPos]->isList())
 353         {
 354             types::List* pList = in[iPos]->getAs<types::List>();
 355
 356             if (pList->getSize() == 0)
 357             {
 358                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, "(string empty)");
 359                 DifferentialEquation::removeDifferentialEquationFunctions();
 360                 FREE(pdYdotData);
 361                 FREE(pdYData);
 362                 FREE(YSize);
 363                 return types::Function::Error;
 364             }
 365
 366             if (bFuncF && bListInfo)
 367             {
 368                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
 369                 DifferentialEquation::removeDifferentialEquationFunctions();
 370                 FREE(pdYdotData);
 371                 FREE(pdYData);
 372                 FREE(YSize);
 373                 return types::Function::Error;
 374             }
 375
 376             // functions can be passed in a list
 377             if (pList->get(0)->isString())
 378             {
 379                 // list(string,...  function case
 380                 types::String* pStr = pList->get(0)->getAs<types::String>();
 381                 bool bOK = false;
 382                 sizeOfpdYData = *YSize;
 383
 384                 if (bFuncF == false)
 385                 {
 386                     bFuncF = true;
 387                     bOK = deFunctionsManager->setFFunction(pStr);
 388                 }
 389                 else if (bFuncJac == false && pDblNg == NULL)
 390                 {
 391                     bFuncJac = true;
 392                     bOK = deFunctionsManager->setJacFunction(pStr);
 393                 }
 394                 else if (bFuncG == false && pDblNg)
 395                 {
 396                     bFuncG = true;
 397                     bOK = deFunctionsManager->setGFunction(pStr);
 398                 }
 399                 else if (bFuncG && bFuncPsol == false)
 400                 {
 401                     bOK = deFunctionsManager->setPsolFunction(pStr);
 402                     bFuncPsol = true;
 403                 }
 404                 else if (bFuncPsol && bFuncPjac == false)
 405                 {
 406                     bOK = deFunctionsManager->setPjacFunction(pStr);
 407                     bFuncPjac = true;
 408                 }
 409
 410                 if (bOK == false)
 411                 {
 412                     char* pst = wide_string_to_UTF8(pStr->get(0));
 413                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, pst);
 414                     FREE(pst);
 415                     DifferentialEquation::removeDifferentialEquationFunctions();
 416                     FREE(pdYdotData);
 417                     FREE(pdYData);
 418                     FREE(YSize);
 419                     return types::Function::Error;
 420                 }
 421
 422                 int* sizeTemp = YSize;
 423                 int totalSize = sizeOfpdYData;
 424
 425                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
 426                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
 427
 428                 std::vector<types::Double*> vpDbl;
 429                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 430                 {
 431                     if (pList->get(iter + 1)->isDouble() == false)
 432                     {
 433                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "daskr", iPos + 1, iter + 1);
 434                         DifferentialEquation::removeDifferentialEquationFunctions();
 435                         FREE(pdYdotData);
 436                         FREE(pdYData);
 437                         FREE(YSize);
 438                         return types::Function::Error;
 439                     }
 440
 441                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
 442                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
 443                     totalSize += YSize[sizeOfYSize + iter];
 444                 }
 445
 446                 double* pdYDataTemp = pdYData;
 447                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
 448                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &iOne, pdYData, &iOne);
 449
 450                 int position = sizeOfpdYData;
 451                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
 452                 {
 453                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &iOne, &pdYData[position], &iOne);
 454                     position += vpDbl[iter]->getSize();
 455                 }
 456                 vpDbl.clear();
 457                 sizeOfpdYData = totalSize;
 458                 sizeOfYSize += pList->getSize() - 1;
 459                 FREE(pdYDataTemp);
 460                 FREE(sizeTemp);
 461             }
 462             else if (pList->get(0)->isCallable())
 463             {
 464                 // list(macro,...  function case
 465                 if (bFuncF == false)
 466                 {
 467                     bFuncF = true;
 468                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
 469                     for (int iter = 1; iter < pList->getSize(); iter++)
 470                     {
 471                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
 472                     }
 473                 }
 474                 else if (bFuncJac == false && pDblNg == NULL)
 475                 {
 476                     bFuncJac = true;
 477                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
 478                     for (int iter = 1; iter < pList->getSize(); iter++)
 479                     {
 480                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
 481                     }
 482                 }
 483                 else if (bFuncG == false && pDblNg)
 484                 {
 485                     bFuncG = true;
 486                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
 487                     for (int iter = 1; iter < pList->getSize(); iter++)
 488                     {
 489                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
 490                     }
 491                 }
 492                 else if (bFuncG && bFuncPsol == false)
 493                 {
 494                     bFuncPsol = true;
 495                     deFunctionsManager->setPsolFunction(pList->get(0)->getAs<types::Callable>());
 496                     for (int iter = 1; iter < pList->getSize(); iter++)
 497                     {
 498                         deFunctionsManager->setPsolArgs(pList->get(iter)->getAs<types::InternalType>());
 499                     }
 500                 }
 501                 else if (bFuncPsol && bFuncPjac == false)
 502                 {
 503                     bFuncPjac = true;
 504                     deFunctionsManager->setPjacFunction(pList->get(0)->getAs<types::Callable>());
 505                     for (int iter = 1; iter < pList->getSize(); iter++)
 506                     {
 507                         deFunctionsManager->setPjacArgs(pList->get(iter)->getAs<types::InternalType>());
 508                     }
 509                 }
 510             }
 511             else if (pList->get(0)->isDouble() && bFuncF == true)
 512             {
 513                 // list(double,... then this list is the info argument
 514                 if (pList->getSize() != 14)
 515                 {
 516                     Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "daskr", iPos + 1, 7);
 517                     DifferentialEquation::removeDifferentialEquationFunctions();
 518                     FREE(pdYdotData);
 519                     FREE(pdYData);
 520                     FREE(YSize);
 521                     return types::Function::Error;
 522                 }
 523
 524                 for (int i = 0; i < 14; i++) // info = list([],0,[],[],[],0,[],0,[],0,[],[],[],1)
 525                 {
 526                     if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)))
 527                     {
 528                         if (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)
 529                         {
 530                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "daskr", iPos + 1, i);
 531                         }
 532                         else
 533                         {
 534                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "daskr", iPos + 1, i);
 535                         }
 536                         DifferentialEquation::removeDifferentialEquationFunctions();
 537                         FREE(pdYdotData);
 538                         FREE(pdYData);
 539                         FREE(YSize);
 540                         return types::Function::Error;
 541                     }
 542                 }
 543
 544                 //  --   subvariable tstop(info) --
 545                 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
 546                 if (pDblTemp->getSize() != 0)
 547                 {
 548                     info[3] = 1;
 549                     tstop = pDblTemp->get(0);
 550                 }
 551
 552                 //  --   subvariable imode(info) --
 553                 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
 554
 555                 //  --   subvariable band(info) --
 556                 pDblTemp = pList->get(2)->getAs<types::Double>();
 557                 if (pDblTemp->getSize() == 2)
 558                 {
 559                     info[5] = 1;
 560                     ml = (int)pDblTemp->get(0);
 561                     mu = (int)pDblTemp->get(1);
 562                     deFunctionsManager->setMl(ml);
 563                     deFunctionsManager->setMu(mu);
 564                 }
 565                 else if (pDblTemp->getSize() != 0)
 566                 {
 567                     Scierror(267, _("%s: Wrong size for input argument #%d: Argument %d in te list must be of size %d.\n"), "daskr", iPos + 1, 3, 2);
 568                     DifferentialEquation::removeDifferentialEquationFunctions();
 569                     FREE(pdYdotData);
 570                     FREE(pdYData);
 571                     FREE(YSize);
 572                     return types::Function::Error;
 573                 }
 574
 575                 //  --   subvariable maxstep(info) --
 576                 pDblTemp = pList->get(3)->getAs<types::Double>();
 577                 if (pDblTemp->getSize() != 0)
 578                 {
 579                     info[6] = 1;
 580                     maxstep = pDblTemp->get(0);
 581                 }
 582
 583                 //  --   subvariable stepin(info) --
 584                 pDblTemp = pList->get(4)->getAs<types::Double>();
 585                 if (pDblTemp->getSize() != 0)
 586                 {
 587                     info[7] = 1;
 588                     stepin = pDblTemp->get(0);
 589                 }
 590
 591                 //  --   subvariable nonneg(info) --
 592                 info[9]  = (int)pList->get(5)->getAs<types::Double>()->get(0);
 593
 594                 //  --   subvariable consistent(info) --
 595                 pDblTemp = pList->get(6)->getAs<types::Double>();
 596                 if (pDblTemp->getSize() != 0)
 597                 {
 598                     info[10] = 1;
 599                     if (info[9] == 0 || info[9] == 2)
 600                     {
 601                         LID = 40;
 602                     }
 603                     else
 604                     {
 605                         LID = 40 + pDblX0->getRows();
 606                     }
 607                 }
 608                 pDblE7 = pDblTemp;
 609
 610                 //  --   subvariable iteration(info) --
 611                 if (pList->get(7)->getAs<types::Double>()->get(0) == 1)
 612                 {
 613                     info[11] = 1;
 614                 }
 615
 616                 //  --   subvariable defaultKrylov(info) --
 617                 pDblTemp = pList->get(8)->getAs<types::Double>();
 618                 if (pDblTemp->getSize() == 0)
 619                 {
 620                     // maxl and kmp need default values
 621                     maxl = min(5, pDblX0->getRows());
 622                     kmp = maxl;
 623                 }
 624                 else
 625                 {
 626                     // info then looks like list(..., [maxl kmp nrmax epli],...)
 627                     info[12] = 1;
 628                     maxl  = (int)pDblTemp->get(0);
 629                     kmp   = (int)pDblTemp->get(1);
 630                     nrmax = (int)pDblTemp->get(2);
 631                     epli  = pDblTemp->get(3);
 632                 }
 633
 634                 //  --   subvariable justConsistentComp(info) --
 635                 dTemp = pList->get(9)->getAs<types::Double>()->get(0);
 636                 if (dTemp)
 637                 {
 638                     // Check that info(11) = 1, meaning that the provided initial values
 639                     // are not consistent
 640                     if (info[10] == 1)
 641                     {
 642                         info[13] = 1;
 643                     }
 644                 }
 645
 646                 //  --   subvariable psolJac(info) --
 647                 dTemp = pList->get(10)->getAs<types::Double>()->get(0);
 648                 if (dTemp)
 649                 {
 650                     info[14] = 1;
 651                 }
 652
 653                 //  --   subvariable excludeAlgebraic(info) --
 654                 pDblTemp = pList->get(11)->getAs<types::Double>();
 655                 if (pDblTemp->getSize() != 0)
 656                 {
 657                     info[15] = 1;
 658                     if (info[9] == 0 || info[9] == 2)
 659                     {
 660                         LID = 40;
 661                     }
 662                     else
 663                     {
 664                         LID = 40 + pDblX0->getRows();
 665                     }
 666                 }
 667                 pDblE12 = pDblTemp;
 668
 669                 //  --   subvariable defaultHeuristic(info) --
 670                 pDblTemp = pList->get(12)->getAs<types::Double>();
 671                 if (pDblTemp->getSize() != 0)
 672                 {
 673                     // info then looks like list(..., [mxnit mxnj lsoff stptol epinit],...)
 674                     info[16] = 1;
 675                     mxnit   = (int)pDblTemp->get(0);
 676                     mxnj    = (int)pDblTemp->get(1);
 677                     mxnh    = (int)pDblTemp->get(2);
 678                     lsoff   = (int)pDblTemp->get(3);
 679                     stptol  = pDblTemp->get(4);
 680                     epinit  = pDblTemp->get(5);
 681                 }
 682
 683                 //  --   subvariable verbosity(info) --
 684                 dTemp = pList->get(13)->getAs<types::Double>()->get(0);
 685                 if (dTemp == 1 || dTemp == 2)
 686                 {
 687                     info[17] = (int)dTemp;
 688                 }
 689
 690                 bListInfo = true;
 691             }
 692             else
 693             {
 694                 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"), "daskr", iPos + 1);
 695                 DifferentialEquation::removeDifferentialEquationFunctions();
 696                 FREE(pdYdotData);
 697                 FREE(pdYData);
 698                 FREE(YSize);
 699                 return types::Function::Error;
 700             }
 701         }
 702         else
 703         {
 704             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "daskr", iPos + 1);
 705             DifferentialEquation::removeDifferentialEquationFunctions();
 706             FREE(pdYdotData);
 707             FREE(pdYData);
 708             FREE(YSize);
 709             return types::Function::Error;
 710         }
 711     }
 712
 713     if (bFuncF == false)
 714     {
 715         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 3);
 716         DifferentialEquation::removeDifferentialEquationFunctions();
 717         FREE(pdYdotData);
 718         FREE(pdYData);
 719         FREE(YSize);
 720         return types::Function::Error;
 721     }
 722
 723     if (pDblNg == NULL)
 724     {
 725         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 2);
 726         DifferentialEquation::removeDifferentialEquationFunctions();
 727         FREE(pdYdotData);
 728         FREE(pdYData);
 729         FREE(YSize);
 730         return types::Function::Error;
 731     }
 732
 733     if (bFuncG == false)
 734     {
 735         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 1);
 736         DifferentialEquation::removeDifferentialEquationFunctions();
 737         FREE(pdYdotData);
 738         FREE(pdYData);
 739         FREE(YSize);
 740         return types::Function::Error;
 741     }
 742
 743     if (bFuncJac == true)
 744     {
 745         info[4] = 1;
 746     }
 747
 748     // *** Initialization. ***
 749     double t0    = pDblT0->get(0);
 750     int idid     = 0;
 751     int maxord   = 5;
 752
 753     //compute itol and set the tolerances rtol and atol.
 754     double* rtol    = NULL;
 755     double* atol    = NULL;
 756     double rpar[2] = {0, 0 };
 757     int ipar[30];
 758     memset(ipar, 0x00, 30 * sizeof(int));
 759
 760     if (pDblAtol)
 761     {
 762         if (pDblAtol->isScalar())
 763         {
 764             atol  = (double*)MALLOC(sizeof(double));
 765             *atol = pDblAtol->get(0);
 766         }
 767         else
 768         {
 769             atol    = pDblAtol->get();
 770             info[1] = 1;
 771         }
 772     }
 773     else
 774     {
 775         atol  = (double*)MALLOC(sizeof(double));
 776         *atol = 1.e-7;
 777     }
 778
 779     if (pDblRtol)
 780     {
 781         if (pDblRtol->isScalar())
 782         {
 783             rtol  = (double*)MALLOC(sizeof(double));
 784             *rtol = pDblRtol->get(0);
 785         }
 786         else
 787         {
 788             rtol = pDblRtol->get();
 789         }
 790     }
 791     else // if rtol is not given atol will be used as a scalar.
 792     {
 793         if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
 794         {
 795             double dblSrc = 1.e-9;
 796             int iSize = pDblAtol->getSize();
 797             int iZero = 0;
 798
 799             rtol = (double*)MALLOC(iSize * sizeof(double));
 800             C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
 801         }
 802         else
 803         {
 804             rtol    = (double*)MALLOC(sizeof(double));
 805             *rtol   = 1.e-9;
 806         }
 807     }
 808
 809     // Compute rwork, iwork size.
 810     // Create them.
 811     int rworksize   = 0;
 812     int* iwork      = NULL;
 813     double* rwork   = NULL;
 814     int* root       = NULL;
 815
 816     // compute size of rwork
 817     rworksize = 60;
 818     if (info[11] == 0)
 819     {
 820         rworksize += max(maxord + 4, 7) * pDblX0->getRows() + 3 * ng;
 821         if (info[5] == 0)
 822         {
 823             // For the full (dense) JACOBIAN case
 824             rworksize += pDblX0->getRows() * pDblX0->getRows();
 825         }
 826         else if (info[4] == 1)
 827         {
 828             // For the banded user-defined JACOBIAN case
 829             rworksize += (2 * ml + mu + 1) * pDblX0->getRows();
 830         }
 831         else if (info[4] == 0)
 832         {
 833             // For the banded finite-difference-generated JACOBIAN case
 834             rworksize += (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
 835         }
 836     }
 837     else // info[11] == 1
 838     {
 839         // LENWP is the length ot the rwork segment containing
 840         // the matrix elements of the preconditioner P
 841
 842         LENWP = pDblX0->getRows() * pDblX0->getRows();
 843         rworksize += (maxord + 5) * pDblX0->getRows() + 3 * ng
 844                      + (maxl + 3 + min(1, maxl - kmp)) * pDblX0->getRows()
 845                      + (maxl + 3) * maxl + 1 + LENWP;
 846     }
 847
 848     if (info[15] == 1)
 849     {
 850         rworksize += pDblX0->getRows();
 851     }
 852
 853     // compute size of iwork
 854     // liw = 40, then augment size according to the case
 855     int iworksize = 40;
 856     if (info[11] == 0)
 857     {
 858         iworksize += pDblX0->getRows();
 859     }
 860     else // info[11] == 1
 861     {
 862         // LENIWP is the length ot the iwork segment containing
 863         // the matrix indexes of the preconditioner P
 864         // (compressed sparse row format)
 865         LENIWP = 2 * pDblX0->getRows() * pDblX0->getRows();
 866         iworksize += LENIWP;
 867     }
 868
 869     if (info[9] == 1 || info[9] == 3)
 870     {
 871         iworksize += pDblX0->getRows();
 872     }
 873
 874     if (info[10] == 1 || info[15] == 1)
 875     {
 876         iworksize += pDblX0->getRows();
 877     }
 878
 879     iwork = (int*)MALLOC(iworksize * sizeof(int));
 880     rwork = (double*)MALLOC(rworksize * sizeof(double));
 881     root  = (int*)MALLOC(ng * sizeof(int));
 882
 883     if (pDblHd != NULL)
 884     {
 885         if (iworksize + rworksize != pDblHd->getSize())
 886         {
 887             Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "daskr", in.size(), iworksize + rworksize);
 888             DifferentialEquation::removeDifferentialEquationFunctions();
 889             FREE(pdYdotData);
 890             FREE(pdYData);
 891             FREE(YSize);
 892             FREE(iwork);
 893             FREE(rwork);
 894             if (pDblAtol == NULL || pDblAtol->isScalar())
 895             {
 896                 FREE(atol);
 897             }
 898             if (pDblRtol == NULL || pDblRtol->isScalar())
 899             {
 900                 FREE(rtol);
 901             }
 902             return types::Function::Error;
 903         }
 904
 905         C2F(dcopy)(&rworksize, pDblHd->get(), &iOne, rwork, &iOne);
 906
 907         for (int i = 0; i < iworksize; i++)
 908         {
 909             iwork[i] = (int)pDblHd->get(rworksize + i);
 910         }
 911
 912         info[0] = 1;
 913     }
 914
 915     if (info[3] == 1)
 916     {
 917         rwork[0] = tstop;
 918     }
 919
 920     if (info[6] == 1)
 921     {
 922         rwork[1] = maxstep;
 923     }
 924
 925     if (info[7] == 1)
 926     {
 927         rwork[2] = stepin;
 928     }
 929
 930     if (info[5] == 1)
 931     {
 932         iwork[0] = ml;
 933         iwork[1] = mu;
 934     }
 935
 936     if (info[10] == 1)
 937     {
 938         for (int i = 0; i < pDblX0->getRows(); i++)
 939         {
 940             iwork[LID + i] = pDblE7->get(i);
 941         }
 942     }
 943
 944     if (info[11] == 1)
 945     {
 946         iwork[26] = LENWP;
 947         iwork[27] = LENIWP;
 948     }
 949
 950     if (info[12] == 1)
 951     {
 952         iwork[23] = maxl;
 953         iwork[24] = kmp;
 954         iwork[25] = nrmax;
 955         rwork[10] = epli;
 956     }
 957
 958     if (info[14] == 1)
 959     {
 960         // Set ipar and rper
 961         ipar[0] = pDblX0->getRows();
 962         ipar[1] = pDblX0->getRows();
 963         ipar[2] = 2;
 964         ipar[3] = 2;
 965         ipar[4] = 1;
 966         rpar[0] = 0.005;
 967         rpar[1] = 0.05;
 968     }
 969
 970     if (info[15] == 1)
 971     {
 972         for (int i = 0; i < pDblX0->getRows(); i++)
 973         {
 974             iwork[LID + i] = pDblE12->get(i);
 975         }
 976     }
 977
 978     if (info[16] == 1)
 979     {
 980         iwork[31] = mxnit;
 981         iwork[32] = mxnj;
 982         iwork[33] = mxnh;
 983         iwork[34] = lsoff;
 984         rwork[13] = stptol;
 985         rwork[14] = epinit;
 986     }
 987
 988     iwork[16] = rworksize;
 989     iwork[17] = iworksize;
 990
 991     // *** Perform operation. ***
 992     std::list<types::Double*> lpDblOut;
 993     int size = pDblX0->getRows();
 994     int rowsOut = 1 + pDblX0->getRows() * 2;
 995     int iret = 0;
 996     int ididtmp = 0;
 997
 998     for (int i = 0; i < pDblT->getSize(); i++)
 999     {
1000         types::Double* pDblOut = new types::Double(rowsOut, 1);
1001         lpDblOut.push_back(pDblOut);
1002
1003         double t = pDblT->get(i);
1004         int pos  = 0;
1005         pDblOut->set(pos, t);
1006
1007         if (t == t0)
1008         {
1009             pos++;
1010             C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1011             pos += pDblX0->getRows();
1012             C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1013
1014             continue;
1015         }
1016
1017         try
1018         {
1019             // SUBROUTINE DDASKR(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
1020             //                  IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
1021             //                  RT, NRT, JROOT)
1022             C2F(ddaskr)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t,
1023                         info, rtol, atol, &idid, rwork, &rworksize,
1024                         iwork, &iworksize, rpar, ipar,
1025                         info[14] == 1 ? (void*)daskr_pjac : (void*)dassl_jac,
1026                         daskr_psol, dasrt_g, &ng, root);
1027
1028             // values of idid says the same things that in ddasrt function,
1029             // except these two values.
1030             ididtmp = idid;
1031             // The integration to TSTOP was successfully completed (T=TSTOP)
1032             if (idid == 4)
1033             {
1034                 idid = 2;
1035             }
1036
1037             // One or more root found (eq to idid = 4 for dassl)
1038             if (idid == 5)
1039             {
1040                 idid = 4;
1041             }
1042
1043             if (idid == 6)
1044             {
1045                 idid = 3;
1046             }
1047
1048             // use the same error function that dasrt
1049             iret = checkError(idid, "daskr");
1050             if (iret == 1) // error
1051             {
1052                 Scierror(999, _("%s: ddaskr return with state %d.\n"), "daskr", ididtmp);
1053             }
1054         }
1055         catch (ast::ScilabError &e)
1056         {
1057             char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1058             sciprint(_("%s: exception caught in '%s' subroutine.\n"), "daskr", "ddaskr");
1059             Scierror(999, pstrMsg);
1060             // set iret to 1 for free allocated memory
1061             iret = 1;
1062         }
1063
1064         // free memory when error occurred
1065         if (iret == 1)
1066         {
1067             lpDblOut.clear();
1068             DifferentialEquation::removeDifferentialEquationFunctions();
1069             FREE(pdYdotData);
1070             FREE(pdYData);
1071             FREE(YSize);
1072             FREE(iwork);
1073             FREE(rwork);
1074             FREE(root);
1075             if (pDblAtol == NULL || pDblAtol->isScalar())
1076             {
1077                 FREE(atol);
1078             }
1079             if (pDblRtol == NULL || pDblRtol->isScalar())
1080             {
1081                 FREE(rtol);
1082             }
1083             return types::Function::Error;
1084         }
1085
1086         pos++;
1087         C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1088         pos += size;
1089         C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1090
1091         if (iret == 2) // warning
1092         {
1093             if (ididtmp == -1 || ididtmp == 5 || ididtmp == 6)
1094             {
1095                 pDblOut->set(0, t0);
1096             }
1097             break;
1098         }
1099
1100         // iret == 0, check idid flag
1101         if (idid == 1)
1102         {
1103             pDblOut->set(0, t0);
1104             i--;
1105         }
1106         else if (idid == -2)
1107         {
1108             t0 = t;
1109             i--;
1110         }
1111         else
1112         {
1113             t0 = t;
1114         }
1115
1116         info[0] = 1;
1117     }
1118
1119     // *** Return result in Scilab. ***
1120     types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
1121
1122     int sizeOfList = (int)lpDblOut.size();
1123     for (int i = 0; i < sizeOfList; i++)
1124     {
1125         int pos = i * rowsOut;
1126         C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &iOne, pDblOut->get() + pos, &iOne);
1127         lpDblOut.pop_front();
1128     }
1129     out.push_back(pDblOut);
1130
1131     int sizeOfRoot = 1;
1132     for (int i = 0; i < ng; i++)
1133     {
1134         if (root[i])
1135         {
1136             sizeOfRoot++;
1137         }
1138     }
1139     types::Double* pDblRoot = new types::Double(1, sizeOfRoot);
1140     pDblRoot->set(0, t0);
1141     int j = 0;
1142     for (int i = 0; i < ng; i++)
1143     {
1144         if (root[i])
1145         {
1146             j++;
1147             pDblRoot->set(j, i + 1);
1148         }
1149     }
1150     out.push_back(pDblRoot);
1151
1152     if (_iRetCount == 3)
1153     {
1154         types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
1155         C2F(dcopy)(&rworksize, rwork, &iOne, pDblHdOut->get(), &iOne);
1156
1157         for (int i = 0; i < iworksize; i++)
1158         {
1159             pDblHdOut->set(rworksize + i, (double)iwork[i]);
1160         }
1161
1162         out.push_back(pDblHdOut);
1163     }
1164
1165     // *** FREE memory. ***
1166     if (pDblAtol == NULL || pDblAtol->isScalar())
1167     {
1168         FREE(atol);
1169     }
1170
1171     if (pDblRtol == NULL || pDblRtol->isScalar())
1172     {
1173         FREE(rtol);
1174     }
1175
1176     FREE(pdYData);
1177     FREE(pdYdotData);
1178     FREE(YSize);
1179     FREE(rwork);
1180     FREE(iwork);
1181     FREE(root);
1182
1183     DifferentialEquation::removeDifferentialEquationFunctions();
1184
1185     return types::Function::OK;
1186 }
1187