Coverity #1321214, #1321215, #1321217, #1321220, #1321223, #1321224, #1321221, #13212...
[scilab.git] / scilab / modules / optimization / sci_gateway / cpp / sci_fsolve.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2013 - Scilab Enterprises - 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 "optimization_gw.hxx"
18 #include "function.hxx"
19 #include "double.hxx"
20 #include "string.hxx"
21 #include "list.hxx"
22 #include "optimizationfunctions.hxx"
23
24 extern "C"
25 {
26 #include "localization.h"
27 #include "Scierror.h"
28 #include "sciprint.h"
29 #include "scioptimfunctions.h"
30 #include "sci_malloc.h"
31 }
32 /*--------------------------------------------------------------------------*/
33
34 types::Function::ReturnValue sci_fsolve(types::typed_list &in, int _iRetCount, types::typed_list &out)
35 {
36     OptimizationFunctions* opFunctionsManager = NULL;
37
38     types::Double* pDblX    = NULL;
39     types::Double* pDblV    = NULL;
40     types::Double* pDblTol  = NULL;
41
42     int iSizeX      = 0;
43     int iWorkSize   = 0;
44     int iInfo       = 0;
45
46     double* pdblWork = NULL;
47     double* pdblJac  = NULL;
48
49     double dTol = 1.0e-10;
50     bool bJac   = false;
51
52     if (in.size() < 2 || in.size() > 4)
53     {
54         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "fsolve", 2, 4);
55         return types::Function::Error;
56     }
57
58     if (_iRetCount > 3)
59     {
60         Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "fsolve", 1, 3);
61         return types::Function::Error;
62     }
63
64     /*** get inputs arguments ***/
65     // get X
66     if (in[0]->isDouble() == false)
67     {
68         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "fsolve", 1);
69         return types::Function::Error;
70     }
71
72     pDblX = in[0]->clone()->getAs<types::Double>();
73
74     if (pDblX->isComplex())
75     {
76         Scierror(999, _("%s: Wrong type for input argument #%d: A real matrix expected.\n"), "fsolve", 1);
77         return types::Function::Error;
78     }
79
80     iSizeX = pDblX->getSize();
81
82     // get function
83     opFunctionsManager = new OptimizationFunctions(L"fsolve");
84     Optimization::addOptimizationFunctions(opFunctionsManager);
85     opFunctionsManager->setXRows(pDblX->getRows());
86     opFunctionsManager->setXCols(pDblX->getCols());
87
88     if (in[1]->isCallable())
89     {
90         types::Callable* pCall = in[1]->getAs<types::Callable>();
91         opFunctionsManager->setFsolveFctFunction(pCall);
92     }
93     else if (in[1]->isString())
94     {
95         types::String* pStr = in[1]->getAs<types::String>();
96         char* pst = wide_string_to_UTF8(pStr->get(0));
97         bool bOK = opFunctionsManager->setFsolveFctFunction(pStr);
98
99         if (bOK == false)
100         {
101             Scierror(50, _("%s: Subroutine not found: %s\n"), "fsolve", pst);
102             FREE(pst);
103             return types::Function::Error;
104         }
105
106         FREE(pst);
107     }
108     else if (in[1]->isList())
109     {
110         types::List* pList = in[1]->getAs<types::List>();
111         if (pList->getSize() == 0)
112         {
113             Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "fsolve", 2, "(string empty)");
114             return types::Function::Error;
115         }
116
117         if (pList->get(0)->isString())
118         {
119             types::String* pStr = pList->get(0)->getAs<types::String>();
120             char* pst = wide_string_to_UTF8(pStr->get(0));
121             bool bOK = opFunctionsManager->setFsolveFctFunction(pStr);
122
123             if (bOK == false)
124             {
125                 Scierror(50, _("%s: Subroutine not found: %s\n"), "fsolve", pst);
126                 FREE(pst);
127                 return types::Function::Error;
128             }
129
130             FREE(pst);
131         }
132         else if (pList->get(0)->isCallable())
133         {
134             types::Callable* pCall = pList->get(0)->getAs<types::Callable>();
135             opFunctionsManager->setFsolveFctFunction(pCall);
136             for (int iter = 1; iter < pList->getSize(); iter++)
137             {
138                 opFunctionsManager->setFsolveFctArgs(pList->get(iter)->getAs<types::InternalType>());
139             }
140         }
141         else
142         {
143             Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "fsolve", 2);
144             return types::Function::Error;
145         }
146     }
147     else
148     {
149         Scierror(999, _("%s: Wrong type for input argument #%d: A matrix or a function expected.\n"), "fsolve", 2);
150         return types::Function::Error;
151     }
152
153     if (in.size() >= 3)
154     {
155         if (in[2]->isCallable())
156         {
157             types::Callable* pCall = in[2]->getAs<types::Callable>();
158             opFunctionsManager->setFsolveJacFunction(pCall);
159
160             bJac = true;
161         }
162         else if (in[2]->isString())
163         {
164             types::String* pStr = in[2]->getAs<types::String>();
165             char* pst = wide_string_to_UTF8(pStr->get(0));
166             bool bOK = opFunctionsManager->setFsolveJacFunction(pStr);
167
168             if (bOK == false)
169             {
170                 Scierror(50, _("%s: Subroutine not found: %s\n"), "fsolve", pst);
171                 FREE(pst);
172                 return types::Function::Error;
173             }
174
175             bJac = true;
176             FREE(pst);
177         }
178         else if (in[2]->isList())
179         {
180             types::List* pList = in[2]->getAs<types::List>();
181             if (pList->getSize() == 0)
182             {
183                 Scierror(50, _("%s: Argument #%d: Subroutine not found in list: %s\n"), "fsolve", 3, "(string empty)");
184                 return types::Function::Error;
185             }
186
187             if (pList->get(0)->isString())
188             {
189                 types::String* pStr = pList->get(0)->getAs<types::String>();
190                 char* pst = wide_string_to_UTF8(pStr->get(0));
191                 bool bOK = opFunctionsManager->setFsolveJacFunction(pStr);
192
193                 if (bOK == false)
194                 {
195                     Scierror(50, _("%s: Subroutine not found: %s\n"), "fsolve", pst);
196                     return types::Function::Error;
197                 }
198
199                 bJac = true;
200                 FREE(pst);
201             }
202             else if (pList->get(0)->isCallable())
203             {
204                 types::Callable* pCall = pList->get(0)->getAs<types::Callable>();
205                 opFunctionsManager->setFsolveJacFunction(pCall);
206                 for (int iter = 1; iter < pList->getSize(); iter++)
207                 {
208                     opFunctionsManager->setFsolveJacArgs(pList->get(iter)->getAs<types::InternalType>());
209                 }
210
211                 bJac = true;
212             }
213             else
214             {
215                 Scierror(999, _("%s: Wrong type for input argument #%d: The first argument in the list must be a string or a function.\n"), "fsolve", 3);
216                 return types::Function::Error;
217             }
218         }
219         else if (in[2]->isDouble() && in.size() == 3)
220         {
221             pDblTol = in[2]->getAs<types::Double>();
222             if (pDblTol->isScalar() == false)
223             {
224                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "fsolve", 3);
225                 return types::Function::Error;
226             }
227
228             dTol = pDblTol->get(0);
229         }
230         else
231         {
232             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar or a function expected.\n"), "fsolve", 3);
233             return types::Function::Error;
234         }
235     }
236
237     if (in.size() == 4)
238     {
239         if (in[3]->isDouble())
240         {
241             pDblTol = in[3]->getAs<types::Double>();
242             if (pDblTol->isScalar() == false)
243             {
244                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "fsolve", 4);
245                 return types::Function::Error;
246             }
247
248             dTol = pDblTol->get(0);
249         }
250         else
251         {
252             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar or a function expected.\n"), "fsolve", 4);
253             return types::Function::Error;
254         }
255     }
256
257     /*** perform operations ***/
258     // alloc working table
259     if (bJac)
260     {
261         iWorkSize = (iSizeX * (iSizeX + 13)) / 2;
262     }
263     else
264     {
265         iWorkSize = (iSizeX * (3 * iSizeX + 13)) / 2;
266     }
267
268     pdblWork = new double[iWorkSize];
269
270     // alloc output data
271     pDblV = new types::Double(pDblX->getDims(), pDblX->getDimsArray());
272
273     char const * pstrFunc = "fct";
274     try
275     {
276         if (bJac)
277         {
278             pstrFunc = "jac";
279             pdblJac = new double[iSizeX * iSizeX];
280             C2F(hybrj1)(jac, &iSizeX, pDblX->get(), pDblV->get(), pdblJac, &iSizeX, &dTol, &iInfo, pdblWork, &iWorkSize);
281         }
282         else
283         {
284             C2F(hybrd1)(fct, &iSizeX, pDblX->get(), pDblV->get(), &dTol, &iInfo, pdblWork, &iWorkSize);
285         }
286     }
287     catch (const ast::InternalError &e)
288     {
289         char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
290         sciprint(_("%s: exception caught in '%s' subroutine.\n"), "fsolve", pstrFunc);
291         Scierror(999, pstrMsg);
292         FREE(pstrMsg);
293         delete[] pdblWork;
294         delete pDblX;
295         if (pdblJac)
296         {
297             delete[] pdblJac;
298         }
299
300         return types::Function::Error;
301     }
302
303     delete[] pdblWork;
304     if (pdblJac)
305     {
306         delete[] pdblJac;
307     }
308
309     /*** return output arguments ***/
310     out.push_back(pDblX);
311
312     if (_iRetCount >= 2)
313     {
314         out.push_back(pDblV);
315     }
316     else
317     {
318         delete pDblV;
319     }
320
321     // info = 0   improper input parameters.
322     // info = 1   algorithm estimates that the relative error
323     // between x and the solution is at most tol.
324     // info = 2   number of calls to fcn with iflag = 1 has
325     // reached 100*(n+1).
326     // info = 3   tol is too small. no further improvement in
327     // the approximate solution x is possible.
328     // info = 4   iteration is not making good progress.
329     if (_iRetCount == 3)
330     {
331         out.push_back(new types::Double((double)iInfo));
332     }
333
334     return types::Function::OK;
335 }
336