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