bvode corrected.
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_intg2.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 "overload.hxx"
18 #include "execvisitor.hxx"
19 #include "callable.hxx"
20
21 extern "C"
22 {
23 #include "localization.h"
24 #include "Scierror.h"
25 #include "scifunctions.h"
26 #include "sci_warning.h"
27 #include "sciprint.h"
28 }
29
30 /*--------------------------------------------------------------------------*/
31 types::Function::ReturnValue sci_intg(types::typed_list &in, int _iRetCount, types::typed_list &out)
32 {
33     double pdA    = 0;
34     double pdB    = 0;
35     double pdEpsR = 1.0e-8;
36     double pdEpsA = 1.0e-14;
37
38     double result = 0;
39     double abserr = 0;
40
41     DifferentialEquationFunctions* deFunctionsManager = new DifferentialEquationFunctions(L"intg");
42     ConfigVariable::addDifferentialEquationFunctions(deFunctionsManager);
43
44     // *** check the minimal number of input args. ***
45     if (in.size() < 3 || in.size() > 5)
46     {
47         ScierrorW(77, _W("%ls: Wrong number of input argument(s): %d expected.\n"), L"intg", 3);
48         return types::Function::Error;
49     }
50
51     // *** check number of output args ***
52     if (_iRetCount > 2)
53     {
54         ScierrorW(78, _W("%ls: Wrong number of output argument(s): %d expected.\n"), L"intg", 2);
55         return types::Function::Error;
56     }
57
58     // *** check type of input args and get it. ***
59     // A
60     if (in[0]->isDouble() == false)
61     {
62         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"intg", 1);
63         return types::Function::Error;
64     }
65
66     types::Double* pDblA = in[0]->getAs<types::Double>();
67
68     if (pDblA->isScalar() == false)
69     {
70         ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A scalar expected.\n"), L"intg", 1);
71         return types::Function::Error;
72     }
73
74     pdA = pDblA->get(0);
75
76     // B
77     if (in[1]->isDouble() == false)
78     {
79         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"intg", 2);
80         return types::Function::Error;
81     }
82
83     types::Double* pDblB = in[1]->getAs<types::Double>();
84
85     if (pDblB->isScalar() == false)
86     {
87         ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A scalar expected.\n"), L"intg", 2);
88         return types::Function::Error;
89     }
90
91     pdB = pDblB->get(0);
92
93     // function
94     if (in[2]->isCallable())
95     {
96         types::Callable* pCall = in[2]->getAs<types::Callable>();
97         deFunctionsManager->setFFunction(pCall);
98
99     }
100     else if (in[2]->isString())
101     {
102         bool bOK = false;
103         types::String* pStr = in[2]->getAs<types::String>();
104         bOK = deFunctionsManager->setFFunction(pStr);
105
106         if (bOK == false)
107         {
108             ScierrorW(50, _W("%ls: Subroutine not found: %ls\n"), L"intg", pStr->get(0));
109             return types::Function::Error;
110         }
111     }
112     /*    else if(in[2]->isList())
113         {
114             types::List* pList = in[2]->getAs<types::List>();
115
116             if(pList->getSize() == 0)
117             {
118                 ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"intg", 3, L"(string empty)");
119                 return types::Function::Error;
120             }
121
122             if(pList->get(0)->isString())
123             {
124                 bool bOK = false;
125                 types::String* pStr = pList->get(0)->getAs<types::String>();
126                 bOK = deFunctionsManager->setOdeFFunction(pStr);
127
128                 if(bOK == false)
129                 {
130                     ScierrorW(50,_W("%ls: Argument #%d : Subroutine not found in list: %ls\n"), L"intg", 3, pStr->get(0));
131                     return types::Function::Error;
132                 }
133
134                 // realloc YSize to store the arguments size
135                 int totalSize = *YSize;
136                 free(YSize);
137                 YSize = (int*)malloc((pList->getSize() - 1) * sizeof(int));
138                 YSize[0] = totalSize;
139
140                 // get arguments and store their size in YSize
141                 std::vector<types::Double*> vpDbl;
142                 for(int iter = 0; iter < pList->getSize() - 1; iter++)
143                 {
144                     if(pList->get(iter + 1)->isDouble() == false)
145                     {
146                         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : Argument %d in the list must be a matrix.\n"), L"intg", 3, iter+1);
147                         return types::Function::Error;
148                     }
149
150                     vpDbl.push_back(pList->get(iter+1)->getAs<types::Double>());
151                     YSize[iter + 1] = vpDbl[iter]->getSize();
152                     totalSize += YSize[iter + 1];
153                 }
154
155                 // realloc Y to store arguments data in pdYData
156                 double* pdYDataTemp = pdYData;
157                 pdYData = (double*)malloc(totalSize * sizeof(double));
158                 C2F(dcopy)(&YSize[0], pdYDataTemp, &one, pdYData, &one);
159                 free(pdYDataTemp);
160
161                 int position = YSize[0];
162                 for(int iter = 0; iter < pList->getSize()-1; iter++)
163                 {
164                     C2F(dcopy)(&YSize[iter + 1], vpDbl[iter]->get(), &one, &pdYData[position], &one);
165                     position += vpDbl[iter]->getSize();
166                 }
167                 vpDbl.clear();
168             }
169             else if(pList->get(0)->isCallable())
170             {
171                 deFunctionsManager->setOdeFFunction(pList->get(0)->getAs<types::Callable>());
172                 for(int iter = 1; iter < pList->getSize(); iter++)
173                 {
174                     deFunctionsManager->setOdeFArgs(pList->get(iter)->getAs<types::InternalType>());
175                 }
176             }
177             else
178             {
179                 ScierrorW(999, _W("%ls: Wrong type for input argument #%d : The first argument in the list must be a string or a function.\n"), L"intg", 3);
180                 return types::Function::Error;
181             }
182         }
183     */
184     else
185     {
186         ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A function expected.\n"), L"intg", 3);
187         return types::Function::Error;
188     }
189
190     if (in.size() > 3)
191     {
192         if (in[3]->isDouble() == false)
193         {
194             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"intg", 4);
195             return types::Function::Error;
196         }
197
198         types::Double* pDblEpsR = in[3]->getAs<types::Double>();
199
200         if (pDblEpsR->isScalar() == false)
201         {
202             ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A scalar expected.\n"), L"intg", 4);
203             return types::Function::Error;
204         }
205
206         pdEpsR = pDblEpsR->get(0);
207     }
208
209     if (in.size() == 5)
210     {
211         if (in[4]->isDouble() == false)
212         {
213             ScierrorW(999, _W("%ls: Wrong type for input argument #%d : A matrix expected.\n"), L"intg", 5);
214             return types::Function::Error;
215         }
216
217         types::Double* pDblEpsA = in[4]->getAs<types::Double>();
218
219         if (pDblEpsA->isScalar() == false)
220         {
221             ScierrorW(999, _W("%ls: Wrong size for input argument #%d : A scalar expected.\n"), L"intg", 5);
222             return types::Function::Error;
223         }
224         pdEpsA = pDblEpsA->get(0);
225     }
226
227     // *** Create working table. ***
228     int limit = 3000 / 4;
229     int workSize = 3000;
230     int iworkSize = (3000 / 8 + 2) * 2 + 1;
231
232     double* work = (double*)malloc(workSize * sizeof(double));
233     int* iwork = (int*)malloc(iworkSize * sizeof(int));
234     double epsabs = fabs(pdEpsA);
235     double epsrel = fabs(pdEpsR);
236
237     // *** Perform operation. ***
238     //C2F(dqag0)(intg_f, &pdA, &pdB, &pdEpsa, &pdEpsr, val, abserr, stk(lpal), lw, stk(lpali), liw, ifail);
239
240     double* alist = work;
241     double* blist = work  + limit;
242     double* elist = blist + limit;
243     double* rlist = elist + limit;
244
245     int ier = 0;
246
247     C2F(dqags)(intg_f, &pdA, &pdB, &epsabs, &epsrel, alist, blist, elist, rlist, &limit, iwork, &iworkSize, &result, &abserr, &ier);
248
249     if (ier)
250     {
251         switch (ier)
252         {
253             case 1 :
254             {
255                 ScierrorW(999, _W("%ls: Maximum number of subdivisions allowed has been achieved.\n"), L"intg");
256                 return types::Function::Error;
257             }
258             case 2 :
259             {
260                 if (getWarningMode())
261                 {
262                     sciprintW(_W("%ls: Warning : The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be under-estimated.\n"), L"intg");
263                 }
264                 break;
265             }
266             case 3 :
267             {
268                 ScierrorW(999, _W("%ls: Extremely bad integrand behaviour occurs at some points of the integration interval.\n"), L"intg");
269                 return types::Function::Error;
270             }
271             case 4 :
272             {
273                 if (getWarningMode())
274                 {
275                     sciprintW(_W("%ls: Warning : The algorithm does not converge. Roundoff error is detected in the extrapolation table. It is presumed that the requested tolerance cannot be achieved, and that the returned result is the best which can be obtained.\n"), L"intg");
276                 }
277                 break;
278             }
279             case 5 :
280             {
281                 ScierrorW(999, _W("%ls: The integral is probably divergent, or slowly convergent.\n"), L"intg");
282                 return types::Function::Error;
283             }
284         }
285     }
286
287     // *** Return result in Scilab. ***
288     types::Double* pDblOut = new types::Double(result);
289     out.push_back(pDblOut);
290
291     if (_iRetCount == 2)
292     {
293         types::Double* pDblErrOut = new types::Double(abserr);
294         out.push_back(pDblErrOut);
295     }
296
297     // *** free. ***
298     free(work);
299     free(iwork);
300
301     ConfigVariable::removeDifferentialEquationFunctions();
302
303     return types::Function::OK;
304 }
305 /*--------------------------------------------------------------------------*/
306