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