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