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