Management of errors in differential equations functions.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_intg.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - 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 #include <cmath>
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
23 extern "C"
24 {
25 #include "MALLOC.h"
26 #include "localization.h"
27 #include "Scierror.h"
28 #include "scifunctions.h"
29 #include "warningmode.h"
30 #include "sciprint.h"
31 #include "matrix_division.h"
32 #include "vfinite.h"
33 }
34
35 /* ISNAN overloading for Mac OS X */
36 #ifdef __APPLE__
37 #undef ISNAN
38 #define ISNAN std::isnan
39 #endif
40
41 /*--------------------------------------------------------------------------*/
42 types::Function::ReturnValue sci_intg(types::typed_list &in, int _iRetCount, types::typed_list &out)
43 {
44     double pdA    = 0;
45     double pdB    = 0;
46     double pdEpsR = 1.0e-8;
47     double pdEpsA = 1.0e-14;
48
49     double result = 0;
50     double abserr = 0;
51
52     int iOne = 1;
53
54     // *** check the minimal number of input args. ***
55     if (in.size() < 3 || in.size() > 5)
56     {
57         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "intg", 3);
58         return types::Function::Error;
59     }
60
61     // *** check number of output args ***
62     if (_iRetCount > 3)
63     {
64         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "intg", 3);
65         return types::Function::Error;
66     }
67
68     // *** check type of input args and get it. ***
69     // A
70     if (in[0]->isDouble() == false)
71     {
72         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 1);
73         return types::Function::Error;
74     }
75
76     types::Double* pDblA = in[0]->getAs<types::Double>();
77
78     if (pDblA->isScalar() == false)
79     {
80         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 1);
81         return types::Function::Error;
82     }
83
84     pdA = pDblA->get(0);
85
86     if (ISNAN(pdA) || C2F(vfinite)(&iOne , &pdA) == false)
87     {
88         Scierror(264, _("%s: Wrong type for input argument #%d: Must not contain NaN or Inf.\n"), "intg", 1);
89         return types::Function::Error;
90     }
91
92     // B
93     if (in[1]->isDouble() == false)
94     {
95         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 2);
96         return types::Function::Error;
97     }
98
99     types::Double* pDblB = in[1]->getAs<types::Double>();
100
101     if (pDblB->isScalar() == false)
102     {
103         Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 2);
104         return types::Function::Error;
105     }
106
107     pdB = pDblB->get(0);
108
109     if (ISNAN(pdB) || C2F(vfinite)(&iOne , &pdB) == false)
110     {
111         Scierror(264, _("%s: Wrong type for input argument #%d: Must not contain NaN or Inf.\n"), "intg", 1);
112         return types::Function::Error;
113     }
114
115     // function
116     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"intg");
117     DifferentialEquation::addDifferentialEquationFunctions(deFunctionsManager);
118
119     if (in[2]->isCallable())
120     {
121         types::Callable* pCall = in[2]->getAs<types::Callable>();
122         deFunctionsManager->setFFunction(pCall);
123
124         // check function
125         double t = 1;
126         double ret = intg_f(&t);
127         if (ret == 0)
128         {
129             Scierror(50, _("%s: Argument #%d: Variable returned by scilab argument function is incorrect.\n"), "intg", 3);
130             DifferentialEquation::removeDifferentialEquationFunctions();
131             return types::Function::Error;
132         }
133     }
134     else if (in[2]->isString())
135     {
136         bool bOK = false;
137         types::String* pStr = in[2]->getAs<types::String>();
138         bOK = deFunctionsManager->setFFunction(pStr);
139
140         if (bOK == false)
141         {
142             char* pst = wide_string_to_UTF8(pStr->get(0));
143             Scierror(50, _("%s: Subroutine not found: %s\n"), "intg", pst);
144             FREE(pst);
145             DifferentialEquation::removeDifferentialEquationFunctions();
146             return types::Function::Error;
147         }
148     }
149     else if (in[2]->isList())
150     {
151         types::List* pList = in[2]->getAs<types::List>();
152
153         if (pList->getSize() == 0)
154         {
155             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "intg", 3, "(string empty)");
156             DifferentialEquation::removeDifferentialEquationFunctions();
157             return types::Function::Error;
158         }
159
160         if (pList->get(0)->isCallable())
161         {
162             deFunctionsManager->setFFunction(pList->get(0)->getAs<types::Callable>());
163             for (int iter = 1; iter < pList->getSize(); iter++)
164             {
165                 deFunctionsManager->setFArgs(pList->get(iter)->getAs<types::InternalType>());
166             }
167         }
168         else
169         {
170             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "intg", 3);
171             DifferentialEquation::removeDifferentialEquationFunctions();
172             return types::Function::Error;
173         }
174     }
175     else
176     {
177         Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "intg", 3);
178         DifferentialEquation::removeDifferentialEquationFunctions();
179         return types::Function::Error;
180     }
181
182     if (in.size() > 3)
183     {
184         if (in[3]->isDouble() == false)
185         {
186             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 4);
187             DifferentialEquation::removeDifferentialEquationFunctions();
188             return types::Function::Error;
189         }
190
191         types::Double* pDblEpsR = in[3]->getAs<types::Double>();
192
193         if (pDblEpsR->isScalar() == false)
194         {
195             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 4);
196             DifferentialEquation::removeDifferentialEquationFunctions();
197             return types::Function::Error;
198         }
199
200         pdEpsR = pDblEpsR->get(0);
201     }
202
203     if (in.size() == 5)
204     {
205         if (in[4]->isDouble() == false)
206         {
207             Scierror(999, _("%s: Wrong type for input argument #%d: A matrix expected.\n"), "intg", 5);
208             DifferentialEquation::removeDifferentialEquationFunctions();
209             return types::Function::Error;
210         }
211
212         types::Double* pDblEpsA = in[4]->getAs<types::Double>();
213
214         if (pDblEpsA->isScalar() == false)
215         {
216             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "intg", 5);
217             DifferentialEquation::removeDifferentialEquationFunctions();
218             return types::Function::Error;
219         }
220         pdEpsA = pDblEpsA->get(0);
221     }
222
223     // *** Create working table. ***
224     int limit   = 750;
225     int neval   = 0;
226     int last    = 0;
227     int lenw    = 4 * limit;
228
229     double* dwork   = (double*)MALLOC(lenw * sizeof(double));
230     int* iwork      = (int*)MALLOC(limit * sizeof(int));
231
232     double epsabs   = fabs(pdEpsA);
233     double epsrel   = fabs(pdEpsR);
234
235     // *** Perform operation. ***
236     int ier = 0;
237     try
238     {
239         C2F(dqags)(intg_f, &pdA, &pdB, &epsabs, &epsrel,
240                    &result, &abserr, &neval, &ier,
241                    &limit, &lenw, &last, iwork, dwork);
242     }
243     catch (ScilabError &e)
244     {
245         char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
246         sciprint(_("%s: exception caught in '%s' subroutine.\n"), "intg", "dqags");
247         Scierror(999, pstrMsg);
248         FREE(dwork);
249         FREE(iwork);
250         DifferentialEquation::removeDifferentialEquationFunctions();
251         return types::Function::Error;
252     }
253
254     FREE(dwork);
255     FREE(iwork);
256     DifferentialEquation::removeDifferentialEquationFunctions();
257
258     if (ier)
259     {
260         char* msg = NULL;
261         switch (ier)
262         {
263             case 1 :
264             {
265                 msg = _("%s: Maximum number of subdivisions achieved. Splitting the interval might help.\n");
266                 break;
267             }
268             case 2 :
269             {
270                 msg = _("%s: Round-off error detected, the requested tolerance (or default) cannot be achieved. Try using bigger tolerances.\n");
271                 break;
272             }
273             case 3 :
274             {
275                 msg = _("%s: Bad integrand behavior occurs at some points of the integration interval.\n");
276                 break;
277             }
278             case 4 :
279             {
280                 msg = _("%s: Convergence problem, round-off error detected. Try using bigger tolerances.\n");
281                 break;
282             }
283             case 5 :
284             {
285                 msg = _("%s: The integral is probably divergent, or slowly convergent.\n");
286                 break;
287             }
288             case 6 :
289             {
290                 msg = _("%s: Invalid input, absolute tolerance <= 0 and relative tolerance < 2.e-14.\n");
291                 break;
292             }
293             default :
294                 msg = _("%s: Convergence problem...\n");
295         }
296
297         if (_iRetCount == 3)
298         {
299             if (getWarningMode())
300             {
301                 sciprint(msg, "intg: Warning");
302             }
303         }
304         else
305         {
306             Scierror(999, msg, "intg: Error");
307             return types::Function::Error;
308         }
309     }
310
311     // *** Return result in Scilab. ***
312     types::Double* pDblOut = new types::Double(result);
313     out.push_back(pDblOut);
314
315     if (_iRetCount > 1)
316     {
317         out.push_back(new types::Double(abserr));
318     }
319
320     if (_iRetCount == 3)
321     {
322         out.push_back(new types::Double((double)ier));
323     }
324
325     return types::Function::OK;
326 }
327 /*--------------------------------------------------------------------------*/
328