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