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