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