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