Management of errors in differential equations functions.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_ode.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 "localization.h"
26 #include "Scierror.h"
27 #include "scifunctions.h"
28 #include "elem_common.h"
29 #include "warningmode.h"
30 #include "sciprint.h"
31 #include "common_structure.h"
32 #include "checkodeerror.h"
33 #include "MALLOC.h"
34 }
35 /*--------------------------------------------------------------------------*/
36
37 types::Function::ReturnValue sci_ode(types::typed_list &in, int _iRetCount, types::typed_list &out)
38 {
39     // Methode
40     types::String* pStrType     = NULL;
41     wchar_t* wcsType            = L"lsoda";
42     int meth                    = 0;
43
44     // y0
45     //    types::Polynom* pPolyY0     = NULL;
46     types::Double* pDblY0       = NULL;
47     double* pdYData             = NULL; // contain y0 following by all args data in list case.
48     int sizeOfpdYData           = 0;
49
50     // Other input args
51     types::Double* pDblT0       = NULL;
52     types::Double* pDblT        = NULL;
53     types::Double* pDblRtol     = NULL;
54     types::Double* pDblAtol     = NULL;
55     types::Double* pDblNg       = NULL;
56     types::Double* pDblW        = NULL;
57     types::Double* pDblIw       = NULL;
58
59     // %ODEOPTIONS
60     types::Double* pDblOdeOptions = NULL;
61
62     // Result
63     types::Double* pDblYOut     = NULL;
64     //    types::Double* pDblTOut     = NULL;
65     //    types::Polynom* pPolyYOut   = NULL;
66
67     // Indicate if the function is given.
68     bool bFuncF     = false;
69     bool bFuncJac   = false;
70     bool bFuncG     = false;
71
72     int iPos        = 0; // Position in types::typed_list in
73
74     int sizeOfYSize = 1;
75     int* YSize      = NULL;    // YSize(1) = size of y0,
76     // YSize(n) = size of Args(n) in list case.
77
78     C2F(eh0001).mesflg = 1; // flag to control printing of error messages in lapack routine.
79     // 1 means print, 0 means no printing.
80     C2F(eh0001).lunit = 6;  // 6 = stdout
81
82     int one = 1; // used in C2F(dcopy)
83
84     // For root methode
85     int* jroot = NULL;
86
87     // *** check the minimal number of input args. ***
88     if (in.size() < 4)
89     {
90         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 4);
91         return types::Function::Error;
92     }
93
94     // *** Get the methode. ***
95     if (in[0]->isString())
96     {
97         pStrType = in[0]->getAs<types::String>();
98         wcsType = pStrType->get(0);
99         iPos++;
100     }
101
102     if (iPos)
103     {
104         if (wcscmp(wcsType, L"adams") == 0)
105         {
106             meth = 1;
107         }
108         else if (wcscmp(wcsType, L"stiff") == 0)
109         {
110             meth = 2;
111         }
112         else if (wcscmp(wcsType, L"root") == 0)
113         {
114             meth = 3;
115         }
116         else if (wcscmp(wcsType, L"discrete") == 0)
117         {
118             meth = 4;
119         }
120         else if (wcscmp(wcsType, L"rk") == 0)
121         {
122             meth = 5;
123         }
124         else if (wcscmp(wcsType, L"rkf") == 0)
125         {
126             meth = 6;
127         }
128         else if (wcscmp(wcsType, L"fix") == 0)
129         {
130             meth = 7;
131         }
132         else
133         {
134             Scierror(999, _("%s: Wrong value for input argument #%d: It must be one of the following strings: adams, stiff, rk, rkf, fix, root or discrete.\n"), "ode", 1);
135             return types::Function::Error;
136         }
137     }
138
139     // *** check number of output args according the methode. ***
140     if (meth < 3)
141     {
142         if (_iRetCount != 1 && _iRetCount != 3)
143         {
144             Scierror(78, _("%s: Wrong number of output argument(s): %d or %d expected.\n"), "ode", 1, 3);
145             return types::Function::Error;
146         }
147     }
148     else if (meth == 3)
149     {
150         if (_iRetCount == 3 || _iRetCount > 4)
151         {
152             Scierror(78, _("%s: Wrong number of output argument(s): %d, %d or %d expected.\n"), "ode", 1, 2, 4);
153             return types::Function::Error;
154         }
155     }
156     else // meth > 3
157     {
158         if (_iRetCount != 1)
159         {
160             Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "ode", 1);
161             return types::Function::Error;
162         }
163     }
164
165     // *** check type of input args and get it. ***
166     // y0
167     if (in[iPos]->isDouble())
168     {
169         pDblY0 = in[iPos]->getAs<types::Double>();
170         if (pDblY0->isComplex())
171         {
172             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "ode", iPos + 1);
173             return types::Function::Error;
174         }
175     }
176     /*
177         else if(in[iPos]->isPoly())
178         {
179             pPolyY0 = in[iPos]->getAs<types::Polynom>();
180             if(pPolyY0->isComplex())
181             {
182                 Scierror(999, _("%s: Wrong type for input argument #%d: A real polynom expected.\n"), "ode", iPos+1);
183                 return types::Function::Error;
184             }
185
186             if(pPolyY0->isScalar() == false)
187             {
188                 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar polynom expected.\n"), "ode", iPos+1);
189                 return types::Function::Error;
190             }
191
192             double* dbl = (double*)MALLOC(pPolyY0->getCoef()->getSize() * sizeof(double));
193             vTransposeRealMatrix(   pPolyY0->getCoef()->get(),
194                                     pPolyY0->getCoef()->getRows(),
195                                     pPolyY0->getCoef()->getCols(),
196                                     dbl);
197
198             pDblY0 = new types::Double(pPolyY0->getCoef()->getCols(), pPolyY0->getCoef()->getRows());
199             pDblY0->set(dbl);
200         }
201     */
202     else
203     {
204         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
205         return types::Function::Error;
206     }
207
208     // t0
209     iPos++;
210     if (in[iPos]->isDouble() == false)
211     {
212         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
213         return types::Function::Error;
214     }
215
216     pDblT0 = in[iPos]->getAs<types::Double>();
217
218     if (pDblT0->isScalar() == false)
219     {
220         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "ode", iPos + 1);
221         return types::Function::Error;
222     }
223
224     // t
225     iPos++;
226     if (in[iPos]->isDouble() == false)
227     {
228         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
229         return types::Function::Error;
230     }
231
232     pDblT = in[iPos]->getAs<types::Double>();
233
234     // get next inputs
235     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"ode");
236     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
237     deFunctionsManager->setOdeYRows(pDblY0->getRows());
238     deFunctionsManager->setOdeYCols(pDblY0->getCols());
239
240     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
241     *YSize = pDblY0->getSize();
242     pdYData = (double*)MALLOC(pDblY0->getSize() * sizeof(double));
243     C2F(dcopy)(YSize, pDblY0->get(), &one, pdYData, &one);
244
245     if (meth == 4)
246     {
247         if (in.size() != 5)
248         {
249             Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", 5);
250             DifferentialEquation::removeDifferentialEquationFunctions();
251             FREE(pdYData);
252             FREE(YSize);
253             return types::Function::Error;
254         }
255
256         if (in[4]->isCallable() == false && in[4]->isString() == false && in[4]->isList() == false)
257         {
258             Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", 5);
259             DifferentialEquation::removeDifferentialEquationFunctions();
260             FREE(pdYData);
261             FREE(YSize);
262             return types::Function::Error;
263         }
264     }
265
266     for (iPos++; iPos < (int)in.size(); iPos++)
267     {
268         if (in[iPos]->isDouble())
269         {
270             if (pDblRtol == NULL && bFuncF == false)
271             {
272                 pDblRtol = in[iPos]->getAs<types::Double>();
273                 if (pDblRtol->getSize() != pDblY0->getSize() && pDblRtol->isScalar() == false)
274                 {
275                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
276                     DifferentialEquation::removeDifferentialEquationFunctions();
277                     FREE(pdYData);
278                     FREE(YSize);
279                     return types::Function::Error;
280                 }
281             }
282             else if (pDblAtol == NULL && bFuncF == false)
283             {
284                 pDblAtol = in[iPos]->getAs<types::Double>();
285                 if (pDblAtol->getSize() != pDblY0->getSize() && pDblAtol->isScalar() == false)
286                 {
287                     Scierror(267, _("%s: Arg %d and arg %d must have equal dimensions.\n"), "ode", pStrType ? 2 : 1, iPos + 1);
288                     DifferentialEquation::removeDifferentialEquationFunctions();
289                     FREE(pdYData);
290                     FREE(YSize);
291                     return types::Function::Error;
292                 }
293             }
294             else if (pDblNg == NULL && bFuncF == true && meth == 3)
295             {
296                 pDblNg = in[iPos]->getAs<types::Double>();
297             }
298             else if (pDblW == NULL && bFuncF == true && (bFuncG == true || meth != 3))
299             {
300                 if ((int)in.size() == iPos + 2)
301                 {
302                     if (in[iPos + 1]->isDouble() == false)
303                     {
304                         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 2);
305                         DifferentialEquation::removeDifferentialEquationFunctions();
306                         FREE(pdYData);
307                         FREE(YSize);
308                         return types::Function::Error;
309                     }
310
311                     pDblW = in[iPos]->getAs<types::Double>();
312                     pDblIw = in[iPos + 1]->getAs<types::Double>();
313                     iPos++;
314                 }
315                 else
316                 {
317                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", iPos + 2);
318                     DifferentialEquation::removeDifferentialEquationFunctions();
319                     FREE(pdYData);
320                     FREE(YSize);
321                     return types::Function::Error;
322                 }
323             }
324             else
325             {
326                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "ode", iPos + 1);
327                 DifferentialEquation::removeDifferentialEquationFunctions();
328                 FREE(pdYData);
329                 FREE(YSize);
330                 return types::Function::Error;
331             }
332         }
333         else if (in[iPos]->isCallable())
334         {
335             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
336             if (bFuncF == false)
337             {
338                 deFunctionsManager->setFFunction(pCall);
339                 bFuncF = true;
340             }
341             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
342             {
343                 deFunctionsManager->setJacFunction(pCall);
344                 bFuncJac = true;
345             }
346             else if (bFuncG == false && meth == 3)
347             {
348                 deFunctionsManager->setGFunction(pCall);
349                 bFuncG = true;
350             }
351             else
352             {
353                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
354                 DifferentialEquation::removeDifferentialEquationFunctions();
355                 FREE(pdYData);
356                 FREE(YSize);
357                 return types::Function::Error;
358             }
359         }
360         else if (in[iPos]->isString())
361         {
362             types::String* pStr = in[iPos]->getAs<types::String>();
363             bool bOK = false;
364
365             if (bFuncF == false)
366             {
367                 bOK = deFunctionsManager->setFFunction(pStr);
368                 bFuncF = true;
369             }
370             else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
371             {
372                 bOK = deFunctionsManager->setJacFunction(pStr);
373                 bFuncJac = true;
374             }
375             else if (bFuncG == false && meth == 3)
376             {
377                 bOK = deFunctionsManager->setGFunction(pStr);
378                 bFuncG = true;
379             }
380             else
381             {
382                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
383                 DifferentialEquation::removeDifferentialEquationFunctions();
384                 FREE(pdYData);
385                 FREE(YSize);
386                 return types::Function::Error;
387             }
388
389             if (bOK == false)
390             {
391                 char* pst = wide_string_to_UTF8(pStr->get(0));
392                 Scierror(50, _("%s: Subroutine not found: %s\n"), "ode", pst);
393                 FREE(pst);
394                 DifferentialEquation::removeDifferentialEquationFunctions();
395                 FREE(pdYData);
396                 FREE(YSize);
397                 return types::Function::Error;
398             }
399         }
400         else if (in[iPos]->isList())
401         {
402             types::List* pList = in[iPos]->getAs<types::List>();
403
404             if (pList->getSize() == 0)
405             {
406                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, "(string empty)");
407                 DifferentialEquation::removeDifferentialEquationFunctions();
408                 FREE(pdYData);
409                 FREE(YSize);
410                 return types::Function::Error;
411             }
412
413             if (bFuncF && (bFuncJac || pDblNg) && (bFuncG || meth != 3))
414             {
415                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "ode", iPos + 1);
416                 DifferentialEquation::removeDifferentialEquationFunctions();
417                 FREE(pdYData);
418                 FREE(YSize);
419                 return types::Function::Error;
420             }
421
422             if (pList->get(0)->isString())
423             {
424                 types::String* pStr = pList->get(0)->getAs<types::String>();
425                 bool bOK = false;
426
427                 if (bFuncF == false)
428                 {
429                     bFuncF = true;
430                     bOK = deFunctionsManager->setFFunction(pStr);
431                     sizeOfpdYData = *YSize;
432                 }
433                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
434                 {
435                     bFuncJac = true;
436                     bOK = deFunctionsManager->setJacFunction(pStr);
437                     if (sizeOfpdYData == 0)
438                     {
439                         sizeOfpdYData = *YSize;
440                     }
441                 }
442                 else if (bFuncG == false && meth == 3)
443                 {
444                     bFuncG = true;
445                     bOK = deFunctionsManager->setGFunction(pStr);
446                     if (sizeOfpdYData == 0)
447                     {
448                         sizeOfpdYData = *YSize;
449                     }
450                 }
451
452                 if (bOK == false)
453                 {
454                     char* pst = wide_string_to_UTF8(pStr->get(0));
455                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "ode", iPos + 1, pst);
456                     FREE(pst);
457                     DifferentialEquation::removeDifferentialEquationFunctions();
458                     FREE(pdYData);
459                     FREE(YSize);
460                     return types::Function::Error;
461                 }
462
463                 int* sizeTemp = YSize;
464                 int totalSize = sizeOfpdYData;
465
466                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
467                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
468
469                 std::vector<types::Double*> vpDbl;
470                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
471                 {
472                     if (pList->get(iter + 1)->isDouble() == false)
473                     {
474                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "ode", iPos + 1, iter + 1);
475                         DifferentialEquation::removeDifferentialEquationFunctions();
476                         FREE(pdYData);
477                         FREE(YSize);
478                         return types::Function::Error;
479                     }
480
481                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
482                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
483                     totalSize += YSize[sizeOfYSize + iter];
484                 }
485
486                 double* pdYDataTemp = pdYData;
487                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
488                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
489
490                 int position = sizeOfpdYData;
491                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
492                 {
493                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
494                     position += vpDbl[iter]->getSize();
495                 }
496                 vpDbl.clear();
497                 sizeOfpdYData = totalSize;
498                 sizeOfYSize += pList->getSize() - 1;
499                 FREE(pdYDataTemp);
500                 FREE(sizeTemp);
501             }
502             else if (pList->get(0)->isCallable())
503             {
504                 if (bFuncF == false)
505                 {
506                     bFuncF = true;
507                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
508                     for (int iter = 1; iter < pList->getSize(); iter++)
509                     {
510                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
511                     }
512                 }
513                 else if (bFuncJac == false && (pDblNg == NULL || meth != 3))
514                 {
515                     bFuncJac = true;
516                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
517                     for (int iter = 1; iter < pList->getSize(); iter++)
518                     {
519                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
520                     }
521                 }
522                 else if (bFuncG == false && meth == 3)
523                 {
524                     bFuncG = true;
525                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
526                     for (int iter = 1; iter < pList->getSize(); iter++)
527                     {
528                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
529                     }
530                 }
531             }
532             else
533             {
534                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "ode", iPos + 1);
535                 DifferentialEquation::removeDifferentialEquationFunctions();
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"), "ode", iPos + 1);
544             DifferentialEquation::removeDifferentialEquationFunctions();
545             FREE(pdYData);
546             FREE(YSize);
547             return types::Function::Error;
548         }
549     }
550
551     if (bFuncF == false)
552     {
553         int val = (meth == 3) ? 3 : 1;
554         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + val);
555         DifferentialEquation::removeDifferentialEquationFunctions();
556         FREE(pdYData);
557         FREE(YSize);
558         return types::Function::Error;
559     }
560     if (pDblNg == NULL && meth == 3)
561     {
562         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 2);
563         DifferentialEquation::removeDifferentialEquationFunctions();
564         FREE(pdYData);
565         FREE(YSize);
566         return types::Function::Error;
567     }
568     if (bFuncG == false && meth == 3)
569     {
570         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "ode", in.size() + 1);
571         DifferentialEquation::removeDifferentialEquationFunctions();
572         FREE(pdYData);
573         FREE(YSize);
574         return types::Function::Error;
575     }
576
577     // *** Initialization. ***
578
579     int itol        = 1;
580     int itask       = 1; // itask = 1 for normal computation of output values of y at t = tout.
581     int istate      = 1; // istate = integer flag (input and output).  set istate = 1.
582     int iopt        = 0; // iopt   = 0 to indicate no optional inputs used.
583     int jt          = 2; // jacobian type indicator.  set jt = 2.
584     int ml          = -1;
585     int mu          = -1;
586
587     // work tab
588     double* rwork   = NULL;
589     int* iwork      = NULL;
590     int rworkSize   = 0;
591     int iworkSize   = 0;
592
593     // contain ls0001, lsa001 and eh0001 structures
594     double* dStructTab  = NULL;
595     int* iStructTab     = NULL;
596     int dStructTabSize  = 0;
597     int iStructTabSize  = 0;
598
599     int rwSize  = 0; // rwSize = dStructTab + rworkSize
600     int iwSize  = 0; // iwSize = iStructTab + iworkSize
601
602     // structures used by lsoda and lsode
603     double* ls0001d = &(C2F(ls0001).tret);
604     int* ls0001i    = &(C2F(ls0001).illin);
605     double* lsa001d = &(C2F(lsa001).tsw);
606     int* lsa001i    = &(C2F(lsa001).insufr);
607     double* lsr001d = &(C2F(lsr001).rownr3[0]);
608     int* lsr001i    = &(C2F(lsr001).lg0);
609     int* eh0001i    = &(C2F(eh0001).mesflg);
610
611     // get %ODEOPTIONS
612     types::InternalType* pIT = symbol::Context::getInstance()->get(symbol::Symbol(L"%ODEOPTIONS"));
613     if (pIT != NULL && pIT->isDouble())
614     {
615         pDblOdeOptions = pIT->getAs<types::Double>();
616         if (pDblOdeOptions->getSize() == 12)
617         {
618             iopt    = 1;
619             itask   = (int)pDblOdeOptions->get(0);
620             jt      = (int)pDblOdeOptions->get(5);
621             ml      = (int)pDblOdeOptions->get(10);
622             mu      = (int)pDblOdeOptions->get(11);
623         }
624     }
625
626     if (iopt == 1 && (pDblOdeOptions->get(4) > pDblOdeOptions->get(3))) // hmin > hmax ?
627     {
628         Scierror(9999, _("%s: Wrong value of hmin and hmax: hmin = %d is greater than hmax = %d.\n"), "ode", pDblOdeOptions->get(4), pDblOdeOptions->get(3));
629         DifferentialEquation::removeDifferentialEquationFunctions();
630         FREE(pdYData);
631         FREE(YSize);
632         return types::Function::Error;
633     }
634
635     if (jt < 0 || jt > 5)
636     {
637         Scierror(9999, _("%s: Wrong value of Jacobian type: A number between %d and %d expected.\n"), "ode", 0, 5);
638         DifferentialEquation::removeDifferentialEquationFunctions();
639         FREE(pdYData);
640         FREE(YSize);
641         return types::Function::Error;
642     }
643
644     if (iopt == 0 && bFuncJac)
645     {
646         jt = 1;
647     }
648
649     if (bFuncJac && (jt == 2 || jt == 5) && getWarningMode())
650     {
651         sciprint(_("%s: Warning: Jacobian is given, but not used.\n"), "ode");
652     }
653
654     if (bFuncJac == false && (jt == 1 || jt == 4))
655     {
656         if (getWarningMode())
657         {
658             sciprint(_("%s: Warning: No Jacobian external given, but one is required by %ODEOPTIONS(6) value. Jacobian will be estimated.\n"), "ode");
659         }
660
661         jt = 2;
662     }
663
664     //compute itol and set the tolerances rtol and atol.
665     double* rtol = NULL;
666     double* atol = NULL;
667
668     if (pDblRtol)
669     {
670         if (pDblRtol->isScalar())
671         {
672             rtol = (double*)MALLOC(sizeof(double));
673             *rtol = pDblRtol->get(0);
674         }
675         else
676         {
677             rtol = pDblRtol->get();
678             itol += 2;
679         }
680     }
681     else
682     {
683         rtol = (double*)MALLOC(sizeof(double));
684         if (meth == 6 || meth == 7)
685         {
686             *rtol = 1.e-3;
687         }
688         else
689         {
690             *rtol = 1.0e-7;
691         }
692     }
693
694     if (pDblAtol)
695     {
696         if (pDblAtol->isScalar())
697         {
698             atol = (double*)MALLOC(sizeof(double));
699             *atol = pDblAtol->get(0);
700         }
701         else
702         {
703             atol = pDblAtol->get();
704             itol ++;
705         }
706     }
707     else
708     {
709         atol = (double*)MALLOC(sizeof(double));
710         if (meth == 6 || meth == 7)
711         {
712             *atol = 1.e-4;
713         }
714         else
715         {
716             *atol = 1.0e-9;
717         }
718     }
719
720     // Compute rwork, iwork size.
721     // Create them.
722     int nyh = (*YSize);
723     if (pDblW) // structure ls0001 have been restored.
724     {
725         nyh = C2F(ls0001).nyh;
726     }
727
728     int mxordn = 12;
729     int mxords = 5;
730
731     if (iopt == 1)
732     {
733         mxordn = (int)pDblOdeOptions->get(7);
734         mxords = (int)pDblOdeOptions->get(8);
735     }
736
737     if (mxordn > 12 || mxords > 5 || mxordn < 1 || mxords < 1)
738     {
739         if (getWarningMode())
740         {
741             sciprint(_("%s: Warning: Wrong value for maximun stiff/non-stiff order allowed :\nAt most %d for mxordn, %d for mxords and no null value for both expected.\nWrong value will be reduced to the default value.\n"), "ode", 12, 5);
742         }
743
744         mxordn = 12;
745         mxords = 5;
746     }
747
748     int maxord  = mxords;
749     int lmat    = 0;
750     int lrn     = 0;
751     int lrs     = 0;
752
753     switch (meth)
754     {
755         case 3: // lsodar (root)
756         {
757             //             lrn = 20 + nyh*(mxordn+1) + 3*neq + 3*ng,
758             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat + 3*ng,
759             //          where
760             //             nyh    = the initial value of neq,
761             //             mxordn = 12, unless a smaller value is given as an
762             //                      optional input,
763             //             mxords = 5, unless a smaller value is given as an
764             //                      optional input,
765             //             lmat   = length of matrix work space..
766             //             lmat   = neq**2 + 2              if jt = 1 or 2,
767             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
768
769             lrn = 3 * (int)pDblNg->get(0);
770             lrs = lrn;
771             dStructTabSize = 246 - 241;
772             iStructTabSize = 59 - 50;
773             jroot = (int*)MALLOC((int)pDblNg->get(0) * sizeof(int));
774         }
775         case 0: // lsoda
776         {
777             //             lrn = 20 + nyh*(mxordn+1) + 3*neq,
778             //             lrs = 20 + nyh*(mxords+1) + 3*neq + lmat,
779             //          where
780             //             nyh    = the initial value of neq,
781             //             mxordn = 12, unless a smaller value is given as an
782             //                      optional input,
783             //             mxords = 5, unless a smaller value is given as an
784             //                      optional input,
785             //             lmat   = length of matrix work space..
786             //             lmat   = neq**2 + 2              if jt = 1 or 2,
787             //             lmat   = (2*ml + mu + 1)*neq + 2 if jt = 4 or 5.
788
789             if (jt == 1 || jt == 2)
790             {
791                 lmat = (*YSize) * (*YSize) + 2;       //  if jt = 1 or 2,
792             }
793             else if (jt == 4 || jt == 5)
794             {
795                 //ml and mu = -1 in all cases
796                 lmat = (2 * ml + mu + 1) * (*YSize) + 2; //  if jt = 4 or 5.
797             }
798
799             lrn += 20 + nyh * (mxordn + 1) + 3 * (*YSize);
800             lrs += 20 + nyh * (mxords + 1) + 3 * (*YSize) + lmat;
801
802             rworkSize   = max(lrn, lrs);
803             iworkSize   = 20 + *YSize;
804
805             dStructTabSize += 241;
806             iStructTabSize += 50;
807
808             break;
809         }
810         case 1: // lsode (adams) (non stiff)
811         {
812             maxord = mxordn;
813         }
814         case 2: // lsode (stiff)
815         {
816             //          20 + nyh*(maxord + 1) + 3*neq + lmat
817             //        where
818             //          nyh    = the initial value of neq,
819             //          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
820             //                   smaller value is given as an optional input),
821             //          lmat   = 0             if miter = 0,
822             //          lmat   = neq**2 + 2    if miter is 1 or 2,
823             //          lmat   = neq + 2       if miter = 3, and
824             //          lmat   = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
825
826             if (jt == 1 || jt == 2)
827             {
828                 lmat = (*YSize) * (*YSize) + 2;
829             }
830             else if (jt == 3)
831             {
832                 lmat = (*YSize) + 2;
833             }
834             else if (jt == 4 || jt == 5)
835             {
836                 lmat = (2 * ml + mu + 1) * (*YSize) + 2;
837             }
838
839             rworkSize = 20 + nyh * (maxord + 1) + 3 * (*YSize) + lmat;
840             iworkSize = 20; // if jt = 0 or 3
841             if (jt == 1 || jt == 2 || jt == 4 || jt == 5) // iSize = 20 + neq  otherwise
842             {
843                 iworkSize += (*YSize);
844             }
845
846             dStructTabSize = 219;
847             iStructTabSize = 41;
848             break;
849         }
850         case 4: // lsdisc (discrete)
851         {
852             rworkSize = *YSize;
853             iworkSize = 1;
854             break;
855         }
856         case 5: // lsrgk (rk)
857         {
858             rworkSize = 9 * (*YSize);
859             iworkSize = 1;
860             break;
861         }
862         case 6: // rkf45 (rkf)
863         {
864             rworkSize = 3 + 8 * (*YSize);
865             iworkSize = 5;
866             break;
867         }
868         case 7: // rksimp (fix)
869         {
870             rworkSize = 3 + 8 * (*YSize);
871             iworkSize = 1;
872             break;
873         }
874     }
875
876     rwSize = rworkSize + dStructTabSize;
877     iwSize = iworkSize + iStructTabSize;
878
879     rwork = (double*)MALLOC(rworkSize * sizeof(double));
880     iwork = (int*)MALLOC(iworkSize * sizeof(int));
881
882     if (meth < 4)
883     {
884         if (pDblW && pDblIw)
885         {
886             if (pDblW->getSize() != rwSize || pDblIw->getSize() != iwSize)
887             {
888                 Scierror(9999, _("%s: Wrong size for w and iw: w = %d and iw = %d expected.\n"), "ode", rwSize, iwSize);
889                 DifferentialEquation::removeDifferentialEquationFunctions();
890                 FREE(pdYData);
891                 FREE(YSize);
892                 FREE(rwork);
893                 FREE(iwork);
894                 if (jroot)
895                 {
896                     FREE(jroot);
897                 }
898                 if (itol == 1 || itol == 3)
899                 {
900                     FREE(atol);
901                 }
902                 if (itol < 3)
903                 {
904                     FREE(rtol);
905                 }
906                 return types::Function::Error;
907             }
908
909             istate = 2; // 1 means this is the first call | 2  means this is not the first call
910
911             dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
912             iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
913
914             C2F(dcopy)(&rworkSize, pDblW->get(), &one, rwork, &one);
915             C2F(dcopy)(&dStructTabSize, pDblW->get() + rworkSize, &one, dStructTab, &one);
916
917             for (int i = 0; i < iworkSize; i++)
918             {
919                 iwork[i] = (int)pDblIw->get(i);
920             }
921
922             for (int i = 0; i < iStructTabSize; i++)
923             {
924                 iStructTab[i] = (int)pDblIw->get(i + iworkSize);
925             }
926         }
927         else
928         {
929             //if iopt = 0 default value are used without set it.
930             if (iopt == 1)
931             {
932                 if (itask >= 4)
933                 {
934                     rwork[0] = pDblOdeOptions->get(1); // tcrit
935                 }
936
937                 rwork[4] = pDblOdeOptions->get(2);      // h0 =0
938                 rwork[5] = pDblOdeOptions->get(3);      // hmax = %inf
939                 rwork[6] = pDblOdeOptions->get(4);      // hmin = 0
940                 iwork[0] = (int)pDblOdeOptions->get(10);// ml = -1
941                 iwork[1] = (int)pDblOdeOptions->get(11);// mu = -1
942                 iwork[4] = (int)pDblOdeOptions->get(9); // ixpr = 0
943                 iwork[5] = (int)pDblOdeOptions->get(6); // mxstep = 500
944                 iwork[6] = 10;  // mxhnil = 10 maximum number of messages printed per problem
945                 iwork[7] = (int)pDblOdeOptions->get(7); // mxordn = 12
946                 iwork[8] = (int)pDblOdeOptions->get(8); // mxords = 5
947             }
948         }
949     }
950
951     if (pDblOdeOptions && pDblOdeOptions->get(9) == 1)
952     {
953         sciprint(_("itask = %d\tmeth = %d\tjactyp = %d\tml = %d\tmu = %d\tiopt = %d\n"), itask, meth, jt, ml, mu, iopt);
954         sciprint(_("tcrit = %lf\th0 = %lf\thmax = %lf\thmin = %lf\n"), pDblOdeOptions->get(1), pDblOdeOptions->get(2), pDblOdeOptions->get(3), pDblOdeOptions->get(4));
955     }
956
957     if (meth < 4)
958     {
959         if (pDblW && pDblIw)
960         {
961             int dPos    = 0;
962             int iPos    = 0;
963             int dSize   = 219;
964
965             //restore ls0001 and eh0001 from w (rwork) and iw (iwork).
966             C2F(dcopy)(&dSize, dStructTab, &one, ls0001d, &one);
967             memcpy(ls0001i, iStructTab, 39 * sizeof(int));
968
969             dPos = dSize;
970             iPos = 39;
971
972             //restore lsa001 from w (rwork) and iw (iwork).
973             if (meth == 0 || meth == 3)
974             {
975                 dSize = 22;
976                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsa001d, &one);
977                 memcpy(lsa001i, &iStructTab[iPos], 9 * sizeof(int));
978
979                 dPos += dSize;
980                 iPos += 9;
981             }
982
983             //restore lsr001
984             if (meth == 3)
985             {
986                 dSize = 5;
987                 C2F(dcopy)(&dSize, dStructTab + dPos, &one, lsr001d, &one);
988                 memcpy(lsr001i, &iStructTab[iPos], 9 * sizeof(int));
989             }
990         }
991     }
992
993     // *** Perform operation. ***
994     int err = 0;
995     double t0 = pDblT0->get(0);
996     bool bOneStep = false;
997
998     if (itask == 5 || itask == 3 || itask == 2)
999     {
1000         bOneStep = true;
1001         if (getWarningMode() && pDblT->isScalar() == false)
1002         {
1003             sciprint(_("itask = %d: At most one value of t is allowed, the last element of t is used.\n"), itask);
1004         }
1005     }
1006
1007     if (bOneStep)
1008     {
1009         std::list<double*> pDblYOutList = std::list<double*>();
1010         std::list<double> pDblTOutList = std::list<double>();
1011         int iLastT = pDblT->getSize() - 1;
1012         double t = pDblT->get(iLastT);
1013         int iDir = t - t0 < 0 ? -1 : 1;
1014
1015         do
1016         {
1017             std::string strMeth;
1018             try
1019             {
1020                 switch (meth)
1021                 {
1022                     case 1: // lsode (adams)
1023                     case 2: // lsode (stiff)
1024                     {
1025                         strMeth = "lsode";
1026                         int jacType = 10 * meth + jt;
1027                         C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1028                         break;
1029                     }
1030                     case 3: // lsodar
1031                     {
1032                         strMeth = "lsodar";
1033                         int ng = (int)pDblNg->get(0);
1034                         C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1035                         break;
1036                     }
1037                     case 4: // lsdisc (discrete)
1038                     {
1039                         strMeth = "lsdisc";
1040                         C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1041                         break;
1042                     }
1043                     case 5: // lsrgk (rk)
1044                     {
1045                         strMeth = "lsrgk";
1046                         C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1047                         break;
1048                     }
1049                     case 6: // rkf45 (rkf)
1050                     {
1051                         strMeth = "rkf45";
1052                         C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1053                         break;
1054                     }
1055                     case 7: // rksimp (fix)
1056                     {
1057                         strMeth = "rksimp";
1058                         C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1059                         break;
1060                     }
1061                     default: // case 0: // lsoda
1062                     {
1063                         strMeth = "lsoda";
1064                         C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1065                     }
1066                 }
1067
1068                 // check error
1069                 err = checkOdeError(meth, istate);
1070                 if (err == 1)
1071                 {
1072                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1073                 }
1074             }
1075             catch (ScilabError &e)
1076             {
1077                 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1078                 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1079                 Scierror(999, pstrMsg);
1080                 FREE(pstrMsg);
1081                 err = 1;
1082             }
1083
1084             // FREE allocated data
1085             if (err == 1) // error case
1086             {
1087                 DifferentialEquation::removeDifferentialEquationFunctions();
1088                 FREE(pdYData);
1089                 FREE(YSize);
1090                 FREE(rwork);
1091                 FREE(iwork);
1092                 if (jroot)
1093                 {
1094                     FREE(jroot);
1095                 }
1096                 if (dStructTab)
1097                 {
1098                     FREE(dStructTab);
1099                 }
1100                 if (iStructTab)
1101                 {
1102                     FREE(iStructTab);
1103                 }
1104                 if (itol == 1 || itol == 3)
1105                 {
1106                     FREE(atol);
1107                 }
1108                 if (itol < 3)
1109                 {
1110                     FREE(rtol);
1111                 }
1112                 return types::Function::Error;
1113             }
1114
1115             pDblYOutList.push_back(pdYData);
1116             pDblTOutList.push_back(t0);
1117
1118             if (err == 2) // warning case
1119             {
1120                 if (getWarningMode())
1121                 {
1122                     sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1123                 }
1124                 break;
1125             }
1126
1127             if (meth == 3 && istate == 3 && getWarningMode())
1128             {
1129                 // istate == 3  means the integration was successful, and one or more
1130                 //              roots were found before satisfying the stop condition
1131                 //              specified by itask.  see jroot.
1132
1133                 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1134                 for (int k = 0; k < pDblNg->get(0); k++)
1135                 {
1136                     sciprint("\t%d", jroot[k]);
1137                 }
1138                 sciprint("\n");
1139                 break;
1140             }
1141         }
1142         while ((t0 - t) * iDir < 0);
1143
1144         int iSizeList = (int)pDblYOutList.size();
1145         /*
1146                 if(pPolyY0) // pPolyY0 is scalar
1147                 {
1148                     // result format = [2 x iSizeList]
1149                     // first lime is the time of result stored in second line.
1150                     int size = 2 * iSizeList;
1151                     int* iRanks = (int*)MALLOC(size * sizeof(int));
1152
1153                     for(int i = 0; i < iSizeList; i++)
1154                     {
1155                         iRanks[i * 2] = 1; // time rank
1156                         iRanks[i * 2 + 1] = pPolyY0->getMaxRank(); // result rank
1157                     }
1158
1159                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 2, iSizeList, iRanks);
1160                     pDblYOut = new types::Double(pDblY0->getRows(), 1);
1161                     pDblTOut = new types::Double(1, 1);
1162
1163                     for(int i = 0; i < iSizeList; i++)
1164                     {
1165                         // pDblY0 contain pPolyY0 (ie: pPolyY0 = 2 + x => pDblY0 = [2 1])
1166                         for(int j = 0; j < pDblY0->getRows(); j++)
1167                         {
1168                             pDblYOut->set(j, pDblYOutList.front()[j]);
1169                         }
1170                         pDblTOut->set(0, pDblTOutList.front());
1171
1172                         pPolyYOut->setCoef(i * 2, pDblTOut);
1173                         pPolyYOut->setCoef(i * 2 + 1, pDblYOut);
1174
1175                         pDblYOutList.pop_front();
1176                         pDblTOutList.pop_front();
1177                     }
1178                 }
1179                 else
1180                 {
1181         */
1182         pDblYOut = new types::Double(pDblY0->getRows() + 1, iSizeList);
1183         for (int i = 0; i < iSizeList; i++)
1184         {
1185             pDblYOut->set(i * (*YSize + 1), pDblTOutList.front());
1186             for (int j = 0; j < *YSize; j++)
1187             {
1188                 pDblYOut->set(i * (*YSize + 1) + (j + 1), pDblYOutList.front()[j]);
1189             }
1190             pDblYOutList.pop_front();
1191             pDblTOutList.pop_front();
1192         }
1193         //        }
1194     }
1195     else
1196     {
1197         // Create result
1198         /*        if(pPolyY0)
1199                 {
1200                     int size = pDblT->getSize();
1201                     int* iRanks = (int*)MALLOC(size * sizeof(int));
1202                     for(int i = 0; i < size; i++)
1203                     {
1204                         iRanks[i] = pPolyY0->getMaxRank();
1205                     }
1206
1207                     pPolyYOut = new types::Polynom(pPolyY0->getVariableName(), 1, pDblT->getSize(), iRanks);
1208
1209                     pDblYOut = new types::Double(pDblY0->getRows(), 1);
1210                 }
1211                 else
1212                 {
1213         */
1214         pDblYOut = new types::Double(pDblY0->getRows(), pDblT->getSize());
1215         //        }
1216         bool bBreak = false;
1217         for (int i = 0; i < pDblT->getSize(); i++)
1218         {
1219             double t = pDblT->get(i);
1220             std::string strMeth;
1221
1222             if (itask >= 4 && t > rwork[0]) // rwork[0] => tcrit
1223             {
1224                 t = rwork[0];
1225                 bBreak = true;
1226             }
1227
1228             try
1229             {
1230                 switch (meth)
1231                 {
1232                     case 1: // lsode (adams)
1233                     case 2: // lsode (stiff)
1234                     {
1235                         strMeth = "lsode";
1236                         int jacType = 10 * meth + jt;
1237                         C2F(lsode)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jacType);
1238                         break;
1239                     }
1240                     case 3: // lsodar
1241                     {
1242                         strMeth = "lsodar";
1243                         int ng = (int)pDblNg->get(0);
1244                         C2F(lsodar)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt, ode_g, &ng, jroot);
1245                         break;
1246                     }
1247                     case 4: // lsdisc (discrete)
1248                     {
1249                         strMeth = "lsdisc";
1250                         C2F(lsdisc)(ode_f, YSize, pdYData, &t0, &t, rwork, &rworkSize, &istate);
1251                         break;
1252                     }
1253                     case 5: // lsrgk (rk)
1254                     {
1255                         strMeth = "lsrgk";
1256                         C2F(lsrgk)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1257                         break;
1258                     }
1259                     case 6: // rkf45 (rkf)
1260                     {
1261                         strMeth = "rkf45";
1262                         C2F(rkf45)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1263                         break;
1264                     }
1265                     case 7: // rksimp (fix)
1266                     {
1267                         strMeth = "rksimp";
1268                         C2F(rksimp)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &meth);
1269                         break;
1270                     }
1271                     default: // case 0: // lsoda
1272                     {
1273                         strMeth = "lsoda";
1274                         C2F(lsoda)(ode_f, YSize, pdYData, &t0, &t, &itol, rtol, atol, &itask, &istate, &iopt, rwork, &rworkSize, iwork, &iworkSize, bFuncJac ? ode_jac : NULL, &jt);
1275                     }
1276                 }
1277
1278                 // check error
1279                 err = checkOdeError(meth, istate);
1280                 if (err == 1)
1281                 {
1282                     Scierror(999, _("%s: %s exit with state %d.\n"), "ode", strMeth.c_str(), istate);
1283                 }
1284             }
1285             catch (ScilabError &e)
1286             {
1287                 char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1288                 sciprint(_("%s: exception caught in '%s' subroutine.\n"), "ode", strMeth.c_str());
1289                 Scierror(999, pstrMsg);
1290                 FREE(pstrMsg);
1291                 err = 1;
1292             }
1293
1294             // FREE allocated data
1295             if (err == 1) // error case
1296             {
1297                 DifferentialEquation::removeDifferentialEquationFunctions();
1298                 FREE(pdYData);
1299                 FREE(YSize);
1300                 FREE(rwork);
1301                 FREE(iwork);
1302                 if (jroot)
1303                 {
1304                     FREE(jroot);
1305                 }
1306                 if (dStructTab)
1307                 {
1308                     FREE(dStructTab);
1309                 }
1310                 if (iStructTab)
1311                 {
1312                     FREE(iStructTab);
1313                 }
1314                 if (itol == 1 || itol == 3)
1315                 {
1316                     FREE(atol);
1317                 }
1318                 if (itol < 3)
1319                 {
1320                     FREE(rtol);
1321                 }
1322                 return types::Function::Error;
1323             }
1324
1325             /*            if(pPolyYOut)
1326                         {
1327                             for(int j = 0; j < pDblY0->getRows(); j++)
1328                             {
1329                                 pDblYOut->set(j, pdYData[j]);
1330                             }
1331                             pPolyYOut->setCoef(i, pDblYOut);
1332                         }
1333                         else
1334                         {
1335             */
1336             for (int j = 0; j < *YSize; j++)
1337             {
1338                 pDblYOut->set(i * (*YSize) + j, pdYData[j]);
1339             }
1340             //            }
1341
1342             if (err == 2) // warning case
1343             {
1344                 if (getWarningMode())
1345                 {
1346                     sciprint(_("Integration was stoped at t = %lf.\n"), t0);
1347                 }
1348
1349                 types::Double* pDblYOutTemp = pDblYOut;
1350                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1351                 for (int k = 0; k < pDblYOut->getSize(); k++)
1352                 {
1353                     pDblYOut->set(k, pDblYOutTemp->get(k));
1354                 }
1355                 delete pDblYOutTemp;
1356
1357                 break;
1358             }
1359
1360             if (meth == 3 && istate == 3 && getWarningMode())
1361             {
1362                 // istate == 3  means the integration was successful, and one or more
1363                 //              roots were found before satisfying the stop condition
1364                 //              specified by itask.  see jroot.
1365
1366                 sciprint(_("%s: Warning: At t = %lf, y is a root, jroot = "), "ode", t0);
1367                 for (int k = 0; k < pDblNg->get(0); k++)
1368                 {
1369                     sciprint("\t%d", jroot[k]);
1370                 }
1371                 sciprint("\n");
1372
1373                 types::Double* pDblYOutTemp = pDblYOut;
1374                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1375                 for (int k = 0; k < pDblYOut->getSize(); k++)
1376                 {
1377                     pDblYOut->set(k, pDblYOutTemp->get(k));
1378                 }
1379                 delete pDblYOutTemp;
1380                 break;
1381             }
1382
1383             if (itask >= 4 && bBreak)
1384             {
1385                 types::Double* pDblYOutTemp = pDblYOut;
1386                 pDblYOut = new types::Double(pDblYOutTemp->getRows(), i + 1);
1387                 for (int k = 0; k < pDblYOut->getSize(); k++)
1388                 {
1389                     pDblYOut->set(k, pDblYOutTemp->get(k));
1390                 }
1391                 delete pDblYOutTemp;
1392                 break;
1393             }
1394         }
1395     }
1396
1397     if (meth < 4)
1398     {
1399         if (_iRetCount > 2) //save ls0001 and eh0001 following w and iw.
1400         {
1401             int dPos    = 0;
1402             int iPos    = 0;
1403             int dSize   = 219;
1404
1405             if (dStructTab == NULL)
1406             {
1407                 dStructTab = (double*)MALLOC(dStructTabSize * sizeof(double));
1408                 iStructTab = (int*)MALLOC(iStructTabSize * sizeof(int));
1409             }
1410
1411             C2F(dcopy)(&dSize, ls0001d, &one, dStructTab, &one);
1412             memcpy(iStructTab, ls0001i, 39 * sizeof(int));
1413
1414             dPos = dSize;
1415             iPos = 39;
1416
1417             //save lsa001
1418             if (meth == 0 || meth == 3)
1419             {
1420                 dSize = 22;
1421                 C2F(dcopy)(&dSize, lsa001d, &one, dStructTab + dPos, &one);
1422                 memcpy(&iStructTab[iPos], lsa001i, 9 * sizeof(int));
1423
1424                 dPos += dSize;
1425                 iPos += 9;
1426             }
1427
1428             //save lsr001
1429             if (meth == 3)
1430             {
1431                 dSize = 5;
1432                 C2F(dcopy)(&dSize, lsr001d, &one, dStructTab + dPos, &one);
1433                 memcpy(&iStructTab[iPos], lsr001i, 9 * sizeof(int));
1434
1435                 iPos += 9;
1436             }
1437
1438             memcpy(&iStructTab[iPos], eh0001i, 2 * sizeof(int));
1439         }
1440     }
1441
1442     // *** Return result in Scilab. ***
1443     /*
1444         if(pPolyYOut)
1445         {
1446             out.push_back(pPolyYOut); // y
1447         }
1448         else
1449         {
1450     */
1451     out.push_back(pDblYOut); // y
1452     //    }
1453
1454     if (meth == 3 && _iRetCount >= 2)
1455     {
1456         int sizeOfRd = 1;
1457         int k = 0;
1458
1459         for (int i = 0; i < pDblNg->get(0); i++)
1460         {
1461             if (jroot[i])
1462             {
1463                 sizeOfRd++;
1464             }
1465         }
1466
1467         types::Double* pDblRd = new types::Double(1, sizeOfRd);
1468         //rd: The first entry contains the stopping time.
1469         pDblRd->set(0, C2F(lsr001).tlast);
1470         for (int i = 0; i < pDblNg->get(0); i++)
1471         {
1472             if (jroot[i])
1473             {
1474                 k++;
1475                 pDblRd->set(k, (double)i + 1);
1476             }
1477         }
1478         out.push_back(pDblRd); // rd
1479     }
1480
1481     if ((meth < 3 && _iRetCount == 3) || (meth == 3 && _iRetCount == 4))
1482     {
1483         types::Double* pDblWOut = new types::Double(1, rwSize);
1484         C2F(dcopy)(&rworkSize, rwork, &one, pDblWOut->get(), &one);
1485         C2F(dcopy)(&dStructTabSize, dStructTab, &one, pDblWOut->get() + rworkSize, &one);
1486
1487         types::Double* pDblIwOut = new types::Double(1, iwSize);
1488         for (int i = 0; i < iworkSize; i++)
1489         {
1490             pDblIwOut->set(i, (double)iwork[i]);
1491         }
1492
1493         for (int i = 0; i < iStructTabSize; i++)
1494         {
1495             pDblIwOut->set(iworkSize + i, (double)iStructTab[i]);
1496         }
1497
1498         out.push_back(pDblWOut);  // w
1499         out.push_back(pDblIwOut); // iw
1500     }
1501
1502     // *** FREE. ***
1503     if (itol == 1 || itol == 3) // atol is scalar
1504     {
1505         FREE(atol);
1506     }
1507
1508     if (itol < 3) // rtol is scalar
1509     {
1510         FREE(rtol);
1511     }
1512
1513     FREE(pdYData);
1514     FREE(YSize);
1515     FREE(rwork);
1516     FREE(iwork);
1517
1518     if (jroot)
1519     {
1520         FREE(jroot);
1521     }
1522
1523     if (dStructTab)
1524     {
1525         FREE(dStructTab);
1526         FREE(iStructTab);
1527     }
1528
1529     DifferentialEquation::removeDifferentialEquationFunctions();
1530
1531     return types::Function::OK;
1532 }
1533