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