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