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