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