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