# HG changeset patch # User jwe # Date 1036124444 0 # Node ID 8c710385c5723d64e30bd23eea707952cc759784 # Parent 303b28a7a7e4d486bdbc9d311cd263b19e37aaa8 [project @ 2002-11-01 04:20:44 by jwe] diff -r 303b28a7a7e4 -r 8c710385c572 liboctave/ChangeLog --- a/liboctave/ChangeLog Fri Nov 01 02:53:14 2002 +0000 +++ b/liboctave/ChangeLog Fri Nov 01 04:20:44 2002 +0000 @@ -1,5 +1,10 @@ 2002-10-31 John W. Eaton + * ODESFunc.h (ODESFunc::ODES_fsub, ODESFunc::ODES_bsub, + ODESFunc::ODES_jsub): Reorder args for consistency with other + solvers. + * ODESSA.cc: Fix all callers. + * mx-inlines.cc (MX_BASE_REDUCTION_OP): Also return scalar MT_RESULT if nr == 1 && nc == 0 && dim == -1 (i.e., sum(zeros(1,0)) returns 0, not [](1x0)). diff -r 303b28a7a7e4 -r 8c710385c572 liboctave/ODESFunc.h --- a/liboctave/ODESFunc.h Fri Nov 01 02:53:14 2002 +0000 +++ b/liboctave/ODESFunc.h Fri Nov 01 04:20:44 2002 +0000 @@ -36,13 +36,13 @@ Matrix *dfdx; }; - typedef ColumnVector (*ODES_fsub) (double, const ColumnVector& x, + typedef ColumnVector (*ODES_fsub) (const ColumnVector& x, double, const ColumnVector& theta); - typedef ColumnVector (*ODES_bsub) (double, const ColumnVector& x, + typedef ColumnVector (*ODES_bsub) (const ColumnVector& x, double, const ColumnVector& theta, int column); - typedef Matrix (*ODES_jsub) (double, const ColumnVector& x, + typedef Matrix (*ODES_jsub) (const ColumnVector& x, double, const ColumnVector& theta); ODESFunc (void) diff -r 303b28a7a7e4 -r 8c710385c572 liboctave/ODESSA.cc --- a/liboctave/ODESSA.cc Fri Nov 01 02:53:14 2002 +0000 +++ b/liboctave/ODESSA.cc Fri Nov 01 04:20:44 2002 +0000 @@ -82,7 +82,7 @@ for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; - ColumnVector tmp_fval = user_fsub (t, tmp_state, tmp_param); + ColumnVector tmp_fval = user_fsub (tmp_state, t, tmp_param); // Return the function value as a C array object @@ -111,7 +111,7 @@ for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; - Matrix tmp_fval = user_jsub (t, tmp_state, tmp_param); + Matrix tmp_fval = user_jsub (tmp_state, t, tmp_param); for (int j = 0; j < n; j++) for (int i = 0; i < nrowpd; i++) @@ -137,7 +137,7 @@ for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; - ColumnVector tmp_fval = user_bsub (t, tmp_state, tmp_param, jpar); + ColumnVector tmp_fval = user_bsub (tmp_state, t, tmp_param, jpar); for (int i = 0; i < n; i++) dfdp[i] = tmp_fval(i); @@ -320,7 +320,7 @@ if (! sanity_checked) { - ColumnVector fval = user_fsub (t, x, theta); + ColumnVector fval = user_fsub (x, t, theta); if (fval.length () != x.length ()) { diff -r 303b28a7a7e4 -r 8c710385c572 src/ChangeLog --- a/src/ChangeLog Fri Nov 01 02:53:14 2002 +0000 +++ b/src/ChangeLog Fri Nov 01 04:20:44 2002 +0000 @@ -1,5 +1,10 @@ 2002-10-31 John W. Eaton + * DLD-FUNCTIONS/odessa.cc (odessa_user_f, odessa_user_j, + odessa_user_b): Reorder args for consistency with other solvers. + (Fodessa): Use extract_function to set args. + (odessa_user_j): Rename from odessa_user_mf. + * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function, fsolve_user_jacobian): Print warning if user returns complex value. * DLD-FUNCTIONS/quad.cc (quad_user_function): Likewise. diff -r 303b28a7a7e4 -r 8c710385c572 src/DLD-FUNCTIONS/odessa.cc --- a/src/DLD-FUNCTIONS/odessa.cc Fri Nov 01 02:53:14 2002 +0000 +++ b/src/DLD-FUNCTIONS/odessa.cc Fri Nov 01 04:20:44 2002 +0000 @@ -60,7 +60,7 @@ static int call_depth = 0; static ColumnVector -odessa_user_f (double t, const ColumnVector& x, const ColumnVector& theta) +odessa_user_f (const ColumnVector& x, double t, const ColumnVector& theta) { ColumnVector retval; @@ -76,14 +76,14 @@ else args(2) = Matrix (); + args(1) = t; + if (n > 1) - args(1) = x; + args(0) = x; else if (n == 1) - args(1) = x(0); + args(0) = x(0); else - args(1) = Matrix (); - - args(0) = t; + args(0) = Matrix (); if (odessa_f) { @@ -116,12 +116,11 @@ } static Matrix -odessa_user_mf (double t, const ColumnVector& x, const ColumnVector& theta, - octave_function *mf) +odessa_user_j (const ColumnVector& x, double t, const ColumnVector& theta) { Matrix retval; - if (mf) + if (odessa_j) { octave_value_list args; @@ -135,16 +134,16 @@ else args(2) = Matrix (); - if (n > 1) - args(1) = x; - else if (n == 1) - args(1) = x(0); - else - args(1) = Matrix (); - args(0) = t; - octave_value_list tmp = mf->do_multi_index_op (1, args); + if (n > 1) + args(0) = x; + else if (n == 1) + args(0) = x(0); + else + args(0) = Matrix (); + + octave_value_list tmp = odessa_j->do_multi_index_op (1, args); if (error_state) { @@ -172,14 +171,8 @@ return retval; } -static Matrix -odessa_user_j (double t, const ColumnVector& x, const ColumnVector& theta) -{ - return odessa_user_mf (t, x, theta, odessa_j); -} - static ColumnVector -odessa_user_b (double t, const ColumnVector& x, +odessa_user_b (const ColumnVector& x, double t, const ColumnVector& theta, int column) { ColumnVector retval; @@ -200,14 +193,14 @@ else args(2) = Matrix (); + args(1) = t; + if (n > 1) - args(1) = x; + args(0) = x; else if (n == 1) - args(1) = x(0); + args(0) = x(0); else - args(1) = Matrix (); - - args(0) = t; + args(0) = Matrix (); octave_value_list tmp = odessa_b->do_multi_index_op (1, args); @@ -331,12 +324,14 @@ The function must have the form\n\ \n\ @example\n\ -@var{xdot} = f (@var{x}, @var{t})\n\ +@var{xdot} = f (@var{x}, @var{t}, @var{p})\n\ @end example\n\ \n\ @noindent\n\ in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\ \n\ +The variable @var{p} is a vector of parameters.\n\ +\n\ The @var{fcn} argument may also be an array of strings\n\ \n\ @example\n\ @@ -351,10 +346,12 @@ The Jacobian function must have the form\n\ \n\ @example\n\ -@var{jac} = j (@var{x}, @var{t})\n\ +@var{jac} = j (@var{x}, @var{t}, @var{p})\n\ @end example\n\ \n\ -in which @var{jac} is the matrix of partial derivatives\n\ +in which @var{x}, @var{t}, and @var{p} have the same meanings as\n\ +above for the function @var{f}, and @var{jac} is the matrix of\n\ +partial derivatives\n\ @tex\n\ $$ J = {\\partial f_i \\over \\partial x_j} $$\n\ @end tex\n\ @@ -371,7 +368,18 @@ The function @var{b} must have the form\n\ \n\ @example\n\ -@var{dfdp} = b (@var{t}, @var{x}, @var{p})\n\ +@var{dfdp} = b (@var{x}, @var{t}, @var{p}, @var{c})\n\ +@end example\n\ +\n\ +in which @var{x}, @var{t}, and @var{p} have the same meanings as\n\ +above for the function @var{f}, @var{c} indicates which partial\n\ +derivatives to return in @var{dfdp}. For example, if @var{c} is 2,\n\ +you should return the vector\n\ +\n\ +@example\n\ + df_i\n\ +dfdp = ----, i = 1:length(x)\n\ + dp_2\n\ @end example\n\ \n\ The second argument, @var{x_0}, specifies the intial state of the system.\n\ @@ -426,47 +434,54 @@ octave_value f_arg = args(0); - if (f_arg.is_string ()) + int nr = f_arg.rows (); + + if (nr == 1) { - string_vector f_str_arg = f_arg.all_strings (); - - int len = f_str_arg.length (); - - if (len > 2) + odessa_f = extract_function + (f_arg, "odessa", "__odessa_fcn__", + "function xdot = __odessa_fcn__ (x, t, p) xdot = ", + "; endfunction"); + } + else if (nr == 2 || nr == 3) + { + string_vector tmp = f_arg.all_strings (); + + if (! error_state) { - std::string t = f_str_arg(2); + odessa_f = extract_function + (tmp(0), "odessa", "__odessa_fcn__", + "function xdot = __odessa_fcn__ (x, t, p) xdot = ", + "; endfunction"); - if (t.length () > 0) + if (odessa_f) { - odessa_b = is_valid_function (t, "odessa", true); - - if (! odessa_b) - ODESSA_ABORT1 - ("expecting b function name as argument 1"); + odessa_j = extract_function + (tmp(1), "odessa", "__odessa_jac__", + "function xdot = __odessa_jac__ (x, t, p) jac = ", + "; endfunction"); + + if (odessa_j && nr == 3) + { + odessa_b = extract_function + (tmp(1), "odessa", "__odessa_b__", + "function dfdp = __odessa_b__ (x, t, p, c) dfdp = ", + "; endfunction"); + + if (! odessa_b) + odessa_j = 0; + } + + if (! odessa_j) + odessa_f = 0; } } + } - if (len > 1) - { - std::string t = f_str_arg(1); - - if (t.length () > 0) - { - odessa_j = is_valid_function (t, "odessa", true); - - if (! odessa_j) - ODESSA_ABORT1 - ("expecting function name as argument 1"); - } - } + if (error_state || ! odessa_f) + ODESSA_ABORT1 + ("expecting function name as argument 1"); - if (len > 0) - odessa_f = is_valid_function (f_str_arg(0), "odessa", true); - - if (! odessa_f) - ODESSA_ABORT1 ("expecting function name as argument 1"); - } - ColumnVector state (args(1).vector_value ()); if (error_state)