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