home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / src / f-fsolve.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  6KB  |  304 lines
  1. // f-fsolve.cc                                           -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1993, 1994, 1995 John W. Eaton
  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
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  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
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <strstream.h>
  29.  
  30. #include "NLEqn.h"
  31.  
  32. #include "tree-const.h"
  33. #include "variables.h"
  34. #include "gripes.h"
  35. #include "error.h"
  36. #include "utils.h"
  37. #include "pager.h"
  38. #include "help.h"
  39. #include "defun-dld.h"
  40.  
  41. // Global pointer for user defined function required by hybrd1.
  42. static tree_fvc *fsolve_fcn;
  43.  
  44. static NLEqn_options fsolve_opts;
  45.  
  46. int
  47. hybrd_info_to_fsolve_info (int info)
  48. {
  49.   switch (info)
  50.     {
  51.     case -1:
  52.       info = -2;
  53.       break;
  54.  
  55.     case 0:
  56.       info = -1;
  57.       break;
  58.  
  59.     case 1:
  60.       break;
  61.  
  62.     case 2:
  63.       info = 4;
  64.       break;
  65.  
  66.     case 3:
  67.     case 4:
  68.     case 5:
  69.       info = 3;
  70.       break;
  71.  
  72.     default:
  73.       panic_impossible ();
  74.       break;
  75.     }
  76.  
  77.   return info;
  78. }
  79.  
  80. ColumnVector
  81. fsolve_user_function (const ColumnVector& x)
  82. {
  83.   ColumnVector retval;
  84.  
  85.   int n = x.capacity ();
  86.  
  87.   Octave_object args;
  88.   args.resize (1);
  89.  
  90.   if (n > 1)
  91.     {
  92.       Matrix m (n, 1);
  93.       for (int i = 0; i < n; i++)
  94.     m (i, 0) = x.elem (i);
  95.       tree_constant vars (m);
  96.       args(0) = vars;
  97.     }
  98.   else
  99.     {
  100.       double d = x.elem (0);
  101.       tree_constant vars (d);
  102.       args(0) = vars;
  103.     }
  104.  
  105.   if (fsolve_fcn)
  106.     {
  107.       Octave_object tmp = fsolve_fcn->eval (0, 1, args);
  108.       if (tmp.length () > 0 && tmp(0).is_defined ())
  109.     {
  110.       retval = tmp(0).vector_value ();
  111.  
  112.       if (error_state || retval.length () <= 0)
  113.         gripe_user_supplied_eval ("fsolve");
  114.     }
  115.       else
  116.     gripe_user_supplied_eval ("fsolve");
  117.     }
  118.  
  119.   return retval;
  120. }
  121.  
  122. DEFUN_DLD_BUILTIN ("fsolve", Ffsolve, Sfsolve, 2, 1,
  123.   "Solve nonlinear equations using Minpack.  Usage:\n\
  124. \n\
  125.   [X, INFO] = fsolve (F, X0)\n\
  126. \n\
  127. Where the first argument is the name of the  function to call to\n\
  128. compute the vector of function values.  It must have the form\n\
  129. \n\
  130.   y = f (x)
  131. \n\
  132. where y and x are vectors.")
  133. {
  134.   Octave_object retval;
  135.  
  136.   int nargin = args.length ();
  137.  
  138.   if (nargin != 2 || nargout > 3)
  139.     {
  140.       print_usage ("fsolve");
  141.       return retval;
  142.     }
  143.  
  144.   fsolve_fcn = is_valid_function (args(0), "fsolve", 1);
  145.   if (! fsolve_fcn || takes_correct_nargs (fsolve_fcn, 1, "fsolve", 1) != 1)
  146.     return retval;
  147.  
  148.   ColumnVector x = args(1).vector_value ();
  149.  
  150.   if (error_state)
  151.     {
  152.       error ("fsolve: expecting vector as second argument");
  153.       return retval;
  154.     }
  155.  
  156.   if (nargin > 2)
  157.     warning ("fsolve: ignoring extra arguments");
  158.  
  159.   if (nargout > 2)
  160.     warning ("fsolve: can't compute path output yet");
  161.  
  162.   NLFunc foo_fcn (fsolve_user_function);
  163.   NLEqn foo (x, foo_fcn);
  164.   foo.copy (fsolve_opts);
  165.  
  166.   int info;
  167.   ColumnVector soln = foo.solve (info);
  168.  
  169.   info = hybrd_info_to_fsolve_info (info);
  170.  
  171.   retval.resize (nargout ? nargout : 1);
  172.   retval(0) = soln, 1;
  173.  
  174.   if (nargout > 1)
  175.     retval(1) = (double) info;
  176.  
  177.   return retval;
  178. }
  179.  
  180. typedef void (NLEqn_options::*d_set_opt_mf) (double);
  181. typedef double (NLEqn_options::*d_get_opt_mf) (void);
  182.  
  183. #define MAX_TOKENS 1
  184.  
  185. struct NLEQN_OPTIONS
  186. {
  187.   const char *keyword;
  188.   const char *kw_tok[MAX_TOKENS + 1];
  189.   int min_len[MAX_TOKENS + 1];
  190.   int min_toks_to_match;
  191.   d_set_opt_mf d_set_fcn;
  192.   d_get_opt_mf d_get_fcn;
  193. };
  194.  
  195. static NLEQN_OPTIONS fsolve_option_table [] =
  196. {
  197.   { "tolerance",
  198.     { "tolerance", 0, },
  199.     { 1, 0, }, 1,
  200.     NLEqn_options::set_tolerance,
  201.     NLEqn_options::tolerance, },
  202.  
  203.   { 0,
  204.     { 0, 0, },
  205.     { 0, 0, }, 0,
  206.     0, 0, },
  207. };
  208.  
  209. static void
  210. print_fsolve_option_list (void)
  211. {
  212.   ostrstream output_buf;
  213.  
  214.   print_usage ("fsolve_options", 1);
  215.  
  216.   output_buf << "\n"
  217.          << "Options for fsolve include:\n\n"
  218.          << "  keyword                                  value\n"
  219.          << "  -------                                  -----\n\n";
  220.  
  221.   NLEQN_OPTIONS *list = fsolve_option_table;
  222.  
  223.   const char *keyword;
  224.   while ((keyword = list->keyword) != 0)
  225.     {
  226.       output_buf.form ("  %-40s ", keyword);
  227.  
  228.       double val = (fsolve_opts.*list->d_get_fcn) ();
  229.       if (val < 0.0)
  230.     output_buf << "computed automatically";
  231.       else
  232.     output_buf << val;
  233.  
  234.       output_buf << "\n";
  235.       list++;
  236.     }
  237.  
  238.   output_buf << "\n" << ends;
  239.   maybe_page_output (output_buf);
  240. }
  241.  
  242. static void
  243. do_fsolve_option (char *keyword, double val)
  244. {
  245.   NLEQN_OPTIONS *list = fsolve_option_table;
  246.  
  247.   while (list->keyword != 0)
  248.     {
  249.       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
  250.                 list->min_toks_to_match, MAX_TOKENS))
  251.     {
  252.       (fsolve_opts.*list->d_set_fcn) (val);
  253.  
  254.       return;
  255.     }
  256.       list++;
  257.     }
  258.  
  259.   warning ("fsolve_options: no match for `%s'", keyword);
  260. }
  261.  
  262. DEFUN_DLD_BUILTIN ("fsolve_options", Ffsolve_options, Sfsolve_options, -1, 1,
  263.   "fsolve_options (KEYWORD, VALUE)\n\
  264. \n\
  265. Set or show options for fsolve.  Keywords may be abbreviated\n\
  266. to the shortest match.")
  267. {
  268.   Octave_object retval;
  269.  
  270.   int nargin = args.length ();
  271.  
  272.   if (nargin == 0)
  273.     {
  274.       print_fsolve_option_list ();
  275.       return retval;
  276.     }
  277.   else if (nargin == 2)
  278.     {
  279.       char *keyword = args(0).string_value ();
  280.  
  281.       if (! error_state)
  282.     {
  283.       double val = args(1).double_value ();
  284.  
  285.       if (! error_state)
  286.         {
  287.           do_fsolve_option (keyword, val);
  288.           return retval;
  289.         }
  290.     }
  291.     }
  292.  
  293.   print_usage ("fsolve_options");
  294.  
  295.   return retval;
  296. }
  297.  
  298. /*
  299. ;;; Local Variables: ***
  300. ;;; mode: C++ ***
  301. ;;; page-delimiter: "^/\\*" ***
  302. ;;; End: ***
  303. */
  304.