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