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