fsolve compatibility with Matlab

[John W. Eaton]

On  9-Dec-2007, Sylvain Pelissier wrote:

| I found that there is a small difference of  fsolve between Octave and
| Matlab.
| In octave the syntax is :
| 
| [x,exitflag,output] = fsolve(fun,x0)

Actually, it is

  [x, info, msg]

| and in Matlab it's :
| 
| [x,fval,exitflag,output] = fsolve(fun,x0)
| 
| it could be better to add the output argument fval which is the value of
| the function at the solution x to have a better compatibility with matlab.

I checked in the following change.  Probably this will break some
existing code, but it is probably better to do that now, before we
release Octave 3.0, rather than later.

Comments?

jwe


liboctave/ChangeLog:

2007-12-10  John W. Eaton  <address@hidden>

        * NLEqn.h (NLEqn::fval): New data member.  Adjust constructors
        and assignment operator.
        (NLEqn::function_value): New function.
        * NLEqn.cc (NLEqn::solve): If solution is successful, compute
        function value.


src/ChangeLog:

2007-12-10  John W. Eaton  <address@hidden>

        * DLD-FUNCTIONS/fsolve.cc (Ffsolve):
        For compatibility, return [x, fval, info] instead of [x, info, msg].
        Move tests here from test/test_nonlin.m.


test/ChangeLog:

2007-12-10  John W. Eaton  <address@hidden>

        * test_nonlin.m: Delete.


Index: NEWS
===================================================================
RCS file: /cvs/octave/NEWS,v
retrieving revision 1.74
diff -u -u -r1.74 NEWS
--- NEWS        26 Nov 2007 20:42:09 -0000      1.74
+++ NEWS        10 Dec 2007 20:55:27 -0000
@@ -174,6 +175,17 @@
 
       [status, output] = system (cmd);
 
+
+ ** For compatibility with Matlab, the output of Octave's fsolve
+    function has been changed from
+
+      [x, info, msg] = fsolve (...);
+
+    to
+
+      [x, fval, info] = fsolve (...);
+
+
  ** For compatibility with Matlab, normcdf, norminv, normpdf, and
     normrnd have been modified to compute distributions using the
     standard deviation instead of the variance.
Index: liboctave/NLEqn.cc
===================================================================
RCS file: /cvs/octave/liboctave/NLEqn.cc,v
retrieving revision 1.41
diff -u -u -r1.41 NLEqn.cc
--- liboctave/NLEqn.cc  12 Oct 2007 21:27:14 -0000      1.41
+++ liboctave/NLEqn.cc  10 Dec 2007 20:55:28 -0000
@@ -178,6 +178,8 @@
 
       if (f77_exception_encountered)
        (*current_liboctave_error_handler) ("unrecoverable error in hybrj1");
+      else
+       fval = ColumnVector (fvec);
     }
   else
     {
@@ -195,6 +197,8 @@
 
       if (f77_exception_encountered)
        (*current_liboctave_error_handler) ("unrecoverable error in hybrd1");
+      else
+       fval = ColumnVector (fvec);
     }
 
   return retval;
Index: liboctave/NLEqn.h
===================================================================
RCS file: /cvs/octave/liboctave/NLEqn.h,v
retrieving revision 1.30
diff -u -u -r1.30 NLEqn.h
--- liboctave/NLEqn.h   30 Nov 2007 18:53:30 -0000      1.30
+++ liboctave/NLEqn.h   10 Dec 2007 20:55:28 -0000
@@ -27,6 +27,7 @@
 #include <cfloat>
 
 #include "NLEqn-opts.h"
+#include "lo-ieee.h"
 #include "lo-math.h"
 
 class
@@ -36,13 +37,15 @@
 public:
 
   NLEqn (void)
-    : NLFunc (), NLEqn_options (), x (), solution_status (0) { }
+    : NLFunc (), NLEqn_options (), x (), fval (),
+      solution_status (0) { }
 
   NLEqn (const ColumnVector& xx, const NLFunc f) 
-    : NLFunc (f), NLEqn_options (), x (xx), solution_status (0) { }
+    : NLFunc (f), NLEqn_options (), x (xx), fval (x.numel (), octave_NaN),
+      solution_status (0) { }
 
   NLEqn (const NLEqn& a)
-    : NLFunc (a.fun, a.jac), NLEqn_options (), x (a.x),
+    : NLFunc (a.fun, a.jac), NLEqn_options (), x (a.x), fval (a.fval),
       solution_status (a.solution_status) { }
 
   NLEqn& operator = (const NLEqn& a)
@@ -53,6 +56,7 @@
          NLEqn_options::operator = (a);
 
          x = a.x;
+         fval = a.fval;
          solution_status = a.solution_status;
        }
       return *this;
@@ -91,11 +95,14 @@
 
   bool solution_ok (void) const { return solution_status == 1; }
 
+  ColumnVector function_value (void) const { return fval; }
+
   std::string error_message (void) const;
 
 private:
 
   ColumnVector x;
+  ColumnVector fval;
   octave_idx_type solution_status;
 
   void error (const char* msg);
Index: src/DLD-FUNCTIONS/fsolve.cc
===================================================================
RCS file: /cvs/octave/src/DLD-FUNCTIONS/fsolve.cc,v
retrieving revision 1.33
diff -u -u -r1.33 fsolve.cc
--- src/DLD-FUNCTIONS/fsolve.cc 12 Oct 2007 21:27:35 -0000      1.33
+++ src/DLD-FUNCTIONS/fsolve.cc 10 Dec 2007 20:55:30 -0000
@@ -223,7 +223,7 @@
 
 DEFUN_DLD (fsolve, args, nargout,
   "-*- texinfo -*-\n\
address@hidden {Loadable Function} address@hidden, @var{info}, @var{msg}] =} 
fsolve (@var{fcn}, @var{x0})\n\
address@hidden {Loadable Function} address@hidden, @var{fval}, @var{info}] =} 
fsolve (@var{fcn}, @var{x0})\n\
 Given @var{fcn}, the name of a function of the form @code{f (@var{x})}\n\
 and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\
 equations such that @code{f(@var{x}) == 0}.\n\
@@ -405,15 +405,15 @@
 
       if (! error_state)
        {
-         std::string msg = nleqn.error_message ();
-
-         retval(2) = msg;
-         retval(1) = static_cast<double> (hybrd_info_to_fsolve_info (info));
-
+         retval(2) = static_cast<double> (hybrd_info_to_fsolve_info (info));
+         retval(1) = nleqn.function_value ();
          retval(0) = soln;
 
          if (! nleqn.solution_ok () && nargout < 2)
-           error ("fsolve: %s", msg.c_str ());
+           {
+             std::string msg = nleqn.error_message ();
+             error ("fsolve: %s", msg.c_str ());
+           }
        }
     }
   else
@@ -425,6 +425,85 @@
 }
 
 /*
+%!function retval = f (p) 
+%!  x = p(1);
+%!  y = p(2);
+%!  z = p(3);
+%!  retval = zeros (3, 1);
+%!  retval(1) = sin(x) + y**2 + log(z) - 7;
+%!  retval(2) = 3*x + 2**y -z**3 + 1;
+%!  retval(3) = x + y + z - 5;
+%!test
+%! x_opt = [ 0.599054;
+%! 2.395931;
+%! 2.005014 ];
+%! tol = 1.0e-5;
+%! [x, fval, info] = fsolve ("f", [ 0.5, 2.0, 2.5 ]);
+%! info_bad = (info != 1);
+%! solution_bad = sum (abs (x - x_opt) > tol);
+%! value_bad = sum (abs (fval) > tol);
+%! if (info_bad)
+%!   printf_assert ("info bad\n");
+%! else
+%!   printf_assert ("info good\n");
+%! endif
+%! if (solution_bad)
+%!   printf_assert ("solution bad\n");
+%! else
+%!   printf_assert ("solution good\n");
+%! endif
+%! if (value_bad)
+%!   printf_assert ("value bad\n");
+%! else
+%!   printf_assert ("value good\n");
+%! endif
+%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
+
+%!function retval = f (p)
+%!  x = p(1);
+%!  y = p(2);
+%!  z = p(3);
+%!  w = p(4);
+%!  retval = zeros (4, 1);
+%!  retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
+%!  retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
+%!  retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
+%!  retval(4) = x^2 + 2*y^3 + z - w - 4;
+%!test
+%! x_opt = [ -0.767297326653401;
+%! 0.590671081117440;
+%! 1.47190018629642;
+%! -1.52719341133957 ];
+%! tol = 1.0e-5;
+%! [x, fval, info] = fsolve ("f", [-1, 1, 2, -1]);
+%! info_bad = (info != 1);
+%! solution_bad = sum (abs (x - x_opt) > tol);
+%! value_bad = sum (abs (fval) > tol);
+%! if (info_bad)
+%!   printf_assert ("info bad\n");
+%! else
+%!   printf_assert ("info good\n");
+%! endif
+%! if (solution_bad)
+%!   printf_assert ("solution bad\n");
+%! else
+%!   printf_assert ("solution good\n");
+%! endif
+%! if (value_bad)
+%!   printf_assert ("value bad\n");
+%! else
+%!   printf_assert ("value good\n");
+%! endif
+%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
+
+%!test
+%! fsolve_options ("tolerance", eps);
+%! assert(fsolve_options ("tolerance") == eps);
+
+%!error <Invalid call to fsolve_options.*> fsolve_options ("foo", 1, 2);
+*/
+
+/*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
 ;;; End: ***
Index: test/test_nonlin.m
===================================================================
RCS file: test/test_nonlin.m
diff -N test/test_nonlin.m
--- test/test_nonlin.m  12 Oct 2007 21:27:37 -0000      1.5
+++ /dev/null   1 Jan 1970 00:00:00 -0000
@@ -1,103 +0,0 @@
-## Copyright (C) 2006, 2007 John W. Eaton
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-%% Automatically generated from DejaGNU files
-
-%% test/octave.test/nonlin/fsolve-1.m
-%!function retval = f (p) 
-%!  x = p(1);
-%!  y = p(2);
-%!  z = p(3);
-%!  retval = zeros (3, 1);
-%!  retval(1) = sin(x) + y**2 + log(z) - 7;
-%!  retval(2) = 3*x + 2**y -z**3 + 1;
-%!  retval(3) = x + y + z - 5;
-%!test
-%! x_opt = [ 0.599054;
-%! 2.395931;
-%! 2.005014 ];
-%! tol = 1.0e-5;
-%! [x, info] = fsolve ("f", [ 0.5, 2.0, 2.5 ]);
-%! val = f (x);
-%! info_bad = (info != 1);
-%! solution_bad = sum (abs (x - x_opt) > tol);
-%! value_bad = sum (abs (val) > tol);
-%! if (info_bad)
-%!   printf_assert ("info bad\n");
-%! else
-%!   printf_assert ("info good\n");
-%! endif
-%! if (solution_bad)
-%!   printf_assert ("solution bad\n");
-%! else
-%!   printf_assert ("solution good\n");
-%! endif
-%! if (value_bad)
-%!   printf_assert ("value bad\n");
-%! else
-%!   printf_assert ("value good\n");
-%! endif
-%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
-
-%% test/octave.test/nonlin/fsolve-2.m
-%!function retval = f (p)
-%!  x = p(1);
-%!  y = p(2);
-%!  z = p(3);
-%!  w = p(4);
-%!  retval = zeros (4, 1);
-%!  retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
-%!  retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
-%!  retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
-%!  retval(4) = x^2 + 2*y^3 + z - w - 4;
-%!test
-%! x_opt = [ -0.767297326653401;
-%! 0.590671081117440;
-%! 1.47190018629642;
-%! -1.52719341133957 ];
-%! tol = 1.0e-5;
-%! [x, info] = fsolve ("f", [-1, 1, 2, -1]);
-%! val = f (x);
-%! info_bad = (info != 1);
-%! solution_bad = sum (abs (x - x_opt) > tol);
-%! value_bad = sum (abs (val) > tol);
-%! if (info_bad)
-%!   printf_assert ("info bad\n");
-%! else
-%!   printf_assert ("info good\n");
-%! endif
-%! if (solution_bad)
-%!   printf_assert ("solution bad\n");
-%! else
-%!   printf_assert ("solution good\n");
-%! endif
-%! if (value_bad)
-%!   printf_assert ("value bad\n");
-%! else
-%!   printf_assert ("value good\n");
-%! endif
-%! assert(prog_output_assert("info good\nsolution good\nvalue good"));
-
-%% test/octave.test/nonlin/fsolve_options-1.m
-%!test
-%! fsolve_options ("tolerance", eps);
-%! assert(fsolve_options ("tolerance") == eps);
-
-%% test/octave.test/nonlin/fsolve_options-3.m
-%!error <Invalid call to fsolve_options.*> fsolve_options ("foo", 1, 2);
-