Mercurial > hg > octave-nkf
annotate src/DLD-FUNCTIONS/fsolve.cc @ 11729:de826649dfa2 release-3-0-x
fix to enable compiling with Intel C++
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Thu, 03 Apr 2008 10:24:18 -0400 |
| parents | 087af2a4ca26 |
| children | 7ca2735d74c2 72830070a17b |
| rev | line source |
|---|---|
| 2928 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, |
| 4 2007 John W. Eaton | |
| 2928 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 2928 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 2928 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <string> | |
| 29 | |
| 3567 | 30 #include <iomanip> |
| 3523 | 31 #include <iostream> |
| 5856 | 32 #include <sstream> |
| 2928 | 33 |
| 34 #include "NLEqn.h" | |
| 35 | |
| 36 #include "defun-dld.h" | |
| 37 #include "error.h" | |
| 38 #include "gripes.h" | |
| 39 #include "oct-obj.h" | |
| 2968 | 40 #include "ov-fcn.h" |
| 5729 | 41 #include "ov-cell.h" |
| 2928 | 42 #include "pager.h" |
| 3243 | 43 #include "unwind-prot.h" |
| 2928 | 44 #include "utils.h" |
| 45 #include "variables.h" | |
| 46 | |
| 3998 | 47 #include "NLEqn-opts.cc" |
| 48 | |
| 2928 | 49 // Global pointer for user defined function required by hybrd1. |
| 2968 | 50 static octave_function *fsolve_fcn; |
| 2928 | 51 |
| 4140 | 52 // Global pointer for optional user defined jacobian function. |
| 53 static octave_function *fsolve_jac; | |
| 54 | |
| 4139 | 55 // Have we warned about imaginary values returned from user function? |
| 4140 | 56 static bool warned_fcn_imaginary = false; |
| 57 static bool warned_jac_imaginary = false; | |
| 4139 | 58 |
| 3243 | 59 // Is this a recursive call? |
| 60 static int call_depth = 0; | |
| 61 | |
| 5275 | 62 octave_idx_type |
| 63 hybrd_info_to_fsolve_info (octave_idx_type info) | |
| 2928 | 64 { |
| 65 switch (info) | |
| 66 { | |
| 67 case -1: | |
| 68 break; | |
| 69 | |
| 70 case 0: | |
| 11611 | 71 info = -2; |
| 2928 | 72 break; |
| 73 | |
| 74 case 1: | |
| 75 break; | |
| 76 | |
| 77 case 2: | |
| 78 info = 4; | |
| 79 break; | |
| 80 | |
| 81 case 3: | |
| 82 case 4: | |
| 83 case 5: | |
| 84 info = 3; | |
| 85 break; | |
| 86 | |
| 87 default: | |
| 5856 | 88 { |
| 89 std::ostringstream buf; | |
| 90 buf << "fsolve: unrecognized value of INFO from MINPACK (= " | |
| 91 << info << ")"; | |
| 92 std::string msg = buf.str (); | |
| 93 warning (msg.c_str ()); | |
| 94 } | |
| 2928 | 95 break; |
| 96 } | |
| 97 | |
| 98 return info; | |
| 99 } | |
| 100 | |
| 101 ColumnVector | |
| 102 fsolve_user_function (const ColumnVector& x) | |
| 103 { | |
| 104 ColumnVector retval; | |
| 105 | |
| 5275 | 106 octave_idx_type n = x.length (); |
| 2928 | 107 |
| 108 octave_value_list args; | |
| 109 args.resize (1); | |
| 110 | |
| 111 if (n > 1) | |
| 112 { | |
| 113 Matrix m (n, 1); | |
| 5275 | 114 for (octave_idx_type i = 0; i < n; i++) |
| 2928 | 115 m (i, 0) = x (i); |
| 116 octave_value vars (m); | |
| 117 args(0) = vars; | |
| 118 } | |
| 119 else | |
| 120 { | |
| 121 double d = x (0); | |
| 122 octave_value vars (d); | |
| 123 args(0) = vars; | |
| 124 } | |
| 125 | |
| 126 if (fsolve_fcn) | |
| 127 { | |
| 3544 | 128 octave_value_list tmp = fsolve_fcn->do_multi_index_op (1, args); |
| 4139 | 129 |
| 2928 | 130 if (tmp.length () > 0 && tmp(0).is_defined ()) |
| 131 { | |
| 4140 | 132 if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) |
| 4139 | 133 { |
| 134 warning ("fsolve: ignoring imaginary part returned from user-supplied function"); | |
| 4140 | 135 warned_fcn_imaginary = true; |
| 4139 | 136 } |
| 137 | |
| 3419 | 138 retval = ColumnVector (tmp(0).vector_value ()); |
| 2928 | 139 |
| 140 if (error_state || retval.length () <= 0) | |
| 141 gripe_user_supplied_eval ("fsolve"); | |
|
11654
087af2a4ca26
fsolve: detect nonsquare systems
John W. Eaton <jwe@octave.org>
parents:
11611
diff
changeset
|
142 else if (retval.length () != x.length ()) |
|
087af2a4ca26
fsolve: detect nonsquare systems
John W. Eaton <jwe@octave.org>
parents:
11611
diff
changeset
|
143 error ("fsolve: unable to solve non-square systems"); |
| 2928 | 144 } |
| 145 else | |
| 146 gripe_user_supplied_eval ("fsolve"); | |
| 147 } | |
| 148 | |
| 149 return retval; | |
| 150 } | |
| 151 | |
| 4140 | 152 Matrix |
| 153 fsolve_user_jacobian (const ColumnVector& x) | |
| 154 { | |
| 155 Matrix retval; | |
| 156 | |
| 5275 | 157 octave_idx_type n = x.length (); |
| 4140 | 158 |
| 159 octave_value_list args; | |
| 160 args.resize (1); | |
| 161 | |
| 162 if (n > 1) | |
| 163 { | |
| 164 Matrix m (n, 1); | |
| 5275 | 165 for (octave_idx_type i = 0; i < n; i++) |
| 4628 | 166 m(i,0) = x(i); |
| 4140 | 167 octave_value vars (m); |
| 168 args(0) = vars; | |
| 169 } | |
| 170 else | |
| 171 { | |
| 4628 | 172 double d = x(0); |
| 4140 | 173 octave_value vars (d); |
| 174 args(0) = vars; | |
| 175 } | |
| 176 | |
| 177 if (fsolve_fcn) | |
| 178 { | |
| 179 octave_value_list tmp = fsolve_jac->do_multi_index_op (1, args); | |
| 180 | |
| 181 if (tmp.length () > 0 && tmp(0).is_defined ()) | |
| 182 { | |
| 183 if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) | |
| 184 { | |
| 185 warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); | |
| 186 warned_fcn_imaginary = true; | |
| 187 } | |
| 188 | |
| 189 retval = tmp(0).matrix_value (); | |
| 190 | |
| 191 if (error_state || retval.length () <= 0) | |
| 192 gripe_user_supplied_eval ("fsolve"); | |
|
11654
087af2a4ca26
fsolve: detect nonsquare systems
John W. Eaton <jwe@octave.org>
parents:
11611
diff
changeset
|
193 else if (! (retval.rows () == x.length () |
|
087af2a4ca26
fsolve: detect nonsquare systems
John W. Eaton <jwe@octave.org>
parents:
11611
diff
changeset
|
194 && retval.columns () == x.length ())) |
|
087af2a4ca26
fsolve: detect nonsquare systems
John W. Eaton <jwe@octave.org>
parents:
11611
diff
changeset
|
195 error ("fsolve: invalid Jacobian matrix dimensions"); |
| 4140 | 196 } |
| 197 else | |
| 198 gripe_user_supplied_eval ("fsolve"); | |
| 199 } | |
| 200 | |
| 201 return retval; | |
| 202 } | |
| 203 | |
| 3323 | 204 #define FSOLVE_ABORT() \ |
| 205 do \ | |
| 206 { \ | |
| 207 unwind_protect::run_frame ("Ffsolve"); \ | |
| 208 return retval; \ | |
| 209 } \ | |
| 210 while (0) | |
| 211 | |
| 212 #define FSOLVE_ABORT1(msg) \ | |
| 213 do \ | |
| 214 { \ | |
| 3747 | 215 ::error ("fsolve: " msg); \ |
| 3323 | 216 FSOLVE_ABORT (); \ |
| 217 } \ | |
| 218 while (0) | |
| 219 | |
| 220 #define FSOLVE_ABORT2(fmt, arg) \ | |
| 221 do \ | |
| 222 { \ | |
| 3747 | 223 ::error ("fsolve: " fmt, arg); \ |
| 3323 | 224 FSOLVE_ABORT (); \ |
| 225 } \ | |
| 226 while (0) | |
| 227 | |
| 2928 | 228 DEFUN_DLD (fsolve, args, nargout, |
| 3368 | 229 "-*- texinfo -*-\n\ |
| 7279 | 230 @deftypefn {Loadable Function} {[@var{x}, @var{fval}, @var{info}] =} fsolve (@var{fcn}, @var{x0})\n\ |
| 3368 | 231 Given @var{fcn}, the name of a function of the form @code{f (@var{x})}\n\ |
| 232 and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\ | |
| 233 equations such that @code{f(@var{x}) == 0}.\n\ | |
| 3964 | 234 \n\ |
| 11611 | 235 On return, @var{fval} contains the value of the function @var{fcn}\n\ |
| 236 evaluated at @var{x}, and @var{info} may be one of the following values:\n\ | |
| 237 \n\ | |
| 238 @table @asis\n\ | |
| 239 \n\ | |
| 240 @item -2\n\ | |
| 241 Invalid input parameters.\n\ | |
| 242 @item -1\n\ | |
| 243 Error in user-supplied function.\n\ | |
| 244 @item 1\n\ | |
| 245 Relative error between two consecutive iterates is at most the\n\ | |
| 246 specified tolerance (see @code{fsolve_options}).\n\ | |
| 247 @item 3\n\ | |
| 248 Algorithm failed to converge.\n\ | |
| 249 @item 4\n\ | |
| 250 Limit on number of function calls reached.\n\ | |
| 251 @end table\n\ | |
| 252 \n\ | |
| 5729 | 253 If @var{fcn} is a two-element string array, or a two element cell array\n\ |
| 254 containing either the function name or inline or function handle. The\n\ | |
| 255 first element names the function @math{f} described above, and the second\n\ | |
| 256 element names a function of the form @code{j (@var{x})} to compute the\n\ | |
| 257 Jacobian matrix with elements\n\ | |
| 4144 | 258 @tex\n\ |
| 259 $$ J = {\\partial f_i \\over \\partial x_j} $$\n\ | |
| 260 @end tex\n\ | |
| 261 @ifinfo\n\ | |
| 262 \n\ | |
| 263 @example\n\ | |
| 4197 | 264 df_i\n\ |
| 265 jac(i,j) = ----\n\ | |
| 266 dx_j\n\ | |
| 4144 | 267 @end example\n\ |
| 268 @end ifinfo\n\ | |
| 269 \n\ | |
| 3964 | 270 You can use the function @code{fsolve_options} to set optional\n\ |
| 271 parameters for @code{fsolve}.\n\ | |
| 3368 | 272 @end deftypefn") |
| 2928 | 273 { |
| 274 octave_value_list retval; | |
| 275 | |
| 4140 | 276 warned_fcn_imaginary = false; |
| 277 warned_jac_imaginary = false; | |
| 4139 | 278 |
| 3243 | 279 unwind_protect::begin_frame ("Ffsolve"); |
| 2928 | 280 |
| 3243 | 281 unwind_protect_int (call_depth); |
| 282 call_depth++; | |
| 283 | |
| 284 if (call_depth > 1) | |
| 3323 | 285 FSOLVE_ABORT1 ("invalid recursive call"); |
| 2928 | 286 |
| 3243 | 287 int nargin = args.length (); |
| 2928 | 288 |
| 3243 | 289 if (nargin == 2 && nargout < 4) |
| 290 { | |
| 5729 | 291 std::string fcn_name, fname, jac_name, jname; |
| 4140 | 292 fsolve_fcn = 0; |
| 293 fsolve_jac = 0; | |
| 294 | |
| 295 octave_value f_arg = args(0); | |
| 296 | |
| 5729 | 297 if (f_arg.is_cell ()) |
| 298 { | |
| 299 Cell c = f_arg.cell_value (); | |
| 300 if (c.length() == 1) | |
| 301 f_arg = c(0); | |
| 302 else if (c.length() == 2) | |
| 303 { | |
| 304 if (c(0).is_function_handle () || c(0).is_inline_function ()) | |
| 305 fsolve_fcn = c(0).function_value (); | |
| 306 else | |
| 307 { | |
| 308 fcn_name = unique_symbol_name ("__fsolve_fcn__"); | |
| 309 fname = "function y = "; | |
| 310 fname.append (fcn_name); | |
| 311 fname.append (" (x) y = "); | |
| 312 fsolve_fcn = extract_function | |
| 313 (c(0), "fsolve", fcn_name, fname, "; endfunction"); | |
| 314 } | |
| 315 | |
| 316 if (fsolve_fcn) | |
| 317 { | |
| 318 if (c(1).is_function_handle () || c(1).is_inline_function ()) | |
| 319 fsolve_jac = c(1).function_value (); | |
| 320 else | |
| 321 { | |
| 322 jac_name = unique_symbol_name ("__fsolve_jac__"); | |
| 323 jname = "function y = "; | |
| 324 jname.append (jac_name); | |
| 325 jname.append (" (x) jac = "); | |
| 326 fsolve_jac = extract_function | |
| 327 (c(1), "fsolve", jac_name, jname, "; endfunction"); | |
| 4140 | 328 |
| 5729 | 329 if (!fsolve_jac) |
| 330 { | |
| 331 if (fcn_name.length()) | |
| 332 clear_function (fcn_name); | |
| 333 fsolve_fcn = 0; | |
| 334 } | |
| 335 } | |
| 336 } | |
| 337 } | |
| 338 else | |
| 339 FSOLVE_ABORT1 ("incorrect number of elements in cell array"); | |
| 340 } | |
| 4140 | 341 |
| 5729 | 342 if (!fsolve_fcn && ! f_arg.is_cell()) |
| 343 { | |
| 344 if (f_arg.is_function_handle () || f_arg.is_inline_function ()) | |
| 345 fsolve_fcn = f_arg.function_value (); | |
| 346 else | |
| 347 { | |
| 348 switch (f_arg.rows ()) | |
| 349 { | |
| 350 case 1: | |
| 351 do | |
| 352 { | |
| 353 fcn_name = unique_symbol_name ("__fsolve_fcn__"); | |
| 354 fname = "function y = "; | |
| 355 fname.append (fcn_name); | |
| 356 fname.append (" (x) y = "); | |
| 357 fsolve_fcn = extract_function | |
| 358 (f_arg, "fsolve", fcn_name, fname, "; endfunction"); | |
| 359 } | |
| 360 while (0); | |
| 361 break; | |
| 4140 | 362 |
| 5729 | 363 case 2: |
| 4140 | 364 { |
| 5729 | 365 string_vector tmp = f_arg.all_strings (); |
| 366 | |
| 367 if (! error_state) | |
| 368 { | |
| 369 fcn_name = unique_symbol_name ("__fsolve_fcn__"); | |
| 370 fname = "function y = "; | |
| 371 fname.append (fcn_name); | |
| 372 fname.append (" (x) y = "); | |
| 373 fsolve_fcn = extract_function | |
| 374 (tmp(0), "fsolve", fcn_name, fname, "; endfunction"); | |
| 4140 | 375 |
| 5729 | 376 if (fsolve_fcn) |
| 377 { | |
| 378 jac_name = unique_symbol_name ("__fsolve_jac__"); | |
| 379 jname = "function y = "; | |
| 380 jname.append (jac_name); | |
| 381 jname.append (" (x) jac = "); | |
| 382 fsolve_jac = extract_function | |
| 383 (tmp(1), "fsolve", jac_name, jname, | |
| 384 "; endfunction"); | |
| 385 | |
| 386 if (!fsolve_jac) | |
| 387 { | |
| 388 if (fcn_name.length()) | |
| 389 clear_function (fcn_name); | |
| 390 fsolve_fcn = 0; | |
| 391 } | |
| 392 } | |
| 393 } | |
| 4140 | 394 } |
| 5729 | 395 } |
| 396 } | |
| 4140 | 397 } |
| 398 | |
| 399 if (error_state || ! fsolve_fcn) | |
| 3323 | 400 FSOLVE_ABORT (); |
| 2928 | 401 |
| 3419 | 402 ColumnVector x (args(1).vector_value ()); |
| 2928 | 403 |
| 3243 | 404 if (error_state) |
| 3323 | 405 FSOLVE_ABORT1 ("expecting vector as second argument"); |
| 3243 | 406 |
| 4132 | 407 if (nargin > 3) |
| 3243 | 408 warning ("fsolve: ignoring extra arguments"); |
| 2928 | 409 |
| 4132 | 410 if (nargout > 3) |
| 3243 | 411 warning ("fsolve: can't compute path output yet"); |
| 2928 | 412 |
| 3971 | 413 NLFunc nleqn_fcn (fsolve_user_function); |
| 4140 | 414 if (fsolve_jac) |
| 415 nleqn_fcn.set_jacobian_function (fsolve_user_jacobian); | |
| 416 | |
| 3971 | 417 NLEqn nleqn (x, nleqn_fcn); |
| 4122 | 418 nleqn.set_options (fsolve_opts); |
| 3243 | 419 |
| 5275 | 420 octave_idx_type info; |
| 3971 | 421 ColumnVector soln = nleqn.solve (info); |
| 2928 | 422 |
| 5729 | 423 if (fcn_name.length()) |
| 424 clear_function (fcn_name); | |
| 425 if (jac_name.length()) | |
| 426 clear_function (jac_name); | |
| 427 | |
| 3971 | 428 if (! error_state) |
| 429 { | |
| 7279 | 430 retval(2) = static_cast<double> (hybrd_info_to_fsolve_info (info)); |
| 431 retval(1) = nleqn.function_value (); | |
| 4132 | 432 retval(0) = soln; |
| 2928 | 433 |
| 4132 | 434 if (! nleqn.solution_ok () && nargout < 2) |
| 7279 | 435 { |
| 436 std::string msg = nleqn.error_message (); | |
| 437 error ("fsolve: %s", msg.c_str ()); | |
| 438 } | |
| 3971 | 439 } |
| 3243 | 440 } |
| 441 else | |
| 5823 | 442 print_usage (); |
| 2928 | 443 |
| 3243 | 444 unwind_protect::run_frame ("Ffsolve"); |
| 2928 | 445 |
| 446 return retval; | |
| 447 } | |
| 448 | |
| 449 /* | |
| 7279 | 450 %!function retval = f (p) |
| 451 %! x = p(1); | |
| 452 %! y = p(2); | |
| 453 %! z = p(3); | |
| 454 %! retval = zeros (3, 1); | |
| 455 %! retval(1) = sin(x) + y**2 + log(z) - 7; | |
| 456 %! retval(2) = 3*x + 2**y -z**3 + 1; | |
| 457 %! retval(3) = x + y + z - 5; | |
| 458 %!test | |
| 459 %! x_opt = [ 0.599054; | |
| 460 %! 2.395931; | |
| 461 %! 2.005014 ]; | |
| 462 %! tol = 1.0e-5; | |
| 463 %! [x, fval, info] = fsolve ("f", [ 0.5, 2.0, 2.5 ]); | |
| 464 %! info_bad = (info != 1); | |
| 465 %! solution_bad = sum (abs (x - x_opt) > tol); | |
| 466 %! value_bad = sum (abs (fval) > tol); | |
| 467 %! if (info_bad) | |
| 468 %! printf_assert ("info bad\n"); | |
| 469 %! else | |
| 470 %! printf_assert ("info good\n"); | |
| 471 %! endif | |
| 472 %! if (solution_bad) | |
| 473 %! printf_assert ("solution bad\n"); | |
| 474 %! else | |
| 475 %! printf_assert ("solution good\n"); | |
| 476 %! endif | |
| 477 %! if (value_bad) | |
| 478 %! printf_assert ("value bad\n"); | |
| 479 %! else | |
| 480 %! printf_assert ("value good\n"); | |
| 481 %! endif | |
| 482 %! assert(prog_output_assert("info good\nsolution good\nvalue good")); | |
| 483 | |
| 484 %!function retval = f (p) | |
| 485 %! x = p(1); | |
| 486 %! y = p(2); | |
| 487 %! z = p(3); | |
| 488 %! w = p(4); | |
| 489 %! retval = zeros (4, 1); | |
| 490 %! retval(1) = 3*x + 4*y + exp (z + w) - 1.007; | |
| 491 %! retval(2) = 6*x - 4*y + exp (3*z + w) - 11; | |
| 492 %! retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20; | |
| 493 %! retval(4) = x^2 + 2*y^3 + z - w - 4; | |
| 494 %!test | |
| 495 %! x_opt = [ -0.767297326653401; | |
| 496 %! 0.590671081117440; | |
| 497 %! 1.47190018629642; | |
| 498 %! -1.52719341133957 ]; | |
| 499 %! tol = 1.0e-5; | |
| 500 %! [x, fval, info] = fsolve ("f", [-1, 1, 2, -1]); | |
| 501 %! info_bad = (info != 1); | |
| 502 %! solution_bad = sum (abs (x - x_opt) > tol); | |
| 503 %! value_bad = sum (abs (fval) > tol); | |
| 504 %! if (info_bad) | |
| 505 %! printf_assert ("info bad\n"); | |
| 506 %! else | |
| 507 %! printf_assert ("info good\n"); | |
| 508 %! endif | |
| 509 %! if (solution_bad) | |
| 510 %! printf_assert ("solution bad\n"); | |
| 511 %! else | |
| 512 %! printf_assert ("solution good\n"); | |
| 513 %! endif | |
| 514 %! if (value_bad) | |
| 515 %! printf_assert ("value bad\n"); | |
| 516 %! else | |
| 517 %! printf_assert ("value good\n"); | |
| 518 %! endif | |
| 519 %! assert(prog_output_assert("info good\nsolution good\nvalue good")); | |
| 520 | |
| 521 %!test | |
| 522 %! fsolve_options ("tolerance", eps); | |
| 523 %! assert(fsolve_options ("tolerance") == eps); | |
| 524 | |
| 525 %!error <Invalid call to fsolve_options.*> fsolve_options ("foo", 1, 2); | |
| 526 */ | |
| 527 | |
| 528 /* | |
| 2928 | 529 ;;; Local Variables: *** |
| 530 ;;; mode: C++ *** | |
| 531 ;;; End: *** | |
| 532 */ |