* Bug #14067 fixed - 3rd argument of fsolve() became mandatory
[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     }
220
221     if (in.size() == 4)
222     {
223         if (in[3]->isDouble())
224         {
225             pDblTol = in[3]->getAs<types::Double>();
226             if (pDblTol->isScalar() == false)
227             {
228                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "fsolve", 4);
229                 return types::Function::Error;
230             }
231
232             dTol = pDblTol->get(0);
233         }
234         else
235         {
236             Scierror(999, _("%s: Wrong type for input argument #%d: A real expected.\n"), "fsolve", 4);
237             return types::Function::Error;
238         }
239     }
240
241     /*** perform operations ***/
242     // alloc working table
243     if (bJac)
244     {
245         iWorkSize = (iSizeX * (iSizeX + 13)) / 2;
246     }
247     else
248     {
249         iWorkSize = (iSizeX * (3 * iSizeX + 13)) / 2;
250     }
251
252     pdblWork = new double[iWorkSize];
253
254     // alloc output data
255     pDblV = new types::Double(pDblX->getDims(), pDblX->getDimsArray());
256
257     char const * pstrFunc = "fct";
258     try
259     {
260         if (bJac)
261         {
262             pstrFunc = "jac";
263             pdblJac = new double[iSizeX * iSizeX];
264             C2F(hybrj1)(jac, &iSizeX, pDblX->get(), pDblV->get(), pdblJac, &iSizeX, &dTol, &iInfo, pdblWork, &iWorkSize);
265         }
266         else
267         {
268             C2F(hybrd1)(fct, &iSizeX, pDblX->get(), pDblV->get(), &dTol, &iInfo, pdblWork, &iWorkSize);
269         }
270     }
271     catch (const ast::InternalError &e)
272     {
273         char* pstrMsg = wide_string_to_UTF8(e.GetErrorMessage().c_str());
274         sciprint(_("%s: exception caught in '%s' subroutine.\n"), "fsolve", pstrFunc);
275         Scierror(999, pstrMsg);
276         FREE(pstrMsg);
277         delete[] pdblWork;
278         delete pDblX;
279         if (pdblJac)
280         {
281             delete[] pdblJac;
282         }
283
284         return types::Function::Error;
285     }
286
287     delete[] pdblWork;
288     if (pdblJac)
289     {
290         delete[] pdblJac;
291     }
292
293     /*** return output arguments ***/
294     out.push_back(pDblX);
295
296     if (_iRetCount >= 2)
297     {
298         out.push_back(pDblV);
299     }
300     else
301     {
302         delete pDblV;
303     }
304
305     // info = 0   improper input parameters.
306     // info = 1   algorithm estimates that the relative error
307     // between x and the solution is at most tol.
308     // info = 2   number of calls to fcn with iflag = 1 has
309     // reached 100*(n+1).
310     // info = 3   tol is too small. no further improvement in
311     // the approximate solution x is possible.
312     // info = 4   iteration is not making good progress.
313     if (_iRetCount == 3)
314     {
315         out.push_back(new types::Double((double)iInfo));
316     }
317
318     return types::Function::OK;
319 }
320