|
 Re: fsolve
On Sun, May 25, 2008 at 8:57 PM, yakup murat mert < ymmert@...> wrote:
> function result=Uval(Us,Ud,Sbt,Uf,Vg)
> result=(0.0256*Us.+0.0256*Sbt.*(exp(Uf).*(exp(-Us).+Us.-1)+exp(-Uf).*(exp(Us.-Ud).-Us.-exp(-Ud))).^0.5).-Vg;
> ...
>
> ...
>
> Us=fsolve(@(Us) Uval(Us,Ud,Sbt,Uf,Vg) , init);
>
> ....
>
> ...
>
> Code above works on matlab. The equation is parametric. And depends on "Ud"
> "Sbt" "Uf" "Vg" parameters. "Us" is unknown (like x) "Ud" is 1*100 matrix
> and rest are constants. In matlab after this code "Us" becomes 1*100
> matrix and has been solved 100 times according to "Ud" .
>
> This code hasnt been accepted by octave. I checked everything but coulndt
> find equivalent...
>
> any suggestion....?
>
Currently, fsolve is implemented to use column vectors exclusively -
you can switch to that as a quick remedy. Matlab is probably more
smart here; it always shapes the
input to user function according to the shape of user's initial guess.
It even allows things like
fsolve (@fcn, ones(3,3,3)) (returns a 3x3x3-array)
The attached two changesets solve this. I noticed that
octave_value::vector_value with force_vector_conversion = true didn't
work on N-d arrays, and generally the XXX_vector_value methods were
sloppy, making unnecessary copies. Hence the first changeset.
The second changeset modifies the fsolve function to behave like
Matlab's. It can be applied without the first one, but then N-d arrays
still won't work.
make check seems to go without problems.
regards,
--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
[ov-xxx-value.diff]
# HG changeset patch
# User Jaroslav Hajek <highegg@...>
# Date 1211803315 -7200
# Node ID 936478aed051aafca8914e36afb89acf004d3807
# Parent edc25a3fb2bce80302402d876e660de57e332c43
simplify & cleanup octave_value::XXX_vector_value functions
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,17 @@
+2008-05-26 Jaroslav Hajek <highegg@...>
+
+ * ov.cc (make_vector_dims): New function.
+ (vector_value, complex_vector_value, float_vector_value,
+ float_complex_vector_value): Query N-d array values and simplify,
+ avoid copying.
+ (column_vector_value, complex_column_vector_value,
+ float_column_vector_value, float_complex_column_vector_value,
+ row_vector_value, complex_row_vector_value,
+ float_row_vector_value, float_complex_row_vector_value):
+ Simplify to trivial wrappers.
+ (int_vector_value): Avoid conversions if integer type, query N-d array
+ value, simplify.
+
2008-05-21 David Bateman <dbateman@...>
* DLD-FUNCTIONS/quad.cc (quad_float_user_function): New function.
diff --git a/src/ov.cc b/src/ov.cc
--- a/src/ov.cc
+++ b/src/ov.cc
@@ -1304,173 +1304,69 @@
return rep->list_value ();
}
-ColumnVector
-octave_value::column_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+static dim_vector
+make_vector_dims (const dim_vector& dv, bool force_vector_conversion,
+ const std::string& my_type, const std::string& wanted_type)
{
- ColumnVector retval;
+ dim_vector retval (dv);
+ retval.chop_trailing_singletons ();
+ octave_idx_type nel = dv.numel ();
- Matrix m = matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nc == 1)
+ if (retval.length () > 2 || (retval(0) != 1 && retval(1) != 1))
{
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real column vector");
+ if (!force_vector_conversion)
+ gripe_implicit_conversion ("Octave:array-as-vector",
+ my_type.c_str (), wanted_type.c_str ());
+ retval = dim_vector (nel);
}
return retval;
}
+ColumnVector
+octave_value::column_vector_value (bool force_string_conv,
+ bool frc_vec_conv) const
+{
+ return ColumnVector (vector_value (force_string_conv,
+ frc_vec_conv));
+}
+
ComplexColumnVector
octave_value::complex_column_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- ComplexColumnVector retval;
-
- ComplexMatrix m = complex_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex column vector");
- }
-
- return retval;
+ return ComplexColumnVector (complex_vector_value (force_string_conv,
+ frc_vec_conv));
}
RowVector
octave_value::row_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- RowVector retval;
-
- Matrix m = matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real row vector");
- }
-
- return retval;
+ return RowVector (vector_value (force_string_conv,
+ frc_vec_conv));
}
ComplexRowVector
octave_value::complex_row_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- ComplexRowVector retval;
-
- ComplexMatrix m = complex_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex row vector");
- }
-
- return retval;
+ return ComplexRowVector (complex_vector_value (force_string_conv,
+ frc_vec_conv));
}
-
-// Sloppy...
Array<double>
octave_value::vector_value (bool force_string_conv,
bool force_vector_conversion) const
{
- Array<double> retval;
-
- Matrix m = matrix_value (force_string_conv);
+ Array<double> retval = array_value (force_string_conv);
if (error_state)
return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
- else if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else if (nr > 0 && nc > 0)
- {
- if (! force_vector_conversion)
- gripe_implicit_conversion ("Octave:array-as-vector",
- type_name (), "real vector");
-
- retval.resize (nr * nc);
- octave_idx_type k = 0;
- for (octave_idx_type j = 0; j < nc; j++)
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
-
- retval (k++) = m (i, j);
- }
- }
else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real vector");
- }
-
- return retval;
+ return retval.reshape (make_vector_dims (retval.dims (),
+ force_vector_conversion,
+ type_name (), "real vector"));
}
Array<int>
@@ -1479,364 +1375,124 @@
{
Array<int> retval;
- Matrix m = matrix_value (force_string_conv);
+ if (is_integer_type ())
+ {
+ // query for the first type that is wide enough
+#if SIZEOF_INT == 2
+ retval = int16_array_value ();
+#elif SIZEOF_INT == 4
+ retval = int32_array_value ();
+#else
+ retval = int64_array_value ();
+#endif
+ }
+ else
+ {
+ const NDArray a = array_value (force_string_conv);
+ if (! error_state)
+ {
+ if (require_int)
+ {
+ retval.resize (a.dims ());
+ for (octave_idx_type i = 0; i < a.numel (); i++)
+ {
+ double ai = a.elem (i);
+ int v = static_cast<int> (ai);
+ if (ai == v)
+ retval.xelem (i) = v;
+ else
+ {
+ error ("conversion to integer value failed");
+ break;
+ }
+ }
+ }
+ else
+ retval = Array<int> (a);
+ }
+ }
+
if (error_state)
return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- {
- OCTAVE_QUIT;
-
- double d = m (0, i);
-
- if (require_int && D_NINT (d) != d)
- {
- error ("conversion to integer value failed");
- return retval;
- }
-
- retval (i) = static_cast<int> (d);
- }
- }
- else if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
-
- double d = m (i, 0);
-
- if (require_int && D_NINT (d) != d)
- {
- error ("conversion to integer value failed");
- return retval;
- }
-
- retval (i) = static_cast<int> (d);
- }
- }
- else if (nr > 0 && nc > 0)
- {
- if (! force_vector_conversion)
- gripe_implicit_conversion ("Octave:array-as-vector",
- type_name (), "real vector");
-
- retval.resize (nr * nc);
- octave_idx_type k = 0;
- for (octave_idx_type j = 0; j < nc; j++)
- {
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
-
- double d = m (i, j);
-
- if (require_int && D_NINT (d) != d)
- {
- error ("conversion to integer value failed");
- return retval;
- }
-
- retval (k++) = static_cast<int> (d);
- }
- }
- }
else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real vector");
- }
-
- return retval;
+ return retval.reshape (make_vector_dims (retval.dims (),
+ force_vector_conversion,
+ type_name (), "integer vector"));
}
Array<Complex>
octave_value::complex_vector_value (bool force_string_conv,
- bool force_vector_conversion) const
+ bool force_vector_conversion) const
{
- Array<Complex> retval;
-
- ComplexMatrix m = complex_matrix_value (force_string_conv);
+ Array<Complex> retval = complex_array_value (force_string_conv);
if (error_state)
return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- {
- OCTAVE_QUIT;
- retval (i) = m (0, i);
- }
- }
- else if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
- retval (i) = m (i, 0);
- }
- }
- else if (nr > 0 && nc > 0)
- {
- if (! force_vector_conversion)
- gripe_implicit_conversion ("Octave:array-as-vector",
- type_name (), "complex vector");
-
- retval.resize (nr * nc);
- octave_idx_type k = 0;
- for (octave_idx_type j = 0; j < nc; j++)
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
- retval (k++) = m (i, j);
- }
- }
else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex vector");
- }
-
- return retval;
+ return retval.reshape (make_vector_dims (retval.dims (),
+ force_vector_conversion,
+ type_name (), "complex vector"));
}
FloatColumnVector
octave_value::float_column_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- FloatColumnVector retval;
-
- FloatMatrix m = float_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real column vector");
- }
-
- return retval;
+ return FloatColumnVector (float_vector_value (force_string_conv,
+ frc_vec_conv));
}
FloatComplexColumnVector
octave_value::float_complex_column_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- FloatComplexColumnVector retval;
-
- FloatComplexMatrix m = float_complex_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex column vector");
- }
-
- return retval;
+ return FloatComplexColumnVector (float_complex_vector_value (force_string_conv,
+ frc_vec_conv));
}
FloatRowVector
octave_value::float_row_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- FloatRowVector retval;
-
- FloatMatrix m = float_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real row vector");
- }
-
- return retval;
+ return FloatRowVector (float_vector_value (force_string_conv,
+ frc_vec_conv));
}
FloatComplexRowVector
octave_value::float_complex_row_vector_value (bool force_string_conv,
- bool /* frc_vec_conv */) const
+ bool frc_vec_conv) const
{
- FloatComplexRowVector retval;
+ return FloatComplexRowVector (float_complex_vector_value (force_string_conv,
+ frc_vec_conv));
+}
- FloatComplexMatrix m = float_complex_matrix_value (force_string_conv);
+Array<float>
+octave_value::float_vector_value (bool force_string_conv,
+ bool force_vector_conversion) const
+{
+ Array<float> retval = float_array_value (force_string_conv);
if (error_state)
return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex row vector");
- }
-
- return retval;
-}
-
-// Sloppy...
-
-Array<float>
-octave_value::float_vector_value (bool force_string_conv,
- bool force_vector_conversion) const
-{
- Array<float> retval;
-
- FloatMatrix m = float_matrix_value (force_string_conv);
-
- if (error_state)
- return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- retval (i) = m (0, i);
- }
- else if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- retval (i) = m (i, 0);
- }
- else if (nr > 0 && nc > 0)
- {
- if (! force_vector_conversion)
- gripe_implicit_conversion ("Octave:array-as-vector",
- type_name (), "real vector");
-
- retval.resize (nr * nc);
- octave_idx_type k = 0;
- for (octave_idx_type j = 0; j < nc; j++)
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
-
- retval (k++) = m (i, j);
- }
- }
- else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "real vector");
- }
-
- return retval;
+ return retval.reshape (make_vector_dims (retval.dims (),
+ force_vector_conversion,
+ type_name (), "real vector"));
}
Array<FloatComplex>
octave_value::float_complex_vector_value (bool force_string_conv,
- bool force_vector_conversion) const
+ bool force_vector_conversion) const
{
- Array<FloatComplex> retval;
-
- FloatComplexMatrix m = float_complex_matrix_value (force_string_conv);
+ Array<FloatComplex> retval = float_complex_array_value (force_string_conv);
if (error_state)
return retval;
-
- octave_idx_type nr = m.rows ();
- octave_idx_type nc = m.columns ();
-
- if (nr == 1)
- {
- retval.resize (nc);
- for (octave_idx_type i = 0; i < nc; i++)
- {
- OCTAVE_QUIT;
- retval (i) = m (0, i);
- }
- }
- else if (nc == 1)
- {
- retval.resize (nr);
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
- retval (i) = m (i, 0);
- }
- }
- else if (nr > 0 && nc > 0)
- {
- if (! force_vector_conversion)
- gripe_implicit_conversion ("Octave:array-as-vector",
- type_name (), "complex vector");
-
- retval.resize (nr * nc);
- octave_idx_type k = 0;
- for (octave_idx_type j = 0; j < nc; j++)
- for (octave_idx_type i = 0; i < nr; i++)
- {
- OCTAVE_QUIT;
- retval (k++) = m (i, j);
- }
- }
else
- {
- std::string tn = type_name ();
- gripe_invalid_conversion (tn.c_str (), "complex vector");
- }
-
- return retval;
+ return retval.reshape (make_vector_dims (retval.dims (),
+ force_vector_conversion,
+ type_name (), "complex vector"));
}
int
[fsolve_xdims.diff]
# HG changeset patch
# User Jaroslav Hajek <highegg@...>
# Date 1211783291 -7200
# Node ID b282f098f4a83f9027ead3acf0bce422e9ec6932
# Parent 936478aed051aafca8914e36afb89acf004d3807
make fsolve pass arguments in a Matlab-compatible way
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -11,6 +11,14 @@
Simplify to trivial wrappers.
(int_vector_value): Avoid conversions if integer type, query N-d array
value, simplify.
+
+2008-05-26 Jaroslav Hajek <highegg@...>
+
+ * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function,
+ fsolve_user_jacobian): Reshape argument to original dims before
+ passing.
+ (Ffsolve): Save original dimensions of the starting guess and reshape
+ on return. Fix tests.
2008-05-21 David Bateman <dbateman@...>
diff --git a/src/DLD-FUNCTIONS/fsolve.cc b/src/DLD-FUNCTIONS/fsolve.cc
--- a/src/DLD-FUNCTIONS/fsolve.cc
+++ b/src/DLD-FUNCTIONS/fsolve.cc
@@ -31,6 +31,7 @@
#include <iostream>
#include <sstream>
+#include "dNDArray.h"
#include "NLEqn.h"
#include "defun-dld.h"
@@ -51,6 +52,9 @@
// Global pointer for optional user defined jacobian function.
static octave_function *fsolve_jac;
+
+// original dimensions of X0
+static dim_vector x_dims;
// Have we warned about imaginary values returned from user function?
static bool warned_fcn_imaginary = false;
@@ -110,9 +114,7 @@
if (n > 1)
{
- Matrix m (n, 1);
- for (octave_idx_type i = 0; i < n; i++)
- m (i, 0) = x (i);
+ NDArray m (ArrayN<double> (x, x_dims));
octave_value vars (m);
args(0) = vars;
}
@@ -135,7 +137,7 @@
warned_fcn_imaginary = true;
}
- retval = ColumnVector (tmp(0).vector_value ());
+ retval = ColumnVector (tmp(0).vector_value (false, true));
if (error_state || retval.length () <= 0)
gripe_user_supplied_eval ("fsolve");
@@ -161,9 +163,7 @@
if (n > 1)
{
- Matrix m (n, 1);
- for (octave_idx_type i = 0; i < n; i++)
- m(i,0) = x(i);
+ NDArray m (ArrayN<double> (x, x_dims));
octave_value vars (m);
args(0) = vars;
}
@@ -399,7 +399,9 @@
if (error_state || ! fsolve_fcn)
FSOLVE_ABORT ();
- ColumnVector x (args(1).vector_value ());
+ NDArray xa = args(1).array_value ();
+ x_dims = xa.dims ();
+ ColumnVector x (xa);
if (error_state)
FSOLVE_ABORT1 ("expecting vector as second argument");
@@ -429,7 +431,7 @@
{
retval(2) = static_cast<double> (hybrd_info_to_fsolve_info (info));
retval(1) = nleqn.function_value ();
- retval(0) = soln;
+ retval(0) = NDArray (ArrayN<double> (soln.reshape (x_dims)));
if (! nleqn.solution_ok () && nargout < 2)
{
@@ -460,7 +462,7 @@
%! 2.395931;
%! 2.005014 ];
%! tol = 1.0e-5;
-%! [x, fval, info] = fsolve ("f", [ 0.5, 2.0, 2.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);
@@ -492,10 +494,7 @@
%! 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 ];
+%! 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);
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave
|
