Management of errors in differential equations functions.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_daskr.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
4 * Copyright (C) 2013 - Scilab Enterprises - Cedric DELAMARRE
5 *
6 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which
8 * you should have received as part of this distribution.  The terms
9 * are also available at
10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13 /*--------------------------------------------------------------------------*/
14
15 #include "differential_equations_gw.hxx"
16 #include "function.hxx"
17 #include "double.hxx"
18 #include "string.hxx"
19 #include "list.hxx"
20 #include "callable.hxx"
21 #include "differentialequationfunctions.hxx"
22 #include "runvisitor.hxx"
23
24 extern "C"
25 {
26 #include "MALLOC.h"
27 #include "localization.h"
28 #include "Scierror.h"
29 #include "scifunctions.h"
30 #include "elem_common.h"
31 #include "checkodeerror.h"
32 }
33
34 /*--------------------------------------------------------------------------*/
35 //  This code was inspired by sci_f_daskr.f from the same folder.
36 //  daskr() has more optional parameters than daskr(), so the gateway
37 //  has to adapt to the new options
38 /*--------------------------------------------------------------------------*/
39 //[y0,nvs[,hotdata]]=daskr(y0,t0,t1[,atol[,rtol]],res[,jac],nh,h[,info[,psol][,pjac]][,hotdata])
40 types::Function::ReturnValue sci_daskr(types::typed_list &in, int _iRetCount, types::typed_list &out)
41 {
42     // input args
43     types::Double* pDblX0   = NULL;
44     types::Double* pDblT0   = NULL;
45     types::Double* pDblT    = NULL;
46     types::Double* pDblRtol = NULL;
47     types::Double* pDblAtol = NULL;
48     types::Double* pDblHd   = NULL;
49     types::Double* pDblNg   = NULL;
50     types::Double* pDblE7   = NULL; // 7th element of list info
51     types::Double* pDblE12  = NULL; // 12th element of list info
52
53
54     // x0 = [y0, ydot0]
55     double* pdYData         = NULL; // contain y0 following by all args data in list case.
56     double* pdYdotData      = NULL;
57     int sizeOfpdYData       = 0;
58
59     int sizeOfYSize = 1;
60     int* YSize      = NULL;    // YSize(1) = size of y0,
61     // YSize(n) = size of Args(n) in list case.
62     int iPos    = 0;
63     int iOne    = 1;
64
65     int info[20]    = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
66
67     double tstop    = 0;
68     double maxstep  = 0;
69     double stepin   = 0;
70     double epli     = 0;
71     double dTemp    = 0;
72     double stptol   = 0;
73     double epinit   = 0;
74
75     int ng      = 0;
76     int mu      = 0;
77     int ml      = 0;
78     int maxl    = 0;
79     int kmp     = 0;
80     int nrmax   = 0;
81     int mxnit   = 0;
82     int mxnj    = 0;
83     int mxnh    = 0;
84     int lsoff   = 0;
85     int LID     = 0;
86     int LENIWP  = 0;
87     int LENWP   = 0;
88
89     // Indicate if the function is given.
90     bool bFuncF     = false;
91     bool bFuncJac   = false;
92     bool bFuncG     = false;
93     bool bFuncPsol  = false;
94     bool bFuncPjac  = false;
95
96     // Indicate if info list is given.
97     bool bListInfo  = false;
98
99     // *** check the minimal number of input args. ***
100     if ((int)in.size() < 6 || (int)in.size() > 13)
101     {
102         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "daskr", 6, 13);
103         return types::Function::Error;
104     }
105
106     // *** check number of output args ***
107     if (_iRetCount != 3 && _iRetCount != 2)
108     {
109         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "daskr", 2, 3);
110         return types::Function::Error;
111     }
112
113     // *** check type of input args and get it. ***
114     // x0 = [y0, yd0]
115     if (in[iPos]->isDouble() == false)
116     {
117         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
118         return types::Function::Error;
119     }
120
121     pDblX0 = in[iPos]->getAs<types::Double>();
122
123     if (pDblX0->isComplex())
124     {
125         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
126         return types::Function::Error;
127     }
128
129     if (pDblX0->getCols() > 2)
130     {
131         Scierror(999, _("%s: Wrong size for input argument #%d: A real matrix with %d to %d colomn(s) expected.\n"), "daskr", iPos + 1, 1, 2);
132         return types::Function::Error;
133     }
134
135     if (pDblX0->getCols() == 1)
136     {
137         info[10] = 1;
138     }
139
140     YSize = (int*)MALLOC(sizeOfYSize * sizeof(int));
141     *YSize = pDblX0->getRows();
142
143     pdYData = (double*)MALLOC(*YSize * sizeof(double));
144     pdYdotData = (double*)MALLOC(*YSize * sizeof(double));
145
146     C2F(dcopy)(YSize, pDblX0->get(), &iOne, pdYData, &iOne);
147     if (pDblX0->getCols() == 2)
148     {
149         C2F(dcopy)(YSize, pDblX0->get() + *YSize, &iOne, pdYdotData, &iOne);
150     }
151     else
152     {
153         memset(pdYdotData, 0x00, *YSize);
154     }
155
156     // t0
157     iPos++;
158     if (in[iPos]->isDouble() == false)
159     {
160         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
161         return types::Function::Error;
162     }
163
164     pDblT0 = in[iPos]->getAs<types::Double>();
165
166     if (pDblT0->isScalar() == false)
167     {
168         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
169         return types::Function::Error;
170     }
171
172     // t
173     iPos++;
174     if (in[iPos]->isDouble() == false)
175     {
176         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
177         return types::Function::Error;
178     }
179
180     pDblT = in[iPos]->getAs<types::Double>();
181
182     if (pDblT->isComplex())
183     {
184         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "daskr", iPos + 1);
185         return types::Function::Error;
186     }
187
188     // get next inputs
189     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"daskr");
190     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
191
192     deFunctionsManager->setOdeYRows(pDblX0->getRows());
193
194     for (iPos++; iPos < in.size(); iPos++)
195     {
196         if (in[iPos]->isDouble())
197         {
198             if (pDblAtol == NULL && bFuncF == false)
199             {
200                 pDblAtol = in[iPos]->getAs<types::Double>();
201                 if (pDblAtol->getSize() != pDblX0->getRows() && pDblAtol->isScalar() == false)
202                 {
203                     Scierror(267, _("%s: Wrong size for input argument #%d: A scalar or a matrix of size %d expected.\n"), "daskr", iPos + 1, pDblX0->getRows());
204                     DifferentialEquation::removeDifferentialEquationFunctions();
205                     FREE(pdYdotData);
206                     FREE(pdYData);
207                     FREE(YSize);
208                     return types::Function::Error;
209                 }
210             }
211             else if (pDblRtol == NULL && bFuncF == false)
212             {
213                 pDblRtol = in[iPos]->getAs<types::Double>();
214                 if (pDblAtol->getSize() != pDblRtol->getSize())
215                 {
216                     Scierror(267, _("%s: Wrong size for input argument #%d: Atol and Rtol must have the same size.\n"), "daskr", iPos + 1, pDblX0->getRows());
217                     DifferentialEquation::removeDifferentialEquationFunctions();
218                     FREE(pdYdotData);
219                     FREE(pdYData);
220                     FREE(YSize);
221                     return types::Function::Error;
222                 }
223             }
224             else if (pDblNg == NULL && bFuncF == true)
225             {
226                 pDblNg = in[iPos]->getAs<types::Double>();
227                 if (pDblNg->isScalar() == false)
228                 {
229                     Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "daskr", iPos + 1);
230                     DifferentialEquation::removeDifferentialEquationFunctions();
231                     FREE(pdYdotData);
232                     FREE(pdYData);
233                     FREE(YSize);
234                     return types::Function::Error;
235                 }
236                 ng = (int)pDblNg->get(0);
237             }
238             else if (pDblHd == NULL && bFuncF == true)
239             {
240                 pDblHd = in[iPos]->getAs<types::Double>();
241                 if ((int)in.size() != iPos + 1)
242                 {
243                     Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", iPos + 1);
244                     DifferentialEquation::removeDifferentialEquationFunctions();
245                     FREE(pdYdotData);
246                     FREE(pdYData);
247                     FREE(YSize);
248                     return types::Function::Error;
249                 }
250             }
251             else
252             {
253                 Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "daskr", iPos + 1);
254                 DifferentialEquation::removeDifferentialEquationFunctions();
255                 FREE(pdYdotData);
256                 FREE(pdYData);
257                 FREE(YSize);
258                 return types::Function::Error;
259             }
260         }
261         else if (in[iPos]->isCallable())
262         {
263             types::Callable* pCall = in[iPos]->getAs<types::Callable>();
264             if (bFuncF == false)
265             {
266                 deFunctionsManager->setFFunction(pCall);
267                 bFuncF = true;
268             }
269             else if (bFuncJac == false && pDblNg == NULL)
270             {
271                 deFunctionsManager->setJacFunction(pCall);
272                 bFuncJac = true;
273             }
274             else if (bFuncG == false && pDblNg)
275             {
276                 deFunctionsManager->setGFunction(pCall);
277                 bFuncG = true;
278             }
279             else if (bFuncG && bFuncPsol == false)
280             {
281                 deFunctionsManager->setPsolFunction(pCall);
282                 bFuncPsol = true;
283             }
284             else if (bFuncPsol && bFuncPjac == false)
285             {
286                 deFunctionsManager->setPjacFunction(pCall);
287                 bFuncPjac = true;
288             }
289             else
290             {
291                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
292                 DifferentialEquation::removeDifferentialEquationFunctions();
293                 FREE(pdYdotData);
294                 FREE(pdYData);
295                 FREE(YSize);
296                 return types::Function::Error;
297             }
298         }
299         else if (in[iPos]->isString())
300         {
301             types::String* pStr = in[iPos]->getAs<types::String>();
302             bool bOK = false;
303
304             if (bFuncF == false)
305             {
306                 bOK = deFunctionsManager->setFFunction(pStr);
307                 bFuncF = true;
308             }
309             else if (bFuncJac == false && pDblNg == NULL)
310             {
311                 bOK = deFunctionsManager->setJacFunction(pStr);
312                 bFuncJac = true;
313             }
314             else if (bFuncG == false && pDblNg)
315             {
316                 bOK = deFunctionsManager->setGFunction(pStr);
317                 bFuncG = true;
318             }
319             else if (bFuncG && bFuncPsol == false)
320             {
321                 bOK = deFunctionsManager->setPsolFunction(pStr);
322                 bFuncPsol = true;
323             }
324             else if (bFuncPsol && bFuncPjac == false)
325             {
326                 bOK = deFunctionsManager->setPjacFunction(pStr);
327                 bFuncPjac = true;
328             }
329             else
330             {
331                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a list expected.\n"), "daskr", iPos + 1);
332                 DifferentialEquation::removeDifferentialEquationFunctions();
333                 FREE(pdYdotData);
334                 FREE(pdYData);
335                 FREE(YSize);
336                 return types::Function::Error;
337             }
338
339             if (bOK == false)
340             {
341                 char* pst = wide_string_to_UTF8(pStr->get(0));
342                 Scierror(50, _("%s: Subroutine not found: %s\n"), "daskr", pst);
343                 FREE(pst);
344                 DifferentialEquation::removeDifferentialEquationFunctions();
345                 FREE(pdYdotData);
346                 FREE(pdYData);
347                 FREE(YSize);
348                 return types::Function::Error;
349             }
350         }
351         else if (in[iPos]->isList())
352         {
353             types::List* pList = in[iPos]->getAs<types::List>();
354
355             if (pList->getSize() == 0)
356             {
357                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, "(string empty)");
358                 DifferentialEquation::removeDifferentialEquationFunctions();
359                 FREE(pdYdotData);
360                 FREE(pdYData);
361                 FREE(YSize);
362                 return types::Function::Error;
363             }
364
365             if (bFuncF && bListInfo)
366             {
367                 Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "daskr", iPos + 1);
368                 DifferentialEquation::removeDifferentialEquationFunctions();
369                 FREE(pdYdotData);
370                 FREE(pdYData);
371                 FREE(YSize);
372                 return types::Function::Error;
373             }
374
375             // functions can be passed in a list
376             if (pList->get(0)->isString())
377             {
378                 // list(string,...  function case
379                 types::String* pStr = pList->get(0)->getAs<types::String>();
380                 bool bOK = false;
381                 sizeOfpdYData = *YSize;
382
383                 if (bFuncF == false)
384                 {
385                     bFuncF = true;
386                     bOK = deFunctionsManager->setFFunction(pStr);
387                 }
388                 else if (bFuncJac == false && pDblNg == NULL)
389                 {
390                     bFuncJac = true;
391                     bOK = deFunctionsManager->setJacFunction(pStr);
392                 }
393                 else if (bFuncG == false && pDblNg)
394                 {
395                     bFuncG = true;
396                     bOK = deFunctionsManager->setGFunction(pStr);
397                 }
398                 else if (bFuncG && bFuncPsol == false)
399                 {
400                     bOK = deFunctionsManager->setPsolFunction(pStr);
401                     bFuncPsol = true;
402                 }
403                 else if (bFuncPsol && bFuncPjac == false)
404                 {
405                     bOK = deFunctionsManager->setPjacFunction(pStr);
406                     bFuncPjac = true;
407                 }
408
409                 if (bOK == false)
410                 {
411                     char* pst = wide_string_to_UTF8(pStr->get(0));
412                     Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "daskr", iPos + 1, pst);
413                     FREE(pst);
414                     DifferentialEquation::removeDifferentialEquationFunctions();
415                     FREE(pdYdotData);
416                     FREE(pdYData);
417                     FREE(YSize);
418                     return types::Function::Error;
419                 }
420
421                 int* sizeTemp = YSize;
422                 int totalSize = sizeOfpdYData;
423
424                 YSize = (int*)MALLOC((sizeOfYSize + pList->getSize() - 1) * sizeof(int));
425                 memcpy(YSize, sizeTemp, sizeOfYSize * sizeof(int));
426
427                 std::vector<types::Double*> vpDbl;
428                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
429                 {
430                     if (pList->get(iter + 1)->isDouble() == false)
431                     {
432                         Scierror(999, _("%s: Wrong type for input argument #%d: Argument %d in the list must be a matrix.\n"), "daskr", iPos + 1, iter + 1);
433                         DifferentialEquation::removeDifferentialEquationFunctions();
434                         FREE(pdYdotData);
435                         FREE(pdYData);
436                         FREE(YSize);
437                         return types::Function::Error;
438                     }
439
440                     vpDbl.push_back(pList->get(iter + 1)->getAs<types::Double>());
441                     YSize[sizeOfYSize + iter] = vpDbl[iter]->getSize();
442                     totalSize += YSize[sizeOfYSize + iter];
443                 }
444
445                 double* pdYDataTemp = pdYData;
446                 pdYData = (double*)MALLOC(totalSize * sizeof(double));
447                 C2F(dcopy)(&sizeOfpdYData, pdYDataTemp, &iOne, pdYData, &iOne);
448
449                 int position = sizeOfpdYData;
450                 for (int iter = 0; iter < pList->getSize() - 1; iter++)
451                 {
452                     C2F(dcopy)(&YSize[sizeOfYSize + iter], vpDbl[iter]->get(), &iOne, &pdYData[position], &iOne);
453                     position += vpDbl[iter]->getSize();
454                 }
455                 vpDbl.clear();
456                 sizeOfpdYData = totalSize;
457                 sizeOfYSize += pList->getSize() - 1;
458                 FREE(pdYDataTemp);
459                 FREE(sizeTemp);
460             }
461             else if (pList->get(0)->isCallable())
462             {
463                 // list(macro,...  function case
464                 if (bFuncF == false)
465                 {
466                     bFuncF = true;
467                     deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
468                     for (int iter = 1; iter < pList->getSize(); iter++)
469                     {
470                         deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
471                     }
472                 }
473                 else if (bFuncJac == false && pDblNg == NULL)
474                 {
475                     bFuncJac = true;
476                     deFunctionsManager->setJacFunction(pList->get(0)->getAs<types::Callable>());
477                     for (int iter = 1; iter < pList->getSize(); iter++)
478                     {
479                         deFunctionsManager->setJacArgs(pList->get(iter)->getAs<types::InternalType>());
480                     }
481                 }
482                 else if (bFuncG == false && pDblNg)
483                 {
484                     bFuncG = true;
485                     deFunctionsManager->setGFunction(pList->get(0)->getAs<types::Callable>());
486                     for (int iter = 1; iter < pList->getSize(); iter++)
487                     {
488                         deFunctionsManager->setGArgs(pList->get(iter)->getAs<types::InternalType>());
489                     }
490                 }
491                 else if (bFuncG && bFuncPsol == false)
492                 {
493                     bFuncPsol = true;
494                     deFunctionsManager->setPsolFunction(pList->get(0)->getAs<types::Callable>());
495                     for (int iter = 1; iter < pList->getSize(); iter++)
496                     {
497                         deFunctionsManager->setPsolArgs(pList->get(iter)->getAs<types::InternalType>());
498                     }
499                 }
500                 else if (bFuncPsol && bFuncPjac == false)
501                 {
502                     bFuncPjac = true;
503                     deFunctionsManager->setPjacFunction(pList->get(0)->getAs<types::Callable>());
504                     for (int iter = 1; iter < pList->getSize(); iter++)
505                     {
506                         deFunctionsManager->setPjacArgs(pList->get(iter)->getAs<types::InternalType>());
507                     }
508                 }
509             }
510             else if (pList->get(0)->isDouble() && bFuncF == true)
511             {
512                 // list(double,... then this list is the info argument
513                 if (pList->getSize() != 14)
514                 {
515                     Scierror(267, _("%s: Wrong size for input argument #%d: A list of size %d expected.\n"), "daskr", iPos + 1, 7);
516                     DifferentialEquation::removeDifferentialEquationFunctions();
517                     FREE(pdYdotData);
518                     FREE(pdYData);
519                     FREE(YSize);
520                     return types::Function::Error;
521                 }
522
523                 for (int i = 0; i < 14; i++) // info = list([],0,[],[],[],0,[],0,[],0,[],[],[],1)
524                 {
525                     if (pList->get(i)->isDouble() == false || (pList->get(i)->getAs<types::Double>()->isScalar() == false && (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)))
526                     {
527                         if (i == 1 || i == 5 || i == 7 || i == 9 || i == 13)
528                         {
529                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a scalar.\n"), "daskr", iPos + 1, i);
530                         }
531                         else
532                         {
533                             Scierror(999, _("%s: Wrong type for input argument #%d: Element %d in the info list must be a matrix.\n"), "daskr", iPos + 1, i);
534                         }
535                         DifferentialEquation::removeDifferentialEquationFunctions();
536                         FREE(pdYdotData);
537                         FREE(pdYData);
538                         FREE(YSize);
539                         return types::Function::Error;
540                     }
541                 }
542
543                 //  --   subvariable tstop(info) --
544                 types::Double* pDblTemp = pList->get(0)->getAs<types::Double>();
545                 if (pDblTemp->getSize() != 0)
546                 {
547                     info[3] = 1;
548                     tstop = pDblTemp->get(0);
549                 }
550
551                 //  --   subvariable imode(info) --
552                 info[2] = (int)pList->get(1)->getAs<types::Double>()->get(0);
553
554                 //  --   subvariable band(info) --
555                 pDblTemp = pList->get(2)->getAs<types::Double>();
556                 if (pDblTemp->getSize() == 2)
557                 {
558                     info[5] = 1;
559                     ml = (int)pDblTemp->get(0);
560                     mu = (int)pDblTemp->get(1);
561                     deFunctionsManager->setMl(ml);
562                     deFunctionsManager->setMu(mu);
563                 }
564                 else if (pDblTemp->getSize() != 0)
565                 {
566                     Scierror(267, _("%s: Wrong size for input argument #%d: Argument %d in te list must be of size %d.\n"), "daskr", iPos + 1, 3, 2);
567                     DifferentialEquation::removeDifferentialEquationFunctions();
568                     FREE(pdYdotData);
569                     FREE(pdYData);
570                     FREE(YSize);
571                     return types::Function::Error;
572                 }
573
574                 //  --   subvariable maxstep(info) --
575                 pDblTemp = pList->get(3)->getAs<types::Double>();
576                 if (pDblTemp->getSize() != 0)
577                 {
578                     info[6] = 1;
579                     maxstep = pDblTemp->get(0);
580                 }
581
582                 //  --   subvariable stepin(info) --
583                 pDblTemp = pList->get(4)->getAs<types::Double>();
584                 if (pDblTemp->getSize() != 0)
585                 {
586                     info[7] = 1;
587                     stepin = pDblTemp->get(0);
588                 }
589
590                 //  --   subvariable nonneg(info) --
591                 info[9]  = (int)pList->get(5)->getAs<types::Double>()->get(0);
592
593                 //  --   subvariable consistent(info) --
594                 pDblTemp = pList->get(6)->getAs<types::Double>();
595                 if (pDblTemp->getSize() != 0)
596                 {
597                     info[10] = 1;
598                     if (info[9] == 0 || info[9] == 2)
599                     {
600                         LID = 40;
601                     }
602                     else
603                     {
604                         LID = 40 + pDblX0->getRows();
605                     }
606                 }
607                 pDblE7 = pDblTemp;
608
609                 //  --   subvariable iteration(info) --
610                 if (pList->get(7)->getAs<types::Double>()->get(0) == 1)
611                 {
612                     info[11] = 1;
613                 }
614
615                 //  --   subvariable defaultKrylov(info) --
616                 pDblTemp = pList->get(8)->getAs<types::Double>();
617                 if (pDblTemp->getSize() == 0)
618                 {
619                     // maxl and kmp need default values
620                     maxl = min(5, pDblX0->getSize());
621                     kmp = maxl;
622                 }
623                 else
624                 {
625                     // info then looks like list(..., [maxl kmp nrmax epli],...)
626                     info[12] = 1;
627                     maxl  = (int)pDblTemp->get(0);
628                     kmp   = (int)pDblTemp->get(1);
629                     nrmax = (int)pDblTemp->get(2);
630                     epli  = pDblTemp->get(3);
631                 }
632
633                 //  --   subvariable justConsistentComp(info) --
634                 dTemp = pList->get(9)->getAs<types::Double>()->get(0);
635                 if (dTemp)
636                 {
637                     // Check that info(11) = 1, meaning that the provided initial values
638                     // are not consistent
639                     if (info[10] == 1)
640                     {
641                         info[13] = 1;
642                     }
643                 }
644
645                 //  --   subvariable psolJac(info) --
646                 dTemp = pList->get(10)->getAs<types::Double>()->get(0);
647                 if (dTemp)
648                 {
649                     info[14] = 1;
650                 }
651
652                 //  --   subvariable excludeAlgebraic(info) --
653                 pDblTemp = pList->get(11)->getAs<types::Double>();
654                 if (pDblTemp->getSize() != 0)
655                 {
656                     info[15] = 1;
657                     if (info[9] == 0 || info[9] == 2)
658                     {
659                         LID = 40;
660                     }
661                     else
662                     {
663                         LID = 40 + pDblX0->getRows();
664                     }
665                 }
666                 pDblE12 = pDblTemp;
667
668                 //  --   subvariable defaultHeuristic(info) --
669                 pDblTemp = pList->get(12)->getAs<types::Double>();
670                 if (pDblTemp->getSize() != 0)
671                 {
672                     // info then looks like list(..., [mxnit mxnj lsoff stptol epinit],...)
673                     info[16] = 1;
674                     mxnit   = (int)pDblTemp->get(0);
675                     mxnj    = (int)pDblTemp->get(1);
676                     mxnh    = (int)pDblTemp->get(2);
677                     lsoff   = (int)pDblTemp->get(3);
678                     stptol  = pDblTemp->get(4);
679                     epinit  = pDblTemp->get(5);
680                 }
681
682                 //  --   subvariable verbosity(info) --
683                 dTemp = pList->get(13)->getAs<types::Double>()->get(0);
684                 if (dTemp == 1 || dTemp == 2)
685                 {
686                     info[17] = (int)dTemp;
687                 }
688
689                 bListInfo = true;
690             }
691             else
692             {
693                 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"), "daskr", iPos + 1);
694                 DifferentialEquation::removeDifferentialEquationFunctions();
695                 FREE(pdYdotData);
696                 FREE(pdYData);
697                 FREE(YSize);
698                 return types::Function::Error;
699             }
700         }
701         else
702         {
703             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "daskr", iPos + 1);
704             DifferentialEquation::removeDifferentialEquationFunctions();
705             FREE(pdYdotData);
706             FREE(pdYData);
707             FREE(YSize);
708             return types::Function::Error;
709         }
710     }
711
712     if (bFuncF == false)
713     {
714         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 3);
715         DifferentialEquation::removeDifferentialEquationFunctions();
716         FREE(pdYdotData);
717         FREE(pdYData);
718         FREE(YSize);
719         return types::Function::Error;
720     }
721
722     if (pDblNg == NULL)
723     {
724         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 2);
725         DifferentialEquation::removeDifferentialEquationFunctions();
726         FREE(pdYdotData);
727         FREE(pdYData);
728         FREE(YSize);
729         return types::Function::Error;
730     }
731
732     if (bFuncG == false)
733     {
734         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "daskr", (int)in.size() + 1);
735         DifferentialEquation::removeDifferentialEquationFunctions();
736         FREE(pdYdotData);
737         FREE(pdYData);
738         FREE(YSize);
739         return types::Function::Error;
740     }
741
742     if (bFuncJac == true)
743     {
744         info[4] = 1;
745     }
746
747     // *** Initialization. ***
748     double t0    = pDblT0->get(0);
749     int idid     = 0;
750     int maxord   = 5;
751
752     //compute itol and set the tolerances rtol and atol.
753     double* rtol    = NULL;
754     double* atol    = NULL;
755     double rpar[2]  = {0, 0};
756     int ipar[30];
757     memset(ipar, 0x00, 30 * sizeof(int));
758
759     if (pDblAtol)
760     {
761         if (pDblAtol->isScalar())
762         {
763             atol  = (double*)MALLOC(sizeof(double));
764             *atol = pDblAtol->get(0);
765         }
766         else
767         {
768             atol    = pDblAtol->get();
769             info[1] = 1;
770         }
771     }
772     else
773     {
774         atol  = (double*)MALLOC(sizeof(double));
775         *atol = 1.e-7;
776     }
777
778     if (pDblRtol)
779     {
780         if (pDblRtol->isScalar())
781         {
782             rtol  = (double*)MALLOC(sizeof(double));
783             *rtol = pDblRtol->get(0);
784         }
785         else
786         {
787             rtol = pDblRtol->get();
788         }
789     }
790     else // if rtol is not given atol will be used as a scalar.
791     {
792         if (pDblAtol && pDblAtol->isScalar() == false) // info[1] == 1
793         {
794             double dblSrc = 1.e-9;
795             int iSize = pDblAtol->getSize();
796             int iZero = 0;
797
798             rtol = (double*)MALLOC(iSize * sizeof(double));
799             C2F(dcopy)(&iSize, &dblSrc, &iZero, rtol, &iOne);
800         }
801         else
802         {
803             rtol    = (double*)MALLOC(sizeof(double));
804             *rtol   = 1.e-9;
805         }
806     }
807
808     // Compute rwork, iwork size.
809     // Create them.
810     int rworksize   = 0;
811     int* iwork      = NULL;
812     double* rwork   = NULL;
813     int* root       = NULL;
814
815     // compute size of rwork
816     rworksize = 60;
817     if (info[11] == 0)
818     {
819         rworksize += max(maxord + 4, 7) * pDblX0->getRows() + 3 * ng;
820         if (info[5] == 0)
821         {
822             // For the full (dense) JACOBIAN case
823             rworksize += pDblX0->getRows() + pDblX0->getRows();
824         }
825         else if (info[4] == 1)
826         {
827             // For the banded user-defined JACOBIAN case
828             rworksize += (2 * ml + mu + 1) * pDblX0->getRows();
829         }
830         else if (info[4] == 0)
831         {
832             // For the banded finite-difference-generated JACOBIAN case
833             rworksize += (2 * ml + mu + 1) * pDblX0->getRows() + 2 * (pDblX0->getRows() / (ml + mu + 1) + 1);
834         }
835     }
836     else // info[11] == 1
837     {
838         // LENWP is the length ot the rwork segment containing
839         // the matrix elements of the preconditioner P
840
841         LENWP = 7 * pDblX0->getRows();
842         rworksize += (maxord + 5) * pDblX0->getRows() + 3 * ng
843                      + (maxl + 3 + min(1, maxl - kmp)) * pDblX0->getRows()
844                      + (maxl + 3) * maxl + 1 + LENWP;
845     }
846
847     if (info[15] == 1)
848     {
849         rworksize += pDblX0->getRows();
850     }
851
852     // compute size of iwork
853     // liw = 40, then augment size according to the case
854     int iworksize = 40;
855     if (info[11] == 0)
856     {
857         iworksize += pDblX0->getRows();
858     }
859     else // info[11] == 1
860     {
861         // LENIWP is the length ot the iwork segment containing
862         // the matrix indexes of the preconditioner P
863         // (compressed sparse row format)
864         LENIWP = 25 * pDblX0->getRows() + 1;
865         iworksize += LENIWP;
866     }
867
868     if (info[9] == 1 || info[9] == 3)
869     {
870         iworksize += pDblX0->getRows();
871     }
872
873     if (info[10] == 1 || info[15] == 1)
874     {
875         iworksize += pDblX0->getRows();
876     }
877
878     iwork = (int*)MALLOC(iworksize * sizeof(int));
879     rwork = (double*)MALLOC(rworksize * sizeof(double));
880     root  = (int*)MALLOC(ng * sizeof(int));
881
882     if (pDblHd != NULL)
883     {
884         if (iworksize + rworksize != pDblHd->getSize())
885         {
886             Scierror(77, _("%s: Wrong size for input argument(s) %d: %d expected.\n"), "daskr", in.size(), iworksize + rworksize);
887             DifferentialEquation::removeDifferentialEquationFunctions();
888             FREE(pdYdotData);
889             FREE(pdYData);
890             FREE(YSize);
891             FREE(iwork);
892             FREE(rwork);
893             if (pDblAtol == NULL || pDblAtol->isScalar())
894             {
895                 FREE(atol);
896             }
897             if (pDblRtol == NULL || pDblRtol->isScalar())
898             {
899                 FREE(rtol);
900             }
901             return types::Function::Error;
902         }
903
904         C2F(dcopy)(&rworksize, pDblHd->get(), &iOne, rwork, &iOne);
905
906         for (int i = 0; i < iworksize; i++)
907         {
908             iwork[i] = (int)pDblHd->get(rworksize + i);
909         }
910
911         info[0] = 1;
912     }
913
914     if (info[3] == 1)
915     {
916         rwork[0] = tstop;
917     }
918
919     if (info[6] == 1)
920     {
921         rwork[1] = maxstep;
922     }
923
924     if (info[7] == 1)
925     {
926         rwork[2] = stepin;
927     }
928
929     if (info[5] == 1)
930     {
931         iwork[0] = ml;
932         iwork[1] = mu;
933     }
934
935     if (info[10] == 1)
936     {
937         for (int i = 0; i < pDblX0->getRows(); i++)
938         {
939             iwork[LID + i] = pDblE7->get(i);
940         }
941     }
942
943     if (info[11] == 1)
944     {
945         iwork[26] = LENWP;
946         iwork[27] = LENIWP;
947     }
948
949     if (info[12] == 1)
950     {
951         iwork[23] = maxl;
952         iwork[24] = kmp;
953         iwork[25] = nrmax;
954         rwork[10] = epli;
955     }
956
957     if (info[14] == 1)
958     {
959         // Set ipar and rper
960         ipar[0] = pDblX0->getRows();
961         ipar[1] = pDblX0->getRows();
962         ipar[2] = 2;
963         ipar[3] = 2;
964         ipar[4] = 1;
965         rpar[0] = 0.005;
966         rpar[1] = 0.005;
967     }
968
969     if (info[15] == 1)
970     {
971         for (int i = 0; i < pDblX0->getRows(); i++)
972         {
973             iwork[LID + i] = pDblE12->get(i);
974         }
975     }
976
977     if (info[16] == 1)
978     {
979         iwork[31] = mxnit;
980         iwork[32] = mxnj;
981         iwork[33] = mxnh;
982         iwork[34] = lsoff;
983         rwork[13] = stptol;
984         rwork[14] = epinit;
985     }
986
987     iwork[16] = rworksize;
988     iwork[17] = iworksize;
989
990     // *** Perform operation. ***
991     std::list<types::Double*> lpDblOut;
992     int size = pDblX0->getRows();
993     int rowsOut = 1 + pDblX0->getRows() * 2;
994     int iret = 0;
995     int ididtmp = 0;
996
997     for (int i = 0; i < pDblT->getSize(); i++)
998     {
999         types::Double* pDblOut = new types::Double(rowsOut, 1);
1000         lpDblOut.push_back(pDblOut);
1001
1002         double t = pDblT->get(i);
1003         int pos  = 0;
1004         pDblOut->set(pos, t);
1005
1006         if (t == t0)
1007         {
1008             pos++;
1009             C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1010             pos += pDblX0->getRows();
1011             C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1012
1013             continue;
1014         }
1015
1016         try
1017         {
1018             // SUBROUTINE DDASKR(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
1019             //                  IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
1020             //                  RT, NRT, JROOT)
1021             C2F(ddaskr)(dassl_f, YSize, &t0, pdYData, pdYdotData, &t,
1022                         info, rtol, atol, &idid, rwork, &rworksize,
1023                         iwork, &iworksize, rpar, ipar,
1024                         info[14] == 1 ? (void*)daskr_pjac : (void*)dassl_jac,
1025                         daskr_psol, dasrt_g, &ng, root);
1026
1027             // values of idid says the same things that in ddasrt function,
1028             // except these two values.
1029             ididtmp = idid;
1030             // The integration to TSTOP was successfully completed (T=TSTOP)
1031             if (idid == 4)
1032             {
1033                 idid = 2;
1034             }
1035
1036             // One or more root found (eq to idid = 4 for dassl)
1037             if (idid == 5)
1038             {
1039                 idid = 4;
1040             }
1041
1042             if (idid == 6)
1043             {
1044                 idid = 3;
1045             }
1046
1047             // use the same error function that dasrt
1048             iret = checkDasslError(idid);
1049             if (iret == 1) // error
1050             {
1051                 Scierror(999, _("%s: ddaskr return with state %d.\n"), "daskr", ididtmp);
1052             }
1053         }
1054         catch (ScilabError &e)
1055         {
1056             char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
1057             sciprint(_("%s: exception caught in '%s' subroutine.\n"), "daskr", "ddaskr");
1058             Scierror(999, pstrMsg);
1059             // set iret to 1 for free allocated memory
1060             iret = 1;
1061         }
1062
1063         // free memory when error occurred
1064         if (iret == 1)
1065         {
1066             lpDblOut.clear();
1067             DifferentialEquation::removeDifferentialEquationFunctions();
1068             FREE(pdYdotData);
1069             FREE(pdYData);
1070             FREE(YSize);
1071             FREE(iwork);
1072             FREE(rwork);
1073             FREE(root);
1074             if (pDblAtol == NULL || pDblAtol->isScalar())
1075             {
1076                 FREE(atol);
1077             }
1078             if (pDblRtol == NULL || pDblRtol->isScalar())
1079             {
1080                 FREE(rtol);
1081             }
1082             return types::Function::Error;
1083         }
1084
1085         pos++;
1086         C2F(dcopy)(&size, pdYData, &iOne, pDblOut->get() + pos, &iOne);
1087         pos += size;
1088         C2F(dcopy)(&size, pdYdotData, &iOne, pDblOut->get() + pos, &iOne);
1089
1090         if (iret == 2) // warning
1091         {
1092             if (ididtmp == -1 || ididtmp == 5 || ididtmp == 6)
1093             {
1094                 pDblOut->set(0, t0);
1095             }
1096             break;
1097         }
1098
1099         // iret == 0, check idid flag
1100         if (idid == 1)
1101         {
1102             pDblOut->set(0, t0);
1103             i--;
1104         }
1105         else if (idid == -2)
1106         {
1107             t0 = t;
1108             i--;
1109         }
1110         else
1111         {
1112             t0 = t;
1113         }
1114
1115         info[0] = 1;
1116     }
1117
1118     // *** Return result in Scilab. ***
1119     types::Double* pDblOut = new types::Double(rowsOut, (int)lpDblOut.size());
1120
1121     int sizeOfList = (int)lpDblOut.size();
1122     for (int i = 0; i < sizeOfList; i++)
1123     {
1124         int pos = i * rowsOut;
1125         C2F(dcopy)(&rowsOut, lpDblOut.front()->get(), &iOne, pDblOut->get() + pos, &iOne);
1126         lpDblOut.pop_front();
1127     }
1128     out.push_back(pDblOut);
1129
1130     int sizeOfRoot = 1;
1131     for (int i = 0; i < ng; i++)
1132     {
1133         if (root[i])
1134         {
1135             sizeOfRoot++;
1136         }
1137     }
1138     types::Double* pDblRoot = new types::Double(1, sizeOfRoot);
1139     pDblRoot->set(0, t0);
1140     int j = 0;
1141     for (int i = 0; i < ng; i++)
1142     {
1143         if (root[i])
1144         {
1145             j++;
1146             pDblRoot->set(j, i + 1);
1147         }
1148     }
1149     out.push_back(pDblRoot);
1150
1151     if (_iRetCount == 3)
1152     {
1153         types::Double* pDblHdOut = new types::Double(rworksize + iworksize, 1);
1154         C2F(dcopy)(&rworksize, rwork, &iOne, pDblHdOut->get(), &iOne);
1155
1156         for (int i = 0; i < iworksize; i++)
1157         {
1158             pDblHdOut->set(rworksize + i, (double)iwork[i]);
1159         }
1160
1161         out.push_back(pDblHdOut);
1162     }
1163
1164     // *** FREE memory. ***
1165     if (pDblAtol == NULL || pDblAtol->isScalar())
1166     {
1167         FREE(atol);
1168     }
1169
1170     if (pDblRtol == NULL || pDblRtol->isScalar())
1171     {
1172         FREE(rtol);
1173     }
1174
1175     FREE(pdYData);
1176     FREE(pdYdotData);
1177     FREE(YSize);
1178     FREE(rwork);
1179     FREE(iwork);
1180     FREE(root);
1181
1182     DifferentialEquation::removeDifferentialEquationFunctions();
1183
1184     return types::Function::OK;
1185 }
1186