differential_equations plugged.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_dassl.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 "checkodeerror.h"
30 #include "xerhlt.h"
31 }
32
33 /*--------------------------------------------------------------------------*/
34 types::Function::ReturnValue sci_dassl(types::typed_list &in, int _iRetCount, types::typed_list &out)
35 {
36     // input args
37     types::Double* pDblX0   = NULL;
38     types::Double* pDblT0   = NULL;
39     types::Double* pDblT    = NULL;
40     types::Double* pDblRtol = NULL;
41     types::Double* pDblAtol = NULL;
42     types::Double* pDblHd   = NULL;
43
44     // x0 = [y0, ydot0]
45     double* pdYData         = NULL; // contain y0 following by all args data in list case.
46     double* pdYdotData      = NULL;
47     int sizeOfpdYData       = 0;
48
49     int iPos    = 0;
50     int one     = 1;
51
52     int info[15]    = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
53
54     double tstop    = 0;
55     double maxstep  = 0;
56     double stepin   = 0;
57     int mu          = 0;
58     int ml          = 0;
59
60     int sizeOfYSize = 1;
61     int* YSize      = NULL;    // YSize(1) = size of y0,
62                                // YSize(n) = size of Args(n) in list case.
63
64     // Indicate if the function is given.
65     bool bFuncF     = false;
66     bool bFuncJac   = false;
67
68     // Indicate if info list is given.
69     bool bListInfo  = false;
70
71 // *** check the minimal number of input args. ***
72     if(in.size() < 4 || in.size() > 9)
73     {
74         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d to %d expected.\n"), L"dassl", 4, 9);
75         return types::Function::Error;
76     }
77
78 // *** check number of output args ***
79     if(_iRetCount > 2)
80     {
81         ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d to %d expected.\n"), L"dassl", 1, 2);
82         return types::Function::Error;
83     }
84
85 // *** check type of input args and get it. ***
86     // x0 = [y0, yd0]
87     if(in[iPos]->isDouble() == false)
88     {
89         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"dassl", iPos+1);
90         return types::Function::Error;
91     }
92
93     pDblX0 = in[iPos]->getAs<types::Double>();
94
95     if(pDblX0->isComplex())
96     {
97         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A real matrix expected.\n"), L"dassl", iPos+1);
98         return types::Function::Error;
99     }
100
101     if(pDblX0->getCols() > 2)
102     {
103         ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A real matrix with %d to %d colomn(s) expected.\n"), L"dassl", iPos+1, 1, 2);
104         return types::Function::Error;
105     }
106
107     if(pDblX0->getCols() == 1)
108     {
109         info[10] = 1;
110     }
111
112     // t0
113     iPos++;
114     if(in[iPos]->isDouble() == false)
115     {
116         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A scalar expected.\n"), L"dassl", iPos+1);
117         return types::Function::Error;
118     }
119
120     pDblT0 = in[iPos]->getAs<types::Double>();
121
122     if(pDblT0->isScalar() == false)
123     {
124         ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A scalar expected.\n"), L"dassl", iPos+1);
125         return types::Function::Error;
126     }
127
128     // t
129     iPos++;
130     if(in[iPos]->isDouble() == false)
131     {
132         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"dassl", iPos+1);
133         return types::Function::Error;
134     }
135
136     pDblT = in[iPos]->getAs<types::Double>();
137
138     if(pDblT->isComplex())
139     {
140         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A real matrix expected.\n"), L"dassl", iPos+1);
141         return types::Function::Error;
142     }
143
144     // get next inputs
145     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"dassl");
146     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
147
148     YSize = (int*)malloc(sizeOfYSize * sizeof(int));
149     *YSize = pDblX0->getRows();
150
151     pdYData = (double*)malloc(*YSize * sizeof(double));
152     pdYdotData = (double*)malloc(*YSize * sizeof(double));
153
154     C2F(dcopy)(YSize, pDblX0->get(), &one, pdYData, &one);
155     if(pDblX0->getCols() == 2)
156     {
157         C2F(dcopy)(YSize, pDblX0->get() + *YSize, &one, pdYdotData, &one);
158     }
159     else
160     {
161         memset(pdYdotData, 0x00, *YSize);
162     }
163
164     deFunctionsManager->setOdeYRows(pDblX0->getRows());
165
166     for(iPos++; iPos < in.size(); iPos++)
167     {
168         if(in[iPos]->isDouble())
169         {
170             if(pDblAtol == NULL && bFuncF == false)
171             {
172                 pDblAtol = in[iPos]->getAs<types::Double>();
173                 if(pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
174                 {
175                     ScierrorW(267,_W("%ls: Wrong size for input argument #%d : A scalar or a matrix of size %d expected.\n"), L"dassl", iPos+1, pDblX0->getRows());
176                     DifferentialEquation::removeDifferentialEquationFunctions();
177                     free(pdYdotData);free(pdYData);free(YSize);
178                     return types::Function::Error;
179                 }
180             }
181             else if(pDblRtol == NULL && bFuncF == false)
182             {
183                 pDblRtol = in[iPos]->getAs<types::Double>();
184                 if(pDblAtol->getSize() != pDblRtol->getSize())
185                 {
186                     ScierrorW(267,_W("%ls: Wrong size for input argument #%d : Atol and Rtol must have the same size.\n"), L"dassl", iPos+1, pDblX0->getRows());
187                     DifferentialEquation::removeDifferentialEquationFunctions();
188                     free(pdYdotData);free(pdYData);free(YSize);
189                     return types::Function::Error;
190                 }
191             }
192             else if(pDblHd == NULL && bFuncF == true)
193             {
194                 pDblHd = in[iPos]->getAs<types::Double>();
195                 if(in.size() != iPos+1)
196                 {
197                     ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"dassl", iPos+1);
198                     DifferentialEquation::removeDifferentialEquationFunctions();
199                     free(pdYdotData);free(pdYData);free(YSize);
200                     return types::Function::Error;
201                 }
202             }
203             else
204             {
205                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A function expected.\n"), L"dassl", iPos+1);
206                 DifferentialEquation::removeDifferentialEquationFunctions();
207                 free(pdYdotData);free(pdYData);free(YSize);
208                 return types::Function::Error;
209             }
210         }
211         else if(in[iPos]->isCallable())
212         {
213             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
214             if(bFuncF == false)
215             {
216                 deFunctionsManager->setFFunction(pCall);
217                 bFuncF = true;
218             }
219             else if(bFuncJac == false)
220             {
221                 deFunctionsManager->setJacFunction(pCall);
222                 bFuncJac = true;
223             }
224             else
225             {
226                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix or a list expected.\n"), L"dassl", iPos+1);
227                 DifferentialEquation::removeDifferentialEquationFunctions();
228                 free(pdYdotData);free(pdYData);free(YSize);
229                 return types::Function::Error;
230             }
231         }
232         else if(in[iPos]->isString())
233         {
234             types::String* pStr = in[iPos]->getAs<types::String>();
235             bool bOK = false;
236
237             if(bFuncF == false)
238             {
239                 bOK = deFunctionsManager->setFFunction(pStr);
240                 bFuncF = true;
241             }
242             else if(bFuncJac == false)
243             {
244                 bOK = deFunctionsManager->setJacFunction(pStr);
245                 bFuncJac = true;
246             }
247             else
248             {
249                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix or a list expected.\n"), L"dassl", iPos+1);
250                 DifferentialEquation::removeDifferentialEquationFunctions();
251                 free(pdYdotData);free(pdYData);free(YSize);
252                 return types::Function::Error;
253             }
254
255             if(bOK == false)
256             {
257                 ScierrorW(50,_W("%ls: Subroutine not found: %ls\n"), L"dassl", pStr->get(0));
258                 DifferentialEquation::removeDifferentialEquationFunctions();
259                 free(pdYdotData);free(pdYData);free(YSize);
260                 return types::Function::Error;
261             }
262         }
263         else if(in[iPos]->isList())
264         {
265             types::List* pList = in[iPos]->getAs<types::List>();
266
267             if(pList->getSize() == 0)
268             {
269                 ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"dassl", iPos+1, L"(string empty)");
270                 DifferentialEquation::removeDifferentialEquationFunctions();
271                 free(pdYdotData);free(pdYData);free(YSize);
272                 return types::Function::Error;
273             }
274
275             if(bFuncF && bListInfo)
276             {
277                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"dassl", iPos+1);
278                 DifferentialEquation::removeDifferentialEquationFunctions();
279                 free(pdYdotData);free(pdYData);free(YSize);
280                 return types::Function::Error;
281             }
282
283             if(pList->get(0)->isString())
284             {
285                 types::String* pStr = pList->get(0)->getAs<types::String>();
286                 bool bOK = false;
287
288                 if(bFuncF == false)
289                 {
290                     bFuncF = true;
291                     bOK = deFunctionsManager->setFFunction(pStr);
292                     sizeOfpdYData = *YSize;
293                 }
294                 else if(bFuncJac == false)
295                 {
296                     bFuncJac = true;
297                     bOK = deFunctionsManager->setJacFunction(pStr);
298                     if(sizeOfpdYData == 0)
299                     {
300                         sizeOfpdYData = *YSize;
301                     }
302                 }
303
304                 if(bOK == false)
305                 {
306                     ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"dassl", iPos+1, pStr->get(0));
307                     DifferentialEquation::removeDifferentialEquationFunctions();
308                     free(pdYdotData);free(pdYData);free(YSize);
309                     return types::Function::Error;
310                 }
311
312                 int* sizeTemp = YSize;
313                 int totalSize = sizeOfpdYData;
314
315                 YSize = (int*)malloc((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
316                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
317
318                 std::vector<types::Double*> vpDbl;
319                 for(int iter = 0; iter < pList->getSize() - 1; iter++)
320                 {
321                     if(pList->get(iter + 1)->isDouble() == false)
322                     {
323                         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : Argument %d in the list must be a matrix.\n"), L"dassl", iPos+1, iter+1);
324                         DifferentialEquation::removeDifferentialEquationFunctions();
325                         free(pdYdotData);free(pdYData);free(YSize);
326                         return types::Function::Error;
327                     }
328
329                     vpDbl.push_back(pList->get(iter+1)->getAs<types::Double>());
330                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
331                     totalSize += YSize[sizeOfYSize + iter];
332                 }
333
334                 double* pdYDataTemp = pdYData;
335                 pdYData = (double*)malloc(totalSize * sizeof(double));
336                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
337
338                 int position = sizeOfpdYData;
339                 for(int iter = 0; iter < pList->getSize()-1; iter++)
340                 {
341                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
342                     position += vpDbl[iter]->getSize();
343                 }
344                 vpDbl.clear();
345                 sizeOfpdYData = totalSize;
346                 sizeOfYSize += pList->getSize() - 1;
347                 free(pdYDataTemp);
348                 free(sizeTemp);
349             }
350             else if(pList->get(0)->isCallable())
351             {
352                 if(bFuncF == false)
353                 {
354                     bFuncF = true;
355                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
356                     for(int iter = 1; iter < pList->getSize(); iter++)
357                     {
358                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
359                     }
360                 }
361                 else if(bFuncJac == false)
362                 {
363                     bFuncJac = true;
364                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
365                     for(int iter = 1; iter < pList->getSize(); iter++)
366                     {
367                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
368                     }
369                 }
370             }
371             else if(pList->get(0)->isDouble() && bFuncF == true)
372             {
373                 if(pList->getSize() != 7)
374                 {
375                     ScierrorW(267,_W("%ls: Wrong size for input argument #%d : A list of size %d expected.\n"), L"dassl", iPos+1, 7);
376                     DifferentialEquation::removeDifferentialEquationFunctions();
377                     free(pdYdotData);free(pdYData);free(YSize);
378                     return types::Function::Error;
379                 }
380
381                 for(int i = 0; i < 7; i++) // info = list([],0,[],[],[],0,0)
382                 {
383                     if(pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 6)))
384                     {
385                         if(i == 1 || i == 5 || i == 6)
386                         {
387                             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : Element %d in the info list must be a scalar.\n"), L"dassl", iPos+1, i);
388                         }
389                         else
390                         {
391                             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : Element %d in the info list must be a matrix.\n"), L"dassl", iPos+1, i);
392                         }
393                         DifferentialEquation::removeDifferentialEquationFunctions();
394                         free(pdYdotData);free(pdYData);free(YSize);
395                         return types::Function::Error;  
396                     }
397                 }
398
399                 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
400                 if(pDblTemp->getSize() != 0)
401                 {
402                     info[3] = 1;
403                     tstop = pDblTemp->get(0);
404                 }
405
406                 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
407
408                 pDblTemp = pList->get(2)->getAs<types::Double>();
409                 if(pDblTemp->getSize() == 2)
410                 {
411                     info[5] = 1;
412                     ml = (int)pDblTemp->get(0);
413                     mu = (int)pDblTemp->get(1);
414                     deFunctionsManager->setMl(ml);
415                     deFunctionsManager->setMu(mu);
416                 }
417                 else if(pDblTemp->getSize() != 0)
418                 {
419                     ScierrorW(267,_W("%ls: Wrong size for input argument #%d : Argument %d in te list must be of size %d.\n"), L"dassl", iPos+1, 3, 2);
420                     DifferentialEquation::removeDifferentialEquationFunctions();
421                     free(pdYdotData);free(pdYData);free(YSize);
422                     return types::Function::Error;
423                 }
424
425                 pDblTemp = pList->get(3)->getAs<types::Double>();
426                 if(pDblTemp->getSize() != 0)
427                 {
428                     info[6] = 1;
429                     maxstep = pDblTemp->get(0);
430                 }
431
432                 pDblTemp = pList->get(4)->getAs<types::Double>();
433                 if(pDblTemp->getSize() != 0)
434                 {
435                     info[7] = 1;
436                     stepin = pDblTemp->get(0);
437                 }
438
439                 info[9]  = (int)pList->get(5)->getAs<types::Double>()->get(0);
440                 if(pList->get(6)->getAs<types::Double>()->get(0) == 1)
441                 {
442                     info[10] = 1;
443                 }
444
445                 bListInfo = true;
446             }
447             else
448             {
449                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : The first argument in the list must be a string, a function or a matrix in case of argument info.\n"), L"dassl", iPos+1);
450                 DifferentialEquation::removeDifferentialEquationFunctions();
451                 free(pdYdotData);free(pdYData);free(YSize);
452                 return types::Function::Error;
453             }
454         }
455         else
456         {
457             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix or a function expected.\n"), L"dassl", iPos+1);
458             DifferentialEquation::removeDifferentialEquationFunctions();
459             free(pdYdotData);free(pdYData);free(YSize);
460             return types::Function::Error;
461         }
462     }
463
464     if(bFuncF == false)
465     {
466         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"dassl", in.size() + 1);
467         DifferentialEquation::removeDifferentialEquationFunctions();
468         free(pdYdotData);free(pdYData);free(YSize);
469         return types::Function::Error;
470     }
471
472     if(bFuncJac == true)
473     {
474         info[4] = 1;
475     }
476
477 // *** Initialization. ***
478     double t0   = pDblT0->get(0);
479     double rpar = 0;
480     int ipar    = 0;
481     int idid    = 0;
482     int maxord  = 5;
483
484     //compute itol and set the tolerances rtol and atol.
485     double* rtol = NULL;
486     double* atol = NULL;
487
488     if(pDblAtol)
489     {
490         if(pDblAtol->isScalar())
491         {
492             atol  = (double*)malloc(sizeof(double));
493             *atol = pDblAtol->get(0);
494         }
495         else
496         {
497             atol    = pDblAtol->get();
498             info[1] = 1;
499         }
500     }
501     else
502     {
503         atol  = (double*)malloc(sizeof(double));
504         *atol = 1.e-7;
505     }
506
507     if(pDblRtol)
508     {
509         if(pDblRtol->isScalar())
510         {
511             rtol  = (double*)malloc(sizeof(double));
512             *rtol = pDblRtol->get(0);
513         }
514         else
515         {
516             rtol = pDblRtol->get();
517         }
518     }
519     else // if rtol is not given atol will be used as a scalar.
520     {
521         if(pDblAtol && pDblAtol->isScalar() == false)// info[1] == 1
522         {
523             double dblSrc = 1.e-9;
524             int iSize = pDblAtol->getSize();
525             int iOne = 1;
526             int iZero = 0;
527
528             rtol = (double*)malloc(iSize * sizeof(double));
529             C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
530         }
531         else
532         {
533             rtol    = (double*)malloc(sizeof(double));
534             *rtol   = 1.e-9;
535         }
536     }
537
538     // Compute rwork, iwork size.
539     // Create them.
540     int iworksize   = 20 + pDblX0->getRows();
541     int rworksize   = 0;
542     int* iwork      = NULL;
543     double* rwork   = NULL;
544
545     if(info[5] == 0)
546     {
547         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + pDblX0->getRows() * pDblX0->getRows();
548     }
549     else if(info[4] == 1)
550     {
551         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows();
552     }
553     else if(info[4] == 0)
554     {
555         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
556     }
557
558     iwork = (int*)malloc(iworksize * sizeof(int));
559     rwork = (double*)malloc(rworksize * sizeof(double));
560
561     if(pDblHd != NULL)
562     {
563         if(iworksize + rworksize != pDblHd->getSize())
564         {
565             ScierrorW(77, _W("%ls: Wrong size for input argument(s) %d: %d expected.\n"), L"dassl", in.size(), iworksize + rworksize);
566             DifferentialEquation::removeDifferentialEquationFunctions();
567             free(pdYdotData);free(pdYData);free(YSize);
568             free(iwork);free(rwork);
569             if(pDblAtol == NULL || pDblAtol->isScalar()) free(atol);
570             if(pDblRtol == NULL || pDblRtol->isScalar()) free(rtol);
571             return types::Function::Error;
572         }
573
574         C2F(dcopy)(&rworksize, pDblHd->get(), &one, rwork, &one);
575
576         for(int i = 0; i < iworksize; i++)
577         {
578             iwork[i] = (int)pDblHd->get(rworksize + i);
579         }
580
581         info[0] = 1;
582     }
583
584     if(info[3] == 1)
585     {
586         rwork[0] = tstop;
587     }
588
589     if(info[6] == 1)
590     {
591         rwork[1] = maxstep;
592     }
593
594     if(info[7] == 1)
595     {
596         rwork[2] = stepin;
597     }
598
599     if(info[5] == 1)
600     {
601         iwork[0] = ml;
602         iwork[1] = mu;
603     }
604
605 // *** Perform operation. ***
606     std::list<types::Double*> lpDblOut;
607     int size = pDblX0->getRows();
608     int rowsOut = 1 + pDblX0->getRows() * 2;
609
610     for(int i = 0; i < pDblT->getSize(); i++)
611     {
612         types::Double* pDblOut = new types::Double(rowsOut, 1);
613         lpDblOut.push_back(pDblOut);
614
615         double t = pDblT->get(i);
616         int pos  = 0;
617         pDblOut->set(pos, t);
618
619         if(t == t0)
620         {
621             pos++;
622             C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
623             pos += pDblX0->getRows();
624             C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
625
626             continue;
627         }
628
629         C2F(dassl)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t, info, rtol, atol, &idid, rwork, &rworksize, iwork, &iworksize, &rpar, &ipar, dassl_jac);
630
631         int iret = checkDasslError(idid);
632
633         if(iret == 1) // error
634         {
635             ScierrorW(999, _W("%ls: dassl return with state %d.\n"), L"dassl", idid);
636             lpDblOut.clear();
637             DifferentialEquation::removeDifferentialEquationFunctions();
638             free(pdYdotData);free(pdYData);free(YSize);
639             free(iwork);free(rwork);
640             if(pDblAtol == NULL || pDblAtol->isScalar()) free(atol);
641             if(pDblRtol == NULL || pDblRtol->isScalar()) free(rtol);
642             return types::Function::Error;
643         }
644
645         pos++;
646         C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
647         pos += size;
648         C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
649
650         if(iret == 2) // warning
651         {
652             break;
653         }
654
655         // iret == 0
656         if(idid == 1)
657         {
658             pDblOut->set(0, t0);
659             i--;
660         }
661         else if(idid == -2)
662         {
663             t0 = t;
664             i--;
665         }
666         else
667         {
668             t0 = t;
669         }
670
671         info[0] = 1;
672     }
673
674 // *** Return result in Scilab. ***
675     types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
676
677     int sizeOfList = (int)lpDblOut.size();
678     for(int i = 0; i < sizeOfList; i++)
679     {
680         int pos = i * rowsOut;
681         C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &one, pDblOut->get() + pos, &one);
682         lpDblOut.pop_front();
683     }
684     out.push_back(pDblOut);
685
686     if(_iRetCount == 2)
687     {
688         types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
689         C2F(dcopy)(&rworksize, rwork, &one, pDblHdOut->get(), &one);
690
691         for(int i = 0; i < iworksize; i++)
692         {
693             pDblHdOut->set(rworksize + i, (double)iwork[i]);
694         }
695
696         out.push_back(pDblHdOut);
697     }
698
699 // *** free. ***
700     if(pDblAtol == NULL || pDblAtol->isScalar())
701     {
702         free(atol);
703     }
704
705     if(pDblRtol == NULL || pDblRtol->isScalar())
706     {
707         free(rtol);
708     }
709
710     free(pdYData);
711     free(pdYdotData);
712     free(YSize);
713     free(rwork);
714     free(iwork);
715
716     DifferentialEquation::removeDifferentialEquationFunctions();
717
718     return types::Function::OK;
719 }
720