Merge remote-tracking branch 'origin/master' into windows
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_bvode.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
25 extern "C"
26 {
27 #include "sci_malloc.h"
28 #include "localization.h"
29 #include "Scierror.h"
30 #include "scifunctions.h"
31 #include "sciprint.h"
32 }
33
34 /*--------------------------------------------------------------------------*/
35 types::Function::ReturnValue sci_bvode(types::typed_list &in, int _iRetCount, types::typed_list &out)
36 {
37     int iPos = 0;
38
39     // input
40     types::Double* pDblXpts     = NULL;
41     types::Double* pDblM        = NULL;
42     types::Double* pDblZeta     = NULL;
43     types::Double* pDblIpar     = NULL;
44     types::Double* pDblLtol     = NULL;
45     types::Double* pDblTol      = NULL;
46     types::Double* pDblFixpnt   = NULL;
47
48     double aleft    = 0;
49     double aright   = 0;
50     int ncomp       = 0;
51     int sumM        = 0;
52     int maxM        = 0;
53     int iflag       = 0;
54     int ipar[11];
55
56     // error message catched
57     std::ostringstream os;
58     bool bCatch = false;
59
60     // *** check the minimal number of input args. ***
61     if (in.size() != 15)
62     {
63         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "bvode", 15);
64         return types::Function::Error;
65     }
66
67     // *** check number of output args ***
68     if (_iRetCount > 1)
69     {
70         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "bvode", 1);
71         return types::Function::Error;
72     }
73
74     // *** check type of input args and get it. ***
75
76     // xpoints
77     if (in[iPos]->isDouble() == false)
78     {
79         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
80         return types::Function::Error;
81     }
82
83     pDblXpts = in[iPos]->getAs<types::Double>();
84     iPos++;
85
86     // ncomp
87     if (in[iPos]->isDouble() == false)
88     {
89         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
90         return types::Function::Error;
91     }
92
93     types::Double* pDblN = in[iPos]->getAs<types::Double>();
94
95     if (pDblN->isScalar() == false)
96     {
97         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
98         return types::Function::Error;
99     }
100
101     ncomp = (int)pDblN->get(0);
102
103     if (ncomp > 20)
104     {
105         Scierror(999, _("%s: Wrong value for input argument #%d: Value at most 20 expected.\n"), "bvode", iPos + 1);
106         return types::Function::Error;
107     }
108     iPos++;
109
110     // m
111     if (in[iPos]->isDouble() == false)
112     {
113         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
114         return types::Function::Error;
115     }
116
117     pDblM = in[iPos]->getAs<types::Double>();
118
119     if (pDblM->getSize() != ncomp)
120     {
121         Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size %d (N) expected.\n"), "bvode", iPos + 1, ncomp);
122         return types::Function::Error;
123     }
124
125     int* M = (int*)MALLOC(pDblM->getSize() * sizeof(int));
126     for (int i = 0; i < pDblM->getSize(); i++)
127     {
128         M[i] = (int)pDblM->get(i);
129         sumM += (int)pDblM->get(i);
130         maxM = std::max(maxM, (int)pDblM->get(i));
131     }
132
133     if (sumM > 40)
134     {
135         Scierror(999, _("%s: Wrong value for input argument #%d: Sum of m (%d) must be less than 40.\n"), "bvode", iPos + 1, sumM);
136         FREE(M);
137         return types::Function::Error;
138     }
139
140     iPos++;
141
142     // aleft
143     if (in[iPos]->isDouble() == false)
144     {
145         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
146         FREE(M);
147         return types::Function::Error;
148     }
149
150     types::Double* pDblXLow = in[iPos]->getAs<types::Double>();
151
152     if (pDblXLow->isScalar() == false)
153     {
154         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
155         FREE(M);
156         return types::Function::Error;
157     }
158
159     aleft = pDblXLow->get(0);
160     iPos++;
161
162     // aright
163     if (in[iPos]->isDouble() == false)
164     {
165         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
166         FREE(M);
167         return types::Function::Error;
168     }
169
170     types::Double* pDblXUp = in[iPos]->getAs<types::Double>();
171
172     if (pDblXUp->isScalar() == false)
173     {
174         Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "bvode", iPos + 1);
175         FREE(M);
176         return types::Function::Error;
177     }
178
179     aright = pDblXUp->get(0);
180     iPos++;
181
182     // zeta
183     if (in[iPos]->isDouble() == false)
184     {
185         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
186         FREE(M);
187         return types::Function::Error;
188     }
189
190     pDblZeta = in[iPos]->getAs<types::Double>();
191
192     for (int i = 0; i < pDblZeta->getSize() - 1; i++)
193     {
194         if (pDblZeta->get(i) > pDblZeta->get(i + 1))
195         {
196             Scierror(999, _("%s: Wrong value for input argument #%d: zeta(j) lower or equal to zeta(j+1) expected.\n"), "bvode", iPos + 1);
197             FREE(M);
198             return types::Function::Error;
199         }
200     }
201
202     iPos++;
203
204     // ipar
205     if (in[iPos]->isDouble() == false)
206     {
207         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
208         FREE(M);
209         return types::Function::Error;
210     }
211
212     pDblIpar = in[iPos]->getAs<types::Double>();
213
214     if (pDblIpar->getSize() != 11)
215     {
216         Scierror(999, _("%s: Wrong size for input argument #%d: A vector of size 11 expected.\n"), "bvode", iPos + 1);
217         FREE(M);
218         return types::Function::Error;
219     }
220
221     for (int i = 0; i < 11; i++)
222     {
223         ipar[i] = (int)pDblIpar->get(i);
224     }
225
226     // check ipar arguments
227
228     // ipar(1) : 0 => linear  1 => nonlinear
229     if (ipar[0] != 0 && ipar[0] != 1)
230     {
231         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d or %d expected.\n"), "bvode", iPos + 1, 1, 0, 1);
232         FREE(M);
233         return types::Function::Error;
234     }
235
236     // ipar(2)=0 then collpnt is set to max(max(m(j))+1, 5-max(m(j))).
237     if (ipar[1] == 0)
238     {
239         ipar[1] = std::max(maxM + 1, 5 - maxM);
240     }
241
242     if (ipar[1] < 0 || ipar[1] > 7)
243     {
244         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d to %d expected.\n"), "bvode", iPos + 1, 2, 0, 7);
245         FREE(M);
246         return types::Function::Error;
247     }
248
249     // ipar(3)
250     if (ipar[2] == 0)
251     {
252         ipar[2] = 5;
253     }
254
255     // 0 < ipar(4) <= M
256     if (ipar[3] < 0 || ipar[3] > sumM)
257     {
258         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %s: %d to sum of m (=%d) expected.\n"), "bvode", iPos + 1, 4, 0, sumM);
259         FREE(M);
260         return types::Function::Error;
261     }
262
263     // bvode: ipar(5)>=nmax*nsizef where
264     // nsizef=4 + 3*M + (5+collpnt*N)*(collpnt*N+M) + (2*M-nrec)*2*M
265     int nmax = ipar[2];
266     int kdm  = ipar[1] * ncomp;
267     int nrec = 0;
268     for (int i = 0; i < sumM; i++)
269     {
270         int ib = sumM + 1 - i;
271         if (pDblZeta->get(ib - 1) >= aright)
272         {
273             nrec = i;
274         }
275     }
276
277     int nsizef = 4 + 3 * sumM + (5 + kdm) * (kdm + sumM) + (2 * sumM - nrec) * 2 * sumM;
278     if (ipar[4] < nmax * nsizef)
279     {
280         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d.\n"), "bvode", iPos + 1, 5);
281         FREE(M);
282         return types::Function::Error;
283     }
284
285     // ipar(6)>=nmax*nsizei where nsizei= 3 + collpnt*N + M.
286     int nsizei = 3 + kdm + sumM;
287     if (ipar[5] < nmax * nsizei)
288     {
289         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d.\n"), "bvode", iPos + 1, 6);
290         FREE(M);
291         return types::Function::Error;
292     }
293
294     // ipar(7) output control
295     if (ipar[6] < -1 || ipar[6] > 1)
296     {
297         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 7, -1, 1);
298         FREE(M);
299         return types::Function::Error;
300     }
301
302     // ipar(8)
303     if (ipar[7] < 0 || ipar[7] > 2)
304     {
305         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 8, 0, 2);
306         FREE(M);
307         return types::Function::Error;
308     }
309
310     // ipar(9)
311     if (ipar[8] < 0 || ipar[8] > 4)
312     {
313         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 9, 0, 4);
314         FREE(M);
315         return types::Function::Error;
316     }
317
318     // ipar(10)
319     if (ipar[9] < 0 || ipar[9] > 2)
320     {
321         Scierror(999, _("%s: Wrong value for input argument #%d: Wrong value for element %d: %d to %d expected.\n"), "bvode", iPos + 1, 9, 0, 2);
322         FREE(M);
323         return types::Function::Error;
324     }
325
326
327     iPos++;
328
329     // ltol
330     if (in[iPos]->isDouble() == false)
331     {
332         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
333         FREE(M);
334         return types::Function::Error;
335     }
336
337     pDblLtol = in[iPos]->getAs<types::Double>();
338
339     if (pDblLtol->getSize() != ipar[3])
340     {
341         Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(4)) expected.\n"), "bvode", iPos + 1, ipar[3]);
342         FREE(M);
343         return types::Function::Error;
344     }
345
346     // verify following cases :
347     // 1 <= ltol(1) < ltol(2) < ... < ltol(NTOL) <= M where M=sum(m)
348     // M is mstar and NTOL is the size of ltol
349     //
350     // 1 <= ltol(1) <= M
351     double* pdblLtol = pDblLtol->get();
352     if (pdblLtol[0] < 1 || pdblLtol[0] > sumM)
353     {
354         Scierror(999, _("%s: Wrong value for input argument #%d: %d to sum of m (=%d) expected.\n"), "bvode", iPos + 1, 1, sumM);
355         FREE(M);
356         return types::Function::Error;
357     }
358
359     // ltol(1) < ltol(2) < ... < ltol(NTOL) <= M
360     for (int i = 1; i < pDblLtol->getSize(); i++)
361     {
362         if (pdblLtol[i - 1] >= pdblLtol[i] || pdblLtol[i - 1] > sumM)
363         {
364             Scierror(999, _("%s: Wrong value for input argument #%d: Bad value for ltol(%d): ltol(1) < ltol(2) < ... < ltol(NTOL) <= M (sum of m) expected.\n"), "bvode", iPos + 1, i);
365             FREE(M);
366             return types::Function::Error;
367         }
368     }
369
370     iPos++;
371
372     // tol
373     if (in[iPos]->isDouble() == false)
374     {
375         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
376         FREE(M);
377         return types::Function::Error;
378     }
379
380     pDblTol = in[iPos]->getAs<types::Double>();
381
382     if (pDblTol->getSize() != ipar[3])
383     {
384         Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(4)) expected.\n"), "bvode", iPos + 1, ipar[3]);
385         FREE(M);
386         return types::Function::Error;
387     }
388
389     iPos++;
390
391     // fixpnt
392     if (in[iPos]->isDouble() == false)
393     {
394         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "bvode", iPos + 1);
395         FREE(M);
396         return types::Function::Error;
397     }
398
399     pDblFixpnt = in[iPos]->getAs<types::Double>();
400
401     if (pDblFixpnt->getSize() != ipar[10] && ipar[10] != 0)
402     {
403         Scierror(999, _("%s: Wrong size for input argument #%d: An array of size %d (ipar(11)) expected.\n"), "bvode", iPos + 1, ipar[10]);
404         FREE(M);
405         return types::Function::Error;
406     }
407
408     iPos++;
409
410     // functions: fsub,dfsub,gsub,dgsub,guess
411     DifferentialEquationFunctions deFunctionsManager("bvode");
412     DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
413     deFunctionsManager.setBvodeM(sumM);
414     deFunctionsManager.setBvodeN(ncomp);
415
416     for (int i = iPos; i < 15; i++)
417     {
418         if (in[i]->isCallable())
419         {
420             types::Callable* pCall = in[i]->getAs<types::Callable>();
421             if (i == 10) // fsub
422             {
423                 deFunctionsManager.setFsubFunction(pCall);
424             }
425             else if (i == 11) // dfsub
426             {
427                 deFunctionsManager.setDfsubFunction(pCall);
428             }
429             else if (i == 12) // gsub
430             {
431                 deFunctionsManager.setGsubFunction(pCall);
432             }
433             else if (i == 13) // dgsub
434             {
435                 deFunctionsManager.setDgsubFunction(pCall);
436             }
437             else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
438             {
439                 deFunctionsManager.setGuessFunction(pCall);
440             }
441         }
442         else if (in[i]->isString())
443         {
444             bool bOK = false;
445             types::String* pStr = in[i]->getAs<types::String>();
446             if (i == 10) // fsub
447             {
448                 bOK = deFunctionsManager.setFsubFunction(pStr);
449             }
450             else if (i == 11) // dfsub
451             {
452                 bOK = deFunctionsManager.setDfsubFunction(pStr);
453             }
454             else if (i == 12) // gsub
455             {
456                 bOK = deFunctionsManager.setGsubFunction(pStr);
457             }
458             else if (i == 13) // dgsub
459             {
460                 bOK = deFunctionsManager.setDgsubFunction(pStr);
461             }
462             else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
463             {
464                 bOK = deFunctionsManager.setGuessFunction(pStr);
465             }
466
467             if (bOK == false)
468             {
469                 const char* pst = pStr->get(0);
470                 Scierror(50, _("%s: Subroutine not found: %s\n"), "bvode", pst);
471                 DifferentialEquation::removeDifferentialEquationFunctions();
472                 FREE(M);
473                 return types::Function::Error;
474             }
475         }
476         else if (in[i]->isList())
477         {
478             types::List* pList = in[i]->getAs<types::List>();
479
480             if (pList->getSize() == 0)
481             {
482                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "bvode", i + 1, "(string empty)");
483                 DifferentialEquation::removeDifferentialEquationFunctions();
484                 FREE(M);
485                 return types::Function::Error;
486             }
487
488             if (pList->get(0)->isCallable())
489             {
490                 if (i == 10) // fsub
491                 {
492                     deFunctionsManager.setFsubFunction(pList->get(0)->getAs<types::Callable>());
493                     for (int iter = 1; iter < pList->getSize(); iter++)
494                     {
495                         deFunctionsManager.setFsubArgs(pList->get(iter)->getAs<types::InternalType>());
496                     }
497                 }
498                 else if (i == 11) // dfsub
499                 {
500                     deFunctionsManager.setDfsubFunction(pList->get(0)->getAs<types::Callable>());
501                     for (int iter = 1; iter < pList->getSize(); iter++)
502                     {
503                         deFunctionsManager.setDfsubArgs(pList->get(iter)->getAs<types::InternalType>());
504                     }
505                 }
506                 else if (i == 12) // gsub
507                 {
508                     deFunctionsManager.setGsubFunction(pList->get(0)->getAs<types::Callable>());
509                     for (int iter = 1; iter < pList->getSize(); iter++)
510                     {
511                         deFunctionsManager.setGsubArgs(pList->get(iter)->getAs<types::InternalType>());
512                     }
513                 }
514                 else if (i == 13) // dgsub
515                 {
516                     deFunctionsManager.setDgsubFunction(pList->get(0)->getAs<types::Callable>());
517                     for (int iter = 1; iter < pList->getSize(); iter++)
518                     {
519                         deFunctionsManager.setDgsubArgs(pList->get(iter)->getAs<types::InternalType>());
520                     }
521                 }
522                 else if (i == 14 && ipar[8] == 1) // guess is needed only if ipar(9) == 1
523                 {
524                     deFunctionsManager.setGuessFunction(pList->get(0)->getAs<types::Callable>());
525                     for (int iter = 1; iter < pList->getSize(); iter++)
526                     {
527                         deFunctionsManager.setGuessArgs(pList->get(iter)->getAs<types::InternalType>());
528                     }
529                 }
530             }
531             else
532             {
533                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "bvode", i + 1);
534                 DifferentialEquation::removeDifferentialEquationFunctions();
535                 FREE(M);
536                 return types::Function::Error;
537             }
538         }
539         else
540         {
541             Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "bvode", i + 1);
542             DifferentialEquation::removeDifferentialEquationFunctions();
543             FREE(M);
544             return types::Function::Error;
545         }
546     }
547
548     // *** Perform operation. ***
549     double* rwork   = (double*)MALLOC(ipar[4] * sizeof(double));
550     int* iwork      = (int*)MALLOC(ipar[5] * sizeof(int));
551     int* ltol       = (int*)MALLOC(pDblLtol->getSize() * sizeof(int));
552
553     for (int i = 0; i < pDblLtol->getSize(); i++)
554     {
555         ltol[i] = (int)pDblLtol->get(i);
556     }
557
558     try
559     {
560         C2F(colnew)(&ncomp, M, &aleft, &aright, pDblZeta->get(), ipar, ltol, pDblTol->get(), pDblFixpnt->get(), iwork, rwork, &iflag, bvode_fsub, bvode_dfsub, bvode_gsub, bvode_dgsub, bvode_guess);
561     }
562     catch (ast::InternalError &ie)
563     {
564         os << ie.GetErrorMessage();
565         bCatch = true;
566     }
567
568     if (bCatch)
569     {
570         FREE(iwork);
571         FREE(rwork);
572         FREE(M);
573         FREE(ltol);
574         DifferentialEquation::removeDifferentialEquationFunctions();
575
576         char szError[bsiz];
577         os_sprintf(szError, bsiz, _("%s: An error occured in '%s' subroutine.\n"), "bvode", "bvode");
578         os << szError;
579         throw ast::InternalError(os.str());
580     }
581
582     if (iflag != 1)
583     {
584         if (iflag == 0)
585         {
586             Scierror(999, _("%s: The collocation matrix is singular.\n"), "bvode");
587         }
588         else if (iflag == -1)
589         {
590             Scierror(999, _("%s: The expected no. of subintervals exceeds storage specifications.\n"), "bvode");
591         }
592         else if (iflag == -2)
593         {
594             Scierror(999, _("%s: The nonlinear iteration has not converged.\n"), "bvode");
595         }
596         else if (iflag == -3)
597         {
598             Scierror(999, _("%s: There is an input data error.\n"), "bvode");
599         }
600
601         FREE(iwork);
602         FREE(rwork);
603         FREE(M);
604         FREE(ltol);
605         DifferentialEquation::removeDifferentialEquationFunctions();
606
607         return types::Function::Error;
608     }
609
610     types::Double* pDblRes = new types::Double(sumM, pDblXpts->getSize());
611     double* res = (double*)MALLOC(pDblXpts->getSize() * sumM * sizeof(double));
612     for (int i = 0; i < pDblXpts->getSize(); i++)
613     {
614         double val = pDblXpts->get(i);
615         C2F(appsln)(&val, &res[i * sumM], rwork, iwork);
616     }
617
618     pDblRes->set(res);
619
620     FREE(iwork);
621     FREE(rwork);
622     FREE(M);
623     FREE(ltol);
624     FREE(res);
625     DifferentialEquation::removeDifferentialEquationFunctions();
626
627     // *** Return result in Scilab. ***
628     out.push_back(pDblRes);
629
630     return types::Function::OK;
631 }
632 /*--------------------------------------------------------------------------*/
633