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