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