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