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