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