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