Management of errors in differential equations functions.
[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);
179                     FREE(pdYData);
180                     FREE(YSize);
181                     return types::Function::Error;
182                 }
183             }
184             else if (pDblRtol == NULL && bFuncF == false)
185             {
186                 pDblRtol = in[iPos]->getAs<types::Double>();
187                 if (pDblAtol->getSize() != pDblRtol->getSize())
188                 {
189                     Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "dassl", iPos + 1, pDblX0->getRows());
190                     DifferentialEquation::removeDifferentialEquationFunctions();
191                     FREE(pdYdotData);
192                     FREE(pdYData);
193                     FREE(YSize);
194                     return types::Function::Error;
195                 }
196             }
197             else if (pDblHd == NULL && bFuncF == true)
198             {
199                 pDblHd = in[iPos]->getAs<types::Double>();
200                 if (in.size() != iPos + 1)
201                 {
202                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", iPos + 1);
203                     DifferentialEquation::removeDifferentialEquationFunctions();
204                     FREE(pdYdotData);
205                     FREE(pdYData);
206                     FREE(YSize);
207                     return types::Function::Error;
208                 }
209             }
210             else
211             {
212                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "dassl", iPos + 1);
213                 DifferentialEquation::removeDifferentialEquationFunctions();
214                 FREE(pdYdotData);
215                 FREE(pdYData);
216                 FREE(YSize);
217                 return types::Function::Error;
218             }
219         }
220         else if (in[iPos]->isCallable())
221         {
222             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
223             if (bFuncF == false)
224             {
225                 deFunctionsManager->setFFunction(pCall);
226                 bFuncF = true;
227             }
228             else if (bFuncJac == false)
229             {
230                 deFunctionsManager->setJacFunction(pCall);
231                 bFuncJac = true;
232             }
233             else
234             {
235                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
236                 DifferentialEquation::removeDifferentialEquationFunctions();
237                 FREE(pdYdotData);
238                 FREE(pdYData);
239                 FREE(YSize);
240                 return types::Function::Error;
241             }
242         }
243         else if (in[iPos]->isString())
244         {
245             types::String* pStr = in[iPos]->getAs<types::String>();
246             bool bOK = false;
247
248             if (bFuncF == false)
249             {
250                 bOK = deFunctionsManager->setFFunction(pStr);
251                 bFuncF = true;
252             }
253             else if (bFuncJac == false)
254             {
255                 bOK = deFunctionsManager->setJacFunction(pStr);
256                 bFuncJac = true;
257             }
258             else
259             {
260                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "dassl", iPos + 1);
261                 DifferentialEquation::removeDifferentialEquationFunctions();
262                 FREE(pdYdotData);
263                 FREE(pdYData);
264                 FREE(YSize);
265                 return types::Function::Error;
266             }
267
268             if (bOK == false)
269             {
270                 char* pst = wide_string_to_UTF8(pStr->get(0));
271                 Scierror(50, _("%s: Subroutine not found: %s\n"), "dassl", pst);
272                 FREE(pst);
273                 DifferentialEquation::removeDifferentialEquationFunctions();
274                 FREE(pdYdotData);
275                 FREE(pdYData);
276                 FREE(YSize);
277                 return types::Function::Error;
278             }
279         }
280         else if (in[iPos]->isList())
281         {
282             types::List* pList = in[iPos]->getAs<types::List>();
283
284             if (pList->getSize() == 0)
285             {
286                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, "(string empty)");
287                 DifferentialEquation::removeDifferentialEquationFunctions();
288                 FREE(pdYdotData);
289                 FREE(pdYData);
290                 FREE(YSize);
291                 return types::Function::Error;
292             }
293
294             if (bFuncF && bListInfo)
295             {
296                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "dassl", iPos + 1);
297                 DifferentialEquation::removeDifferentialEquationFunctions();
298                 FREE(pdYdotData);
299                 FREE(pdYData);
300                 FREE(YSize);
301                 return types::Function::Error;
302             }
303
304             if (pList->get(0)->isString())
305             {
306                 types::String* pStr = pList->get(0)->getAs<types::String>();
307                 bool bOK = false;
308
309                 if (bFuncF == false)
310                 {
311                     bFuncF = true;
312                     bOK = deFunctionsManager->setFFunction(pStr);
313                     sizeOfpdYData = *YSize;
314                 }
315                 else if (bFuncJac == false)
316                 {
317                     bFuncJac = true;
318                     bOK = deFunctionsManager->setJacFunction(pStr);
319                     if (sizeOfpdYData == 0)
320                     {
321                         sizeOfpdYData = *YSize;
322                     }
323                 }
324
325                 if (bOK == false)
326                 {
327                     char* pst = wide_string_to_UTF8(pStr->get(0));
328                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "dassl", iPos + 1, pst);
329                     FREE(pst);
330                     DifferentialEquation::removeDifferentialEquationFunctions();
331                     FREE(pdYdotData);
332                     FREE(pdYData);
333                     FREE(YSize);
334                     return types::Function::Error;
335                 }
336
337                 int* sizeTemp = YSize;
338                 int totalSize = sizeOfpdYData;
339
340                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
341                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
342
343                 std::vector<types::Double*> vpDbl;
344                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
345                 {
346                     if (pList->get(iter + 1)->isDouble() == false)
347                     {
348                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "dassl", iPos + 1, iter + 1);
349                         DifferentialEquation::removeDifferentialEquationFunctions();
350                         FREE(pdYdotData);
351                         FREE(pdYData);
352                         FREE(YSize);
353                         return types::Function::Error;
354                     }
355
356                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
357                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
358                     totalSize += YSize[sizeOfYSize + iter];
359                 }
360
361                 double* pdYDataTemp = pdYData;
362                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
363                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &one, pdYData, &one);
364
365                 int position = sizeOfpdYData;
366                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
367                 {
368                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &one, &pdYData[position], &one);
369                     position += vpDbl[iter]->getSize();
370                 }
371                 vpDbl.clear();
372                 sizeOfpdYData = totalSize;
373                 sizeOfYSize += pList->getSize() - 1;
374                 FREE(pdYDataTemp);
375                 FREE(sizeTemp);
376             }
377             else if (pList->get(0)->isCallable())
378             {
379                 if (bFuncF == false)
380                 {
381                     bFuncF = true;
382                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
383                     for (int iter = 1; iter < pList->getSize(); iter++)
384                     {
385                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
386                     }
387                 }
388                 else if (bFuncJac == false)
389                 {
390                     bFuncJac = true;
391                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
392                     for (int iter = 1; iter < pList->getSize(); iter++)
393                     {
394                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
395                     }
396                 }
397             }
398             else if (pList->get(0)->isDouble() && bFuncF == true)
399             {
400                 if (pList->getSize() != 7)
401                 {
402                     Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "dassl", iPos + 1, 7);
403                     DifferentialEquation::removeDifferentialEquationFunctions();
404                     FREE(pdYdotData);
405                     FREE(pdYData);
406                     FREE(YSize);
407                     return types::Function::Error;
408                 }
409
410                 for (int i = 0; i < 7; i++) // info = list([],0,[],[],[],0,0)
411                 {
412                     if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 6)))
413                     {
414                         if (i == 1 || i == 5 || i == 6)
415                         {
416                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "dassl", iPos + 1, i);
417                         }
418                         else
419                         {
420                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "dassl", iPos + 1, i);
421                         }
422                         DifferentialEquation::removeDifferentialEquationFunctions();
423                         FREE(pdYdotData);
424                         FREE(pdYData);
425                         FREE(YSize);
426                         return types::Function::Error;
427                     }
428                 }
429
430                 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
431                 if (pDblTemp->getSize() != 0)
432                 {
433                     info[3] = 1;
434                     tstop = pDblTemp->get(0);
435                 }
436
437                 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
438
439                 pDblTemp = pList->get(2)->getAs<types::Double>();
440                 if (pDblTemp->getSize() == 2)
441                 {
442                     info[5] = 1;
443                     ml = (int)pDblTemp->get(0);
444                     mu = (int)pDblTemp->get(1);
445                     deFunctionsManager->setMl(ml);
446                     deFunctionsManager->setMu(mu);
447                 }
448                 else if (pDblTemp->getSize() != 0)
449                 {
450                     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);
451                     DifferentialEquation::removeDifferentialEquationFunctions();
452                     FREE(pdYdotData);
453                     FREE(pdYData);
454                     FREE(YSize);
455                     return types::Function::Error;
456                 }
457
458                 pDblTemp = pList->get(3)->getAs<types::Double>();
459                 if (pDblTemp->getSize() != 0)
460                 {
461                     info[6] = 1;
462                     maxstep = pDblTemp->get(0);
463                 }
464
465                 pDblTemp = pList->get(4)->getAs<types::Double>();
466                 if (pDblTemp->getSize() != 0)
467                 {
468                     info[7] = 1;
469                     stepin = pDblTemp->get(0);
470                 }
471
472                 info[9]  = (int)pList->get(5)->getAs<types::Double>()->get(0);
473                 if (pList->get(6)->getAs<types::Double>()->get(0) == 1)
474                 {
475                     info[10] = 1;
476                 }
477
478                 bListInfo = true;
479             }
480             else
481             {
482                 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);
483                 DifferentialEquation::removeDifferentialEquationFunctions();
484                 FREE(pdYdotData);
485                 FREE(pdYData);
486                 FREE(YSize);
487                 return types::Function::Error;
488             }
489         }
490         else
491         {
492             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "dassl", iPos + 1);
493             DifferentialEquation::removeDifferentialEquationFunctions();
494             FREE(pdYdotData);
495             FREE(pdYData);
496             FREE(YSize);
497             return types::Function::Error;
498         }
499     }
500
501     if (bFuncF == false)
502     {
503         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "dassl", in.size() + 1);
504         DifferentialEquation::removeDifferentialEquationFunctions();
505         FREE(pdYdotData);
506         FREE(pdYData);
507         FREE(YSize);
508         return types::Function::Error;
509     }
510
511     if (bFuncJac == true)
512     {
513         info[4] = 1;
514     }
515
516     // *** Initialization. ***
517     double t0   = pDblT0->get(0);
518     double rpar = 0;
519     int ipar    = 0;
520     int idid    = 0;
521     int maxord  = 5;
522
523     //compute itol and set the tolerances rtol and atol.
524     double* rtol = NULL;
525     double* atol = NULL;
526
527     if (pDblAtol)
528     {
529         if (pDblAtol->isScalar())
530         {
531             atol  = (double*)MALLOC(sizeof(double));
532             *atol = pDblAtol->get(0);
533         }
534         else
535         {
536             atol    = pDblAtol->get();
537             info[1] = 1;
538         }
539     }
540     else
541     {
542         atol  = (double*)MALLOC(sizeof(double));
543         *atol = 1.e-7;
544     }
545
546     if (pDblRtol)
547     {
548         if (pDblRtol->isScalar())
549         {
550             rtol  = (double*)MALLOC(sizeof(double));
551             *rtol = pDblRtol->get(0);
552         }
553         else
554         {
555             rtol = pDblRtol->get();
556         }
557     }
558     else // if rtol is not given atol will be used as a scalar.
559     {
560         if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
561         {
562             double dblSrc = 1.e-9;
563             int iSize = pDblAtol->getSize();
564             int iOne = 1;
565             int iZero = 0;
566
567             rtol = (double*)MALLOC(iSize * sizeof(double));
568             C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
569         }
570         else
571         {
572             rtol    = (double*)MALLOC(sizeof(double));
573             *rtol   = 1.e-9;
574         }
575     }
576
577     // Compute rwork, iwork size.
578     // Create them.
579     int iworksize   = 20 + pDblX0->getRows();
580     int rworksize   = 0;
581     int* iwork      = NULL;
582     double* rwork   = NULL;
583
584     if (info[5] == 0)
585     {
586         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + pDblX0->getRows() * pDblX0->getRows();
587     }
588     else if (info[4] == 1)
589     {
590         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows();
591     }
592     else if (info[4] == 0)
593     {
594         rworksize = 40 + (maxord + 4) * pDblX0->getRows() + (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
595     }
596
597     iwork = (int*)MALLOC(iworksize * sizeof(int));
598     rwork = (double*)MALLOC(rworksize * sizeof(double));
599
600     if (pDblHd != NULL)
601     {
602         if (iworksize + rworksize != pDblHd->getSize())
603         {
604             Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "dassl", in.size(), iworksize + rworksize);
605             DifferentialEquation::removeDifferentialEquationFunctions();
606             FREE(pdYdotData);
607             FREE(pdYData);
608             FREE(YSize);
609             FREE(iwork);
610             FREE(rwork);
611             if (pDblAtol == NULL || pDblAtol->isScalar())
612             {
613                 FREE(atol);
614             }
615             if (pDblRtol == NULL || pDblRtol->isScalar())
616             {
617                 FREE(rtol);
618             }
619             return types::Function::Error;
620         }
621
622         C2F(dcopy)(&rworksize, pDblHd->get(), &one, rwork, &one);
623
624         for (int i = 0; i < iworksize; i++)
625         {
626             iwork[i] = (int)pDblHd->get(rworksize + i);
627         }
628
629         info[0] = 1;
630     }
631
632     if (info[3] == 1)
633     {
634         rwork[0] = tstop;
635     }
636
637     if (info[6] == 1)
638     {
639         rwork[1] = maxstep;
640     }
641
642     if (info[7] == 1)
643     {
644         rwork[2] = stepin;
645     }
646
647     if (info[5] == 1)
648     {
649         iwork[0] = ml;
650         iwork[1] = mu;
651     }
652
653     // *** Perform operation. ***
654     std::list<types::Double*> lpDblOut;
655     int size = pDblX0->getRows();
656     int rowsOut = 1 + pDblX0->getRows() * 2;
657     int iret = 0;
658
659     for (int i = 0; i < pDblT->getSize(); i++)
660     {
661         types::Double* pDblOut = new types::Double(rowsOut, 1);
662         lpDblOut.push_back(pDblOut);
663
664         double t = pDblT->get(i);
665         int pos  = 0;
666         pDblOut->set(pos, t);
667
668         if (t == t0)
669         {
670             pos++;
671             C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
672             pos += pDblX0->getRows();
673             C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
674
675             continue;
676         }
677
678         try
679         {
680             C2F(dassl)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t, info, rtol, atol, &idid, rwork, &rworksize, iwork, &iworksize, &rpar, &ipar, dassl_jac);
681
682             iret = checkDasslError(idid);
683             if (iret == 1) // error
684             {
685                 Scierror(999, _("%s: dassl return with state %d.\n"), "dassl", idid);
686             }
687         }
688         catch (ScilabError &e)
689         {
690             char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
691             sciprint(_("%s: exception caught in '%s' subroutine.\n"), "dassl", "dassl");
692             Scierror(999, pstrMsg);
693             // set iret to 1 for FREE allocated memory
694             iret = 1;
695         }
696
697         if (iret == 1) // error
698         {
699             lpDblOut.clear();
700             DifferentialEquation::removeDifferentialEquationFunctions();
701             FREE(pdYdotData);
702             FREE(pdYData);
703             FREE(YSize);
704             FREE(iwork);
705             FREE(rwork);
706             if (pDblAtol == NULL || pDblAtol->isScalar())
707             {
708                 FREE(atol);
709             }
710             if (pDblRtol == NULL || pDblRtol->isScalar())
711             {
712                 FREE(rtol);
713             }
714             return types::Function::Error;
715         }
716
717         pos++;
718         C2F(dcopy)(&size, pdYData, &one, pDblOut->get() + pos, &one);
719         pos += size;
720         C2F(dcopy)(&size, pdYdotData, &one, pDblOut->get() + pos, &one);
721
722         if (iret == 2) // warning
723         {
724             break;
725         }
726
727         // iret == 0
728         if (idid == 1)
729         {
730             pDblOut->set(0, t0);
731             i--;
732         }
733         else if (idid == -2)
734         {
735             t0 = t;
736             i--;
737         }
738         else
739         {
740             t0 = t;
741         }
742
743         info[0] = 1;
744     }
745
746     // *** Return result in Scilab. ***
747     types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
748
749     int sizeOfList = (int)lpDblOut.size();
750     for (int i = 0; i < sizeOfList; i++)
751     {
752         int pos = i * rowsOut;
753         C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &one, pDblOut->get() + pos, &one);
754         lpDblOut.pop_front();
755     }
756     out.push_back(pDblOut);
757
758     if (_iRetCount == 2)
759     {
760         types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
761         C2F(dcopy)(&rworksize, rwork, &one, pDblHdOut->get(), &one);
762
763         for (int i = 0; i < iworksize; i++)
764         {
765             pDblHdOut->set(rworksize + i, (double)iwork[i]);
766         }
767
768         out.push_back(pDblHdOut);
769     }
770
771     // *** FREE. ***
772     if (pDblAtol == NULL || pDblAtol->isScalar())
773     {
774         FREE(atol);
775     }
776
777     if (pDblRtol == NULL || pDblRtol->isScalar())
778     {
779         FREE(rtol);
780     }
781
782     FREE(pdYData);
783     FREE(pdYdotData);
784     FREE(YSize);
785     FREE(rwork);
786     FREE(iwork);
787
788     DifferentialEquation::removeDifferentialEquationFunctions();
789
790     return types::Function::OK;
791 }
792