5a763444e7b24a3814fcf78b7bd0aecf5a08548c
[scilab.git] / 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
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 "checkodeerror.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->getSize());
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 = 7 * 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 = 25 * pDblX0->getRows() + 1;
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.005;
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 = checkDasslError(idid);
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