differential_equations: fix dasrt/dassl memleaks
[scilab.git] / scilab / modules / differential_equations / sci_gateway / cpp / sci_int3d.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 "elem_common.h"
29 #include "localization.h"
30 #include "Scierror.h"
31 #include "scifunctions.h"
32 #include "configvariable_interface.h"
33 #include "sciprint.h"
34 }
35
36 /*--------------------------------------------------------------------------*/
37 types::Function::ReturnValue sci_int3d(types::typed_list &in, int _iRetCount, types::typed_list &out)
38 {
39     //input
40     types::Double* pDblX = NULL;
41     types::Double* pDblY = NULL;
42     types::Double* pDblZ = NULL;
43
44     int nf          = 1;
45     int minpts      = 0;
46     int maxpts      = 1000;
47     double epsabs   = 0.0;
48     double epsrel   = 1.0e-5;
49     int irestar     = 0;
50
51     // output
52     int ifail = 0;
53     int nevals = 0;
54
55     // error message catched
56     std::wostringstream os;
57     bool bCatch = false;
58
59     // *** check the minimal number of input args. ***
60     if (in.size() < 4 || in.size() > 6)
61     {
62         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "int3d", 4, 6);
63         return types::Function::Error;
64     }
65
66     // *** check number of output args according the methode. ***
67     if (_iRetCount > 3)
68     {
69         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "int3d", 2);
70         return types::Function::Error;
71     }
72
73     // *** check type of input args and get it. ***
74     // X
75     if (in[0]->isDouble() == false)
76     {
77         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 1);
78         return types::Function::Error;
79     }
80     pDblX = in[0]->getAs<types::Double>();
81     if (pDblX->isComplex())
82     {
83         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 1);
84         return types::Function::Error;
85     }
86
87     if (pDblX->getRows() != 4)
88     {
89         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 1);
90         return types::Function::Error;
91     }
92
93     // Y
94     if (in[1]->isDouble() == false)
95     {
96         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 2);
97         return types::Function::Error;
98     }
99     pDblY = in[1]->getAs<types::Double>();
100     if (pDblY->isComplex())
101     {
102         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 2);
103         return types::Function::Error;
104     }
105
106     if (pDblY->getRows() != 4)
107     {
108         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 2);
109         return types::Function::Error;
110     }
111
112     if (pDblY->getCols() != pDblX->getCols())
113     {
114         Scierror(999, _("%s: Wrong size for input argument #%d: Same size of input argument %d expected.\n"), "int3d", 2, 1);
115         return types::Function::Error;
116     }
117
118     // Z
119     if (in[2]->isDouble() == false)
120     {
121         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 3);
122         return types::Function::Error;
123     }
124     pDblZ = in[2]->getAs<types::Double>();
125     if (pDblZ->isComplex())
126     {
127         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 3);
128         return types::Function::Error;
129     }
130
131     if (pDblZ->getRows() != 4)
132     {
133         Scierror(999, _("%s: Wrong size for input argument #%d: A 4 by N matrix expected.\n"), "int3d", 4);
134         return types::Function::Error;
135     }
136
137     if (pDblZ->getCols() != pDblX->getCols())
138     {
139         Scierror(999, _("%s: Wrong size for input argument #%d: Same size of input argument %d expected.\n"), "int3d", 3, 1);
140         return types::Function::Error;
141     }
142
143     // function
144     DifferentialEquationFunctions deFunctionsManager(L"int3d");
145     DifferentialEquation::addDifferentialEquationFunctions(&deFunctionsManager);
146
147     if (in[3]->isCallable())
148     {
149         types::Callable* pCall = in[3]->getAs<types::Callable>();
150         deFunctionsManager.setFFunction(pCall);
151     }
152     else if (in[3]->isString())
153     {
154         bool bOK = false;
155         types::String* pStr = in[3]->getAs<types::String>();
156         bOK = deFunctionsManager.setFFunction(pStr);
157
158         if (bOK == false)
159         {
160             char* pst = wide_string_to_UTF8(pStr->get(0));
161             Scierror(50, _("%s: Subroutine not found: %s\n"), "int3d", pst);
162             FREE(pst);
163             DifferentialEquation::removeDifferentialEquationFunctions();
164             return types::Function::Error;
165         }
166     }
167     else if (in[3]->isList())
168     {
169         types::List* pList = in[3]->getAs<types::List>();
170
171         if (pList->getSize() == 0)
172         {
173             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "int3d", 4, "(string empty)");
174             DifferentialEquation::removeDifferentialEquationFunctions();
175             return types::Function::Error;
176         }
177
178         if (pList->get(0)->isCallable())
179         {
180             deFunctionsManager.setFFunction(pList->get(0)->getAs<types::Callable>());
181             for (int iter = 1; iter < pList->getSize(); iter++)
182             {
183                 deFunctionsManager.setFArgs(pList->get(iter)->getAs<types::InternalType>());
184             }
185         }
186         else
187         {
188             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a Scilab function.\n"), "int3d", 4);
189             DifferentialEquation::removeDifferentialEquationFunctions();
190             return types::Function::Error;
191         }
192     }
193     else
194     {
195         Scierror(999, _("%s: Wrong type for input argument #%d: A function expected.\n"), "int3d", 4);
196         DifferentialEquation::removeDifferentialEquationFunctions();
197         return types::Function::Error;
198     }
199
200     // nf (optional)
201     if (in.size() > 4)
202     {
203         if (in[4]->isDouble() == false)
204         {
205             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 5);
206             DifferentialEquation::removeDifferentialEquationFunctions();
207             return types::Function::Error;
208         }
209         types::Double* pDblNf = in[4]->getAs<types::Double>();
210         if (pDblNf->isComplex())
211         {
212             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 5);
213             DifferentialEquation::removeDifferentialEquationFunctions();
214             return types::Function::Error;
215         }
216
217         if (pDblNf->isScalar() == false)
218         {
219             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "int3d", 5);
220             DifferentialEquation::removeDifferentialEquationFunctions();
221             return types::Function::Error;
222         }
223         nf = (int)pDblNf->get(0);
224
225         if (nf < 1)
226         {
227             Scierror(999, _("%s: Wrong value for input argument #%d: A positive value expected.\n"), "int3d", 5);
228             DifferentialEquation::removeDifferentialEquationFunctions();
229             return types::Function::Error;
230         }
231     }
232
233     // params (optional)
234     if (in.size() == 6)
235     {
236         if (in[5]->isDouble() == false)
237         {
238             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 6);
239             DifferentialEquation::removeDifferentialEquationFunctions();
240             return types::Function::Error;
241         }
242
243         types::Double* pDblParams = in[5]->getAs<types::Double>();
244         if (pDblParams->isComplex())
245         {
246             Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "int3d", 6);
247             DifferentialEquation::removeDifferentialEquationFunctions();
248             return types::Function::Error;
249         }
250
251         if (pDblParams->getSize() != 4)
252         {
253             Scierror(999, _("%s: Wrong size for input argument #%d: %d expected.\n"), "int3d", 6, 4);
254             DifferentialEquation::removeDifferentialEquationFunctions();
255             return types::Function::Error;
256         }
257
258         if (getWarningMode())
259         {
260             for (int i = 0; i < 4; i++)
261             {
262                 if (pDblParams->get(i) < 0)
263                 {
264                     sciprint(_("%ls: Warning: Wrong value for the element %d of argument #%d: The default value will be used.\n"), L"int3d", i + 1, 6);
265                 }
266             }
267         }
268
269         minpts = pDblParams->get(0) < 0 ? minpts : (int)pDblParams->get(0);
270         maxpts = pDblParams->get(1) < 0 ? maxpts : (int)pDblParams->get(1);
271         epsabs = pDblParams->get(2) < 0.0 ? epsabs : pDblParams->get(2);
272         epsrel = pDblParams->get(3) < 0.0 ? epsrel : pDblParams->get(3);
273
274         if (pDblParams->get(2) == 0.0 && pDblParams->get(3) == 0.0)
275         {
276             if (getWarningMode())
277             {
278                 sciprint(_("%ls: Warning: Wrong value for the element %d and %d of argument #%d: The default value will be used.\n"), L"int3d", 3, 4, 6);
279             }
280             epsabs = 0.0;
281             epsrel = 1.0e-5;
282         }
283
284         if (minpts > maxpts)
285         {
286             Scierror(999, _("%s: Wrong value for input argument #%d: minpts smaller than maxpts expected.\n"), "int3d", 6);
287             DifferentialEquation::removeDifferentialEquationFunctions();
288             return types::Function::Error;
289         }
290     }
291
292     maxpts = std::max(maxpts, 43 * pDblX->getCols());
293
294     // *** Perform operation. ***
295     int one     = 1;
296     int three   = 3;
297     int maxsub  = 7 * (maxpts - 43 * pDblX->getCols()) / (8 * 43) + pDblX->getCols();
298     int size    = 4 * pDblX->getCols();
299
300     // input data: 3 rows, 4*maxsub cols
301     int dataSize = 3 * 4 * maxsub;
302     double* pdData = (double*)MALLOC(dataSize * sizeof(double));
303     C2F(dcopy)(&size, pDblX->get(), &one, pdData, &three);
304     C2F(dcopy)(&size, pDblY->get(), &one, pdData + 1, &three);
305     C2F(dcopy)(&size, pDblZ->get(), &one, pdData + 2, &three);
306
307     // output result
308     double* pdResult = (double*)MALLOC(nf * sizeof(double));
309
310     // output err
311     double* pdErr = (double*)MALLOC(nf * sizeof(double));
312
313     // workspace
314     int mdiv = 1; // default value, but can be changed on parallel computers. (see dcutet.f)
315     int dworkSize   = maxsub * (2 * nf + 1) + 7 * std::max(8 * mdiv, pDblX->getCols()) * nf + 1;
316     int iworkSize   = maxsub + mdiv;
317     double* dwork   = (double*)MALLOC(dworkSize * sizeof(double));
318     int* iwork      = (int*)MALLOC(iworkSize * sizeof(int));
319
320     int cols = pDblX->getCols();
321     try
322     {
323         C2F(dcutet)(int3d_f, &nf, pdData, &cols, &minpts, &maxpts, &epsabs, &epsrel, &maxsub, &dworkSize, &irestar, pdResult, pdErr, &nevals, &ifail, dwork, iwork);
324     }
325     catch (ast::InternalError &ie)
326     {
327         os << ie.GetErrorMessage();
328         bCatch = true;
329     }
330
331     FREE(pdData);
332     FREE(dwork);
333     FREE(iwork);
334     DifferentialEquation::removeDifferentialEquationFunctions();
335
336     if (bCatch)
337     {
338         wchar_t szError[bsiz];
339         os_swprintf(szError, bsiz, _W("%s: An error occurred in '%s' subroutine.\n").c_str(), "int3d", "dcutet");
340         os << szError;
341         throw ast::InternalError(os.str());
342     }
343
344     if (ifail)
345     {
346         if (ifail == 1)
347         {
348             if (getWarningMode())
349             {
350                 sciprint(_("%ls: Warning: maxpts was too small to obtain the required accuracy.\n"), L"int3d");
351             }
352         }
353         else if (ifail == 3)
354         {
355             Scierror(999, _("%s: The volume of one of the initially given tetrahedrons is zero.\n"), "int3d");
356             FREE(pdResult);
357             FREE(pdErr);
358             return types::Function::Error;
359         }
360         else // normaly nerver call.
361         {
362             Scierror(999, _("%s: dcutet return with error %d.\n"), "int3d", ifail);
363             FREE(pdResult);
364             FREE(pdErr);
365             return types::Function::Error;
366         }
367     }
368
369     // *** Return result in Scilab. ***
370     types::Double* pDblOut = new types::Double(nf, 1);
371     pDblOut->set(pdResult);
372     out.push_back(pDblOut);
373
374     if (_iRetCount > 1)
375     {
376         types::Double* pDblErrOut = new types::Double(nf, 1);
377         pDblErrOut->set(pdErr);
378         out.push_back(pDblErrOut);
379     }
380
381     if (_iRetCount == 3)
382     {
383         types::Double* pDblNevalsOut = new types::Double((double)nevals);
384         out.push_back(pDblNevalsOut);
385     }
386
387     FREE(pdResult);
388     FREE(pdErr);
389
390     return types::Function::OK;
391 }
392 /*--------------------------------------------------------------------------*/
393