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