*** ./liboctave/NLEqn-opts.in.orig 2004-08-06 12:09:55.000000000 +0200 --- ./liboctave/NLEqn-opts.in 2004-08-06 12:11:12.000000000 +0200 *************** *** 23,25 **** --- 23,37 ---- INIT_VALUE = "::sqrt (DBL_EPSILON)" SET_EXPR = "(val > 0.0) ? val : ::sqrt (DBL_EPSILON)" END_OPTION + + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{fsolve}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/DASPK-opts.in.orig 2004-08-06 12:09:55.000000000 +0200 --- ./liboctave/DASPK-opts.in 2004-08-06 12:11:12.000000000 +0200 *************** *** 296,298 **** --- 296,310 ---- INIT_VALUE = "-1.0" SET_EXPR = "(val >= 0.0) ? val : -1.0" END_OPTION + + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{daspk}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/DAEFunc.h.orig 2004-08-06 12:09:55.000000000 +0200 --- ./liboctave/DAEFunc.h 2004-08-24 11:43:06.000000000 +0200 *************** *** 43,59 **** const ColumnVector& xdot, double t, double cj); DAEFunc (void) ! : fun (0), jac (0), reset (true) { } DAEFunc (DAERHSFunc f) ! : fun (f), jac (0), reset (true) { } DAEFunc (DAERHSFunc f, DAEJacFunc j) ! : fun (f), jac (j), reset (true) { } DAEFunc (const DAEFunc& a) ! : fun (a.fun), jac (a.jac), reset (a.reset) { } DAEFunc& operator = (const DAEFunc& a) { --- 43,68 ---- const ColumnVector& xdot, double t, double cj); + typedef ColumnVector (*DAERHSFuncCached) (const ColumnVector& x, + const ColumnVector& xdot, + double t, int& ires, double cj); + + typedef Matrix (*DAEJacFuncCached) (const ColumnVector& x, + const ColumnVector& xdot, + double t, int& ires, double cj); + DAEFunc (void) ! : fun (0), jac (0), fun_cached (0), jac_cached (0), reset (true) { } DAEFunc (DAERHSFunc f) ! : fun (f), jac (0), fun_cached (0), jac_cached (0), reset (true) { } DAEFunc (DAERHSFunc f, DAEJacFunc j) ! : fun (f), jac (j), fun_cached (0), jac_cached (0), reset (true) { } DAEFunc (const DAEFunc& a) ! : fun (a.fun), jac (a.jac), fun_cached (a.fun_cached), ! jac_cached (a.jac_cached), reset (a.reset) { } DAEFunc& operator = (const DAEFunc& a) { *************** *** 70,75 **** --- 79,86 ---- DAERHSFunc function (void) const { return fun; } + DAERHSFuncCached function_cached (void) const { return fun_cached; } + DAEFunc& set_function (DAERHSFunc f) { fun = f; *************** *** 77,84 **** --- 88,104 ---- return *this; } + DAEFunc& set_function_cached (DAERHSFuncCached f) + { + fun_cached = f; + reset = true; + return *this; + } + DAEJacFunc jacobian_function (void) const { return jac; } + DAEJacFuncCached jacobian_function_cached (void) const { return jac_cached; } + DAEFunc& set_jacobian_function (DAEJacFunc j) { jac = j; *************** *** 86,95 **** --- 106,124 ---- return *this; } + DAEFunc& set_jacobian_function_cached (DAEJacFuncCached j) + { + jac_cached = j; + reset = true; + return *this; + } + protected: DAERHSFunc fun; DAEJacFunc jac; + DAERHSFuncCached fun_cached; + DAEJacFuncCached jac_cached; // This variable is TRUE when this object is constructed, and also // after any internal data has changed. Derived classes may use *** ./liboctave/DASPK.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./liboctave/DASPK.cc 2004-08-06 14:59:12.000000000 +0200 *************** *** 65,75 **** static DAEFunc::DAERHSFunc user_fun; static DAEFunc::DAEJacFunc user_jac; static int nn; static int ddaspk_f (const double& time, const double *state, const double *deriv, ! const double&, double *delta, int& ires, double *, int *) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; --- 65,77 ---- static DAEFunc::DAERHSFunc user_fun; static DAEFunc::DAEJacFunc user_jac; + static DAEFunc::DAERHSFuncCached user_fun_cached; + static DAEFunc::DAEJacFuncCached user_jac_cached; static int nn; static int ddaspk_f (const double& time, const double *state, const double *deriv, ! const double& cj, double *delta, int& ires, double *, int *) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; *************** *** 83,89 **** tmp_state.elem (i) = state [i]; } ! tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires); if (ires >= 0) { --- 85,94 ---- tmp_state.elem (i) = state [i]; } ! if (user_fun_cached) ! tmp_delta = user_fun_cached (tmp_state, tmp_deriv, time, ires, cj); ! else ! tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires); if (ires >= 0) { *************** *** 137,143 **** tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) --- 142,155 ---- tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd; ! if (user_jac_cached) ! { ! int dummy; ! tmp_pd = user_jac_cached (tmp_state, tmp_deriv, time, dummy, cj); ! } ! else ! tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) *************** *** 188,195 **** // DAEFunc ! user_fun = DAEFunc::function (); ! user_jac = DAEFunc::jacobian_function (); if (user_fun) { --- 200,219 ---- // DAEFunc ! if (jacobian ()) ! { ! user_fun_cached = DAEFunc::function_cached (); ! user_jac_cached = DAEFunc::jacobian_function_cached (); ! user_fun = 0; ! user_jac = 0; ! } ! else ! { ! user_fun_cached = 0; ! user_jac_cached = 0; ! user_fun = DAEFunc::function (); ! user_jac = DAEFunc::jacobian_function (); ! } if (user_fun) { *************** *** 206,211 **** --- 230,250 ---- return retval; } } + else if (user_fun_cached) + { + int ires = 0; + double dummy = 0; + ColumnVector res = (*user_fun_cached) (x, xdot, t, ires, dummy); + + if (res.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("daspk: inconsistent sizes for state and residual vectors"); + + integration_error = true; + return retval; + } + } else { (*current_liboctave_error_handler) *************** *** 215,221 **** return retval; } ! info(4) = user_jac ? 1 : 0; DAEFunc::reset = false; --- 254,260 ---- return retval; } ! info(4) = (user_jac || user_jac_cached) ? 1 : 0; DAEFunc::reset = false; *** ./liboctave/ODESFunc.h.orig 2004-08-24 10:54:09.000000000 +0200 --- ./liboctave/ODESFunc.h 2004-08-24 11:15:44.000000000 +0200 *************** *** 45,64 **** typedef Matrix (*ODES_jsub) (const ColumnVector& x, double, const ColumnVector& theta); ODESFunc (void) ! : fsub (0), bsub (0), jsub (0) { } ODESFunc (ODES_fsub f) ! : fsub (f), bsub (0), jsub (0) { } ODESFunc (ODES_fsub f, ODES_bsub b) ! : fsub (f), bsub (b), jsub (0) { } ODESFunc (ODES_fsub f, ODES_bsub b, ODES_jsub j) ! : fsub (f), bsub (b), jsub (j) { } ODESFunc (const ODESFunc& a) ! : fsub (a.fsub), bsub (a.bsub), jsub (a.jsub) { } ODESFunc& operator = (const ODESFunc& a) { --- 45,82 ---- typedef Matrix (*ODES_jsub) (const ColumnVector& x, double, const ColumnVector& theta); + typedef ColumnVector (*ODES_fsub_cached) (const ColumnVector& x, double, + const ColumnVector& theta, + int column); + + typedef ColumnVector (*ODES_bsub_cached) (const ColumnVector& x, double, + const ColumnVector& theta, + int column); + + typedef Matrix (*ODES_jsub_cached) (const ColumnVector& x, double, + const ColumnVector& theta, + int column); + ODESFunc (void) ! : fsub (0), bsub (0), jsub (0), fsub_cached (0), bsub_cached (0), ! jsub_cached (0) { } ODESFunc (ODES_fsub f) ! : fsub (f), bsub (0), jsub (0), fsub_cached (0), bsub_cached (0), ! jsub_cached (0) { } ODESFunc (ODES_fsub f, ODES_bsub b) ! : fsub (f), bsub (b), jsub (0), fsub_cached (0), bsub_cached (0), ! jsub_cached (0) { } ODESFunc (ODES_fsub f, ODES_bsub b, ODES_jsub j) ! : fsub (f), bsub (b), jsub (j), fsub_cached (0), bsub_cached (0), ! jsub_cached (0) { } ODESFunc (const ODESFunc& a) ! : fsub (a.fsub), bsub (a.bsub), jsub (a.jsub), ! fsub_cached (a.fsub_cached), bsub_cached (a.bsub_cached), ! jsub_cached (a.jsub_cached) { } ODESFunc& operator = (const ODESFunc& a) { *************** *** 97,107 **** --- 115,152 ---- return *this; } + ODES_fsub_cached fsub_function_cached (void) const { return fsub_cached; } + + ODESFunc& set_fsub_function_cached (ODES_fsub_cached f) + { + fsub_cached = f; + return *this; + } + + ODES_bsub_cached bsub_function_cached (void) const { return bsub_cached; } + + ODESFunc& set_bsub_function_cached (ODES_bsub_cached b) + { + bsub_cached = b; + return *this; + } + + ODES_jsub_cached jsub_function_cached (void) const { return jsub_cached; } + + ODESFunc& set_jsub_function_cached (ODES_jsub_cached j) + { + jsub_cached = j; + return *this; + } + protected: ODES_fsub fsub; ODES_bsub bsub; ODES_jsub jsub; + ODES_fsub_cached fsub_cached; + ODES_bsub_cached bsub_cached; + ODES_jsub_cached jsub_cached; }; #endif *** ./liboctave/DASRT-opts.in.orig 2004-08-06 14:32:45.000000000 +0200 --- ./liboctave/DASRT-opts.in 2004-08-06 14:33:32.000000000 +0200 *************** *** 104,106 **** --- 104,118 ---- INIT_VALUE = "-1" SET_EXPR = "(val >= 0) ? val : -1" END_OPTION + + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{dasrt}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/DASRT.cc.orig 2004-08-06 19:16:38.000000000 +0200 --- ./liboctave/DASRT.cc 2004-08-08 18:37:49.000000000 +0200 *************** *** 36,41 **** --- 36,42 ---- #include "lo-error.h" #include "lo-sstream.h" #include "quit.h" + #include "lo-ieee.h" typedef int (*dasrt_fcn_ptr) (const double&, const double*, const double*, double*, int&, double*, int*); *************** *** 59,64 **** --- 60,67 ---- static DAEFunc::DAERHSFunc user_fsub; static DAEFunc::DAEJacFunc user_jsub; + static DAEFunc::DAERHSFuncCached user_fsub_cached; + static DAEFunc::DAEJacFuncCached user_jsub_cached; static DAERTFunc::DAERTConstrFunc user_csub; static int nn; *************** *** 78,84 **** tmp_deriv(i) = deriv[i]; } ! ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); if (tmp_fval.length () == 0) ires = -2; --- 81,94 ---- tmp_deriv(i) = deriv[i]; } ! ColumnVector tmp_fval; ! if (user_fsub_cached) ! // XXX FIXME XXX The underlying fortran doesn't give us the value of ! // cj!! Set it to an illegal value, so that the jacobian in the cache ! // will be ignored ! tmp_fval = (*user_fsub_cached) (tmp_state, tmp_deriv, t, ires, octave_Inf); ! else ! tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); if (tmp_fval.length () == 0) ires = -2; *************** *** 110,116 **** tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) --- 120,132 ---- tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd; ! int ires; ! ! if (user_jsub_cached) ! tmp_pd = (*user_jsub_cached) (tmp_state, tmp_deriv, time, ires, cj); ! else ! tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) *************** *** 225,232 **** // DAEFunc ! user_fsub = DAEFunc::function (); ! user_jsub = DAEFunc::jacobian_function (); if (user_fsub) { --- 241,260 ---- // DAEFunc ! if (jacobian ()) ! { ! user_fsub_cached = DAEFunc::function_cached (); ! user_jsub_cached = DAEFunc::jacobian_function_cached (); ! user_fsub = 0; ! user_jsub = 0; ! } ! else ! { ! user_fsub = DAEFunc::function (); ! user_jsub = DAEFunc::jacobian_function (); ! user_fsub_cached = 0; ! user_jsub_cached = 0; ! } if (user_fsub) { *************** *** 243,248 **** --- 271,292 ---- return; } } + else if (user_fsub_cached) + { + int ires = 0; + double dummy = 0; + + ColumnVector fval = (*user_fsub_cached) (x, xdot, t, ires, dummy); + + if (fval.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("dasrt: inconsistent sizes for state and residual vectors"); + + integration_error = true; + return; + } + } else { (*current_liboctave_error_handler) *************** *** 252,258 **** return; } ! info(4) = user_jsub ? 1 : 0; DAEFunc::reset = false; --- 296,302 ---- return; } ! info(4) = (user_jsub || user_jsub_cached) ? 1 : 0; DAEFunc::reset = false; *** ./liboctave/DASSL-opts.in.orig 2004-08-06 14:34:07.000000000 +0200 --- ./liboctave/DASSL-opts.in 2004-08-06 14:34:50.000000000 +0200 *************** *** 131,133 **** --- 131,145 ---- INIT_VALUE = "-1" SET_EXPR = "(val >= 0) ? val : -1" END_OPTION + + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{dassl}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/LSODE-opts.in.orig 2004-08-06 14:34:35.000000000 +0200 --- ./liboctave/LSODE-opts.in 2004-08-06 14:35:01.000000000 +0200 *************** *** 128,130 **** --- 128,142 ---- INIT_VALUE = "100000" SET_EXPR = "val" END_OPTION + + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{lsode}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/ODESSA.cc.orig 2004-08-24 11:06:26.000000000 +0200 --- ./liboctave/ODESSA.cc 2004-08-24 11:29:30.000000000 +0200 *************** *** 67,73 **** static ODESFunc::ODES_fsub user_fsub; static ODESFunc::ODES_bsub user_bsub; static ODESFunc::ODES_jsub user_jsub; ! static int odessa_f (int* neq, const double& t, double *state, --- 67,75 ---- static ODESFunc::ODES_fsub user_fsub; static ODESFunc::ODES_bsub user_bsub; static ODESFunc::ODES_jsub user_jsub; ! static ODESFunc::ODES_fsub_cached user_fsub_cached; ! static ODESFunc::ODES_bsub_cached user_bsub_cached; ! static ODESFunc::ODES_jsub_cached user_jsub_cached; static int odessa_f (int* neq, const double& t, double *state, *************** *** 88,94 **** for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! ColumnVector tmp_fval = user_fsub (tmp_state, t, tmp_param); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); --- 90,103 ---- for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! ColumnVector tmp_fval; ! if (user_fsub_cached) ! // XXX FIXME XXX The underlying fortran doesn't give us the value of ! // jpar!! Set it to an illegal value, so that the inhomogeneity matrix ! // in the cache will be ignored ! tmp_fval = user_fsub_cached (tmp_state, t, tmp_param, -1); ! else ! tmp_fval = user_fsub (tmp_state, t, tmp_param); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); *************** *** 122,128 **** for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! Matrix tmp_fval = user_jsub (tmp_state, t, tmp_param); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); --- 131,144 ---- for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! Matrix tmp_fval; ! if (user_jsub_cached) ! // XXX FIXME XXX The underlying fortran doesn't give us the value of ! // jpar!! Set it to an illegal value, so that the inhomogeneity matrix ! // in the cache will be ignored ! tmp_fval = user_jsub_cached (tmp_state, t, tmp_param, -1); ! else ! tmp_fval = user_jsub (tmp_state, t, tmp_param); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); *************** *** 157,163 **** for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! ColumnVector tmp_fval = user_bsub (tmp_state, t, tmp_param, jpar); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); --- 173,183 ---- for (int i = 0; i < n_par; i++) tmp_param(i) = par[i]; ! ColumnVector tmp_fval; ! if (user_bsub_cached) ! tmp_fval = user_bsub_cached (tmp_state, t, tmp_param, jpar); ! else ! tmp_fval = user_bsub (tmp_state, t, tmp_param, jpar); if (tmp_fval.length () == 0) octave_jump_to_enclosing_context (); *************** *** 269,281 **** integration_error = false; ! user_fsub = ODESFunc::fsub_function (); ! user_bsub = ODESFunc::bsub_function (); ! user_jsub = ODESFunc::jsub_function (); ! int idf; ! if (user_bsub) idf = 1; else idf = 0; --- 289,316 ---- integration_error = false; ! if (jacobian ()) ! { ! user_fsub_cached = ODESFunc::fsub_function_cached (); ! user_jsub_cached = ODESFunc::jsub_function_cached (); ! user_bsub_cached = ODESFunc::bsub_function_cached (); ! user_fsub = 0; ! user_bsub = 0; ! user_jsub = 0; ! } ! else ! { ! user_fsub = ODESFunc::fsub_function (); ! user_bsub = ODESFunc::bsub_function (); ! user_jsub = ODESFunc::jsub_function (); ! user_fsub_cached = 0; ! user_bsub_cached = 0; ! user_jsub_cached = 0; ! } ! int idf; ! if (user_bsub || user_bsub_cached) idf = 1; else idf = 0; *************** *** 292,298 **** if (integration_method () == "stiff") { ! if (user_jsub) method_flag = 21; else method_flag = 22; --- 327,333 ---- if (integration_method () == "stiff") { ! if (user_jsub || user_jsub_cached) method_flag = 21; else method_flag = 22; *************** *** 316,322 **** if (isopt) { ! if (user_jsub) method_flag = 11; else method_flag = 12; --- 351,357 ---- if (isopt) { ! if (user_jsub || user_jsub_cached) method_flag = 11; else method_flag = 12; *************** *** 377,383 **** if (! sanity_checked) { ! ColumnVector fval = user_fsub (x, t, theta); if (fval.length () != x.length ()) { --- 412,430 ---- if (! sanity_checked) { ! ColumnVector fval; ! if (user_fsub) ! fval = user_fsub (x, t, theta); ! else if (user_fsub_cached) ! fval = user_fsub_cached (x, t, theta, -1); ! else ! { ! (*current_liboctave_error_handler) ! ("odessa: no user supplied RHS subroutine!"); ! ! integration_error = true; ! return; ! } if (fval.length () != x.length ()) { *** ./liboctave/ODESSA-opts.in.orig 2004-08-06 14:35:27.000000000 +0200 --- ./liboctave/ODESSA-opts.in 2004-08-23 18:07:33.000000000 +0200 *************** *** 129,131 **** --- 129,154 ---- SET_EXPR = "val" END_OPTION + OPTION + NAME = "jacobian" + DOC_ITEM + Flag whether the user function returns the Jacobian as a second argument. + The Jacobian can also be entered as a seperate equation as a string matrix + to @code{odessa}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION + + OPTION + NAME = "inhomogeneity" + DOC_ITEM + Flag whether the user function returns the inhomogeneity matrix as a third + arguement. The function to calculate the inhomogenity matrix can also entered + as a seperate equation as the third row of a string matrix to @code{odessa}. + END_DOC_ITEM + TYPE = "int" + INIT_VALUE = "0" + SET_EXPR = "val" + END_OPTION *** ./liboctave/DASSL.cc.orig 2004-08-08 17:12:48.000000000 +0200 --- ./liboctave/DASSL.cc 2004-08-08 18:40:32.000000000 +0200 *************** *** 36,41 **** --- 36,42 ---- #include "lo-error.h" #include "lo-sstream.h" #include "quit.h" + #include "lo-ieee.h" typedef int (*dassl_fcn_ptr) (const double&, const double*, const double*, double*, int&, double*, int*); *************** *** 56,61 **** --- 57,64 ---- static DAEFunc::DAERHSFunc user_fun; static DAEFunc::DAEJacFunc user_jac; + static DAEFunc::DAERHSFuncCached user_fun_cached; + static DAEFunc::DAEJacFuncCached user_jac_cached; static int nn; *************** *** 77,83 **** tmp_state.elem (i) = state [i]; } ! tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires); if (ires >= 0) { --- 80,92 ---- tmp_state.elem (i) = state [i]; } ! if (user_fun_cached) ! // XXX FIXME XXX The underlying fortran doesn't give us the value of ! // cj!! Set it to an illegal value, so that the jacobian in the cache ! // will be ignored ! tmp_delta = user_fun_cached (tmp_state, tmp_deriv, time, ires, octave_Inf); ! else ! tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires); if (ires >= 0) { *************** *** 112,118 **** tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) --- 121,133 ---- tmp_state.elem (i) = state [i]; } ! Matrix tmp_pd; ! int ires; ! ! if (user_jac_cached) ! tmp_pd = user_jac_cached (tmp_state, tmp_deriv, time, ires, cj); ! else ! tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) *************** *** 171,178 **** // DAEFunc ! user_fun = DAEFunc::function (); ! user_jac = DAEFunc::jacobian_function (); if (user_fun) { --- 186,205 ---- // DAEFunc ! if (jacobian ()) ! { ! user_fun = 0; ! user_jac = 0; ! user_fun_cached = DAEFunc::function_cached (); ! user_jac_cached = DAEFunc::jacobian_function_cached (); ! } ! else ! { ! user_fun = DAEFunc::function (); ! user_jac = DAEFunc::jacobian_function (); ! user_fun_cached = 0; ! user_jac_cached = 0; ! } if (user_fun) { *************** *** 189,194 **** --- 216,237 ---- return retval; } } + else if (user_fun_cached) + { + int ires = 0; + double dummy = 0; + + ColumnVector res = (*user_fun_cached) (x, xdot, t, ires, dummy); + + if (res.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("dassl: inconsistent sizes for state and residual vectors"); + + integration_error = true; + return retval; + } + } else { (*current_liboctave_error_handler) *************** *** 198,204 **** return retval; } ! info(4) = user_jac ? 1 : 0; DAEFunc::reset = false; --- 241,247 ---- return retval; } ! info(4) = (user_jac || user_jac_cached) ? 1 : 0; DAEFunc::reset = false; *** ./src/DLD-FUNCTIONS/fsolve.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/DLD-FUNCTIONS/fsolve.cc 2004-08-06 18:38:12.000000000 +0200 *************** *** 49,54 **** --- 49,57 ---- // Global pointer for optional user defined jacobian function. static octave_function *fsolve_jac; + // Global pointer to a user function that returns the jacobian as a second arg + static octave_function *fsolve_fcn_jac; + // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; static bool warned_jac_imaginary = false; *************** *** 170,179 **** if (tmp.length () > 0 && tmp(0).is_defined ()) { ! if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) { warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); ! warned_fcn_imaginary = true; } retval = tmp(0).matrix_value (); --- 173,182 ---- if (tmp.length () > 0 && tmp(0).is_defined ()) { ! if (! warned_jac_imaginary && tmp(0).is_complex_type ()) { warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); ! warned_jac_imaginary = true; } retval = tmp(0).matrix_value (); *************** *** 188,193 **** --- 191,347 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_x; + + ColumnVector + fsolve_user_function_cached (const ColumnVector& x) + { + ColumnVector retval; + + if (cached_x == x) + retval = cached_fcn; + else + { + int n = x.length (); + + octave_value_list args; + args.resize (1); + + if (n > 1) + { + Matrix m (n, 1); + for (int i = 0; i < n; i++) + m (i, 0) = x (i); + octave_value vars (m); + args(0) = vars; + } + else + { + double d = x (0); + octave_value vars (d); + args(0) = vars; + } + + if (fsolve_fcn_jac) + { + octave_value_list tmp = fsolve_fcn_jac->do_multi_index_op (2, args); + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("fsolve: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + + if (error_state || cached_fcn.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + + if (error_state || cached_jac.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + + cached_x = x; + retval = cached_fcn; + } + else + gripe_user_supplied_eval ("fsolve"); + } + } + + return retval; + } + + Matrix + fsolve_user_jacobian_cached (const ColumnVector& x) + { + Matrix retval; + + if (cached_x == x) + retval = cached_jac; + else + { + int n = x.length (); + + octave_value_list args; + args.resize (1); + + if (n > 1) + { + Matrix m (n, 1); + for (int i = 0; i < n; i++) + m (i, 0) = x (i); + octave_value vars (m); + args(0) = vars; + } + else + { + double d = x (0); + octave_value vars (d); + args(0) = vars; + } + + if (fsolve_fcn_jac) + { + octave_value_list tmp = fsolve_fcn_jac->do_multi_index_op (2, args); + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("fsolve: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + + if (error_state || cached_fcn.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + + if (error_state || cached_jac.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + + cached_x = x; + retval = cached_jac; + } + else + gripe_user_supplied_eval ("fsolve"); + } + } + + return retval; + } + #define FSOLVE_ABORT() \ do \ { \ *************** *** 215,223 **** DEFUN_DLD (fsolve, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} address@hidden, @var{info}, @var{msg}] =} 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\ \n\ If @var{fcn} is a two-element string array, the first element names\n\ the function @math{f} described above, and the second element names\n\ --- 369,382 ---- DEFUN_DLD (fsolve, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} address@hidden, @var{info}, @var{msg}] =} fsolve (@var{fcn}, @var{x0})\n\ ! Given @var{fcn}, the name of a function, a function handle or an inline\n\ ! function of the form @code{f (@var{x})} and an initial starting point\n\ ! @var{x0}, @code{fsolve} solves the set of equations such that\n\ ! @code{f(@var{x}) == 0}.\n\ ! \n\ ! The function @var{fcn} can also return a second argument, in which can the\n\ ! second argument is the jacobian @code{j (@var{x})}. In this case the\n\ ! @code{jacobian} option of @code{fsolve_options} must be set.\n\ \n\ If @var{fcn} is a two-element string array, the first element names\n\ the function @math{f} described above, and the second element names\n\ *************** *** 256,301 **** if (nargin == 2 && nargout < 4) { fsolve_fcn = 0; fsolve_jac = 0; octave_value f_arg = args(0); ! switch (f_arg.rows ()) { ! case 1: ! fsolve_fcn = extract_function ! (f_arg, "fsolve", "__fsolve_fcn__", ! "function y = __fsolve_fcn__ (x) y = ", ! "; endfunction"); ! break; ! ! case 2: ! { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! fsolve_fcn = extract_function ! (tmp(0), "fsolve", "__fsolve_fcn__", ! "function y = __fsolve_fcn__ (x) y = ", ! "; endfunction"); ! if (fsolve_fcn) { ! fsolve_jac = extract_function ! (tmp(1), "fsolve", "__fsolve_jac__", ! "function jac = __fsolve_jac__ (x) jac = ", ! "; endfunction"); ! ! if (! fsolve_jac) ! fsolve_fcn = 0; } } ! } } ! if (error_state || ! fsolve_fcn) FSOLVE_ABORT (); ColumnVector x (args(1).vector_value ()); --- 415,495 ---- if (nargin == 2 && nargout < 4) { + std::string fcn_name, fname, jac_name, jname; fsolve_fcn = 0; fsolve_jac = 0; + fsolve_fcn_jac = 0; octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) ! { ! if (fsolve_opts.jacobian ()) ! fsolve_fcn_jac = f_arg.function_value (); ! else ! fsolve_fcn = f_arg.function_value (); ! } ! else { ! switch (f_arg.rows ()) ! { ! case 1: ! if (fsolve_opts.jacobian ()) ! { ! fcn_name = unique_symbol_name ("__fsolve_fcn_jac__"); ! fname = "function [y, jac] = "; ! fname.append (fcn_name); ! fname.append (" (x) [y, jac] = "); ! fsolve_fcn_jac = extract_function ! (f_arg, "fsolve", fcn_name, fname, "; endfunction"); ! } ! else ! { ! fcn_name = unique_symbol_name ("__fsolve_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x) y = "); ! fsolve_fcn = extract_function ! (f_arg, "fsolve", fcn_name, fname, "; endfunction"); ! } ! break; ! case 2: { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! fcn_name = unique_symbol_name ("__fsolve_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x) y = "); ! fsolve_fcn = extract_function ! (tmp(0), "fsolve", fcn_name, fname, "; endfunction"); ! ! if (fsolve_fcn) ! { ! jac_name = unique_symbol_name ("__fsolve_jac__"); ! jname = "function jac = "; ! jname.append(jac_name); ! jname.append (" (x) jac = "); ! fsolve_jac = extract_function ! (tmp(1), "fsolve", jac_name, jname, "; endfunction"); ! ! if (! fsolve_jac) ! fsolve_fcn = 0; ! } } } ! break; ! ! default: ! FSOLVE_ABORT1 ! ("first arg should be a string or 2-element string array"); ! } } ! if (error_state || ! (fsolve_fcn || fsolve_fcn_jac)) FSOLVE_ABORT (); ColumnVector x (args(1).vector_value ()); *************** *** 309,317 **** if (nargout > 3) warning ("fsolve: can't compute path output yet"); ! NLFunc nleqn_fcn (fsolve_user_function); ! if (fsolve_jac) ! nleqn_fcn.set_jacobian_function (fsolve_user_jacobian); NLEqn nleqn (x, nleqn_fcn); nleqn.set_options (fsolve_opts); --- 503,521 ---- if (nargout > 3) warning ("fsolve: can't compute path output yet"); ! NLFunc nleqn_fcn; ! ! if (fsolve_fcn) ! { ! nleqn_fcn.set_function (fsolve_user_function); ! if (fsolve_jac) ! nleqn_fcn.set_jacobian_function (fsolve_user_jacobian); ! } ! else ! { ! nleqn_fcn.set_function (fsolve_user_function_cached); ! nleqn_fcn.set_jacobian_function (fsolve_user_jacobian_cached); ! } NLEqn nleqn (x, nleqn_fcn); nleqn.set_options (fsolve_opts); *************** *** 319,324 **** --- 523,533 ---- int info; ColumnVector soln = nleqn.solve (info); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (! error_state) { std::string msg = nleqn.error_message (); *************** *** 331,336 **** --- 540,548 ---- if (! nleqn.solution_ok () && nargout < 2) error ("fsolve: %s", msg.c_str ()); } + + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); } else print_usage ("fsolve"); *** ./src/DLD-FUNCTIONS/quad.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/DLD-FUNCTIONS/quad.cc 2004-08-06 12:11:13.000000000 +0200 *************** *** 130,137 **** "-*- texinfo -*-\n\ @deftypefn {Loadable Function} address@hidden, @var{ier}, @var{nfun}, @var{err}] =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ Integrate a nonlinear function of one variable using Quadpack.\n\ ! The first argument is the name of the function to call to compute the\n\ ! value of the integrand. It must have the form\n\ \n\ @example\n\ y = f (x)\n\ --- 130,138 ---- "-*- texinfo -*-\n\ @deftypefn {Loadable Function} address@hidden, @var{ier}, @var{nfun}, @var{err}] =} quad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ Integrate a nonlinear function of one variable using Quadpack.\n\ ! The first argument is the name of the function, the function handle or\n\ ! the inline function to call to compute the value of the integrand. It\n\ ! must have the form\n\ \n\ @example\n\ y = f (x)\n\ *************** *** 179,187 **** if (nargin > 2 && nargin < 6 && nargout < 5) { ! quad_fcn = extract_function (args(0), "quad", "__quad_fcn__", ! "function y = __quad_fcn__ (x) y = ", ! "; endfunction"); if (! quad_fcn) QUAD_ABORT (); --- 180,192 ---- if (nargin > 2 && nargin < 6 && nargout < 5) { ! if (args(0).is_function_handle () || args(0).is_inline_function ()) ! quad_fcn = args(0).function_value (); ! else ! quad_fcn = extract_function (args(0), "quad", "__quad_fcn__", ! "function y = __quad_fcn__ (x) y = ", ! "; endfunction"); ! if (! quad_fcn) QUAD_ABORT (); *** ./src/DLD-FUNCTIONS/daspk.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/DLD-FUNCTIONS/daspk.cc 2004-08-08 18:36:47.000000000 +0200 *************** *** 49,54 **** --- 49,57 ---- // Global pointer for optional user defined jacobian function. static octave_function *daspk_jac; + // Global pointer to a user function that returns the jacobian as a second arg + static octave_function *daspk_fcn_jac; + // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; static bool warned_jac_imaginary = false; *************** *** 150,155 **** --- 153,324 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_x; + static ColumnVector cached_xdot; + static double cached_t; + static double cached_cj; + static int cached_ires; + + ColumnVector + daspk_user_function_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int& ires, double cj) + { + ColumnVector retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t) + { + retval = cached_fcn; + ires = cached_ires; + } + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (daspk_fcn) + { + octave_value_list tmp = daspk_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("daspk"); + return retval; + } + + int tlen = tmp.length (); + if (tlen > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + retval = cached_fcn; + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + if (tlen > 2) + { + cached_ires = tmp(2).int_value (); + ires = cached_ires; + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("daspk"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + cached_ires = ires; + } + else + gripe_user_supplied_eval ("daspk"); + } + } + + return retval; + } + + Matrix + daspk_user_jacobian_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int&, double cj) + { + Matrix retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t && cached_cj == cj) + retval = cached_jac; + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (daspk_fcn) + { + octave_value_list tmp = daspk_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("daspk"); + return retval; + } + + int tlen = tmp.length (); + if (tlen > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + retval = cached_jac; + } + + if (tlen > 2) + { + cached_ires = tmp(2).int_value (); + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("daspk"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + } + else + gripe_user_supplied_eval ("daspk"); + } + } + + return retval; + } + #define DASPK_ABORT() \ do \ { \ *************** *** 204,212 **** row of the output @var{x} is @var{x_0} and the first row\n\ of the output @var{xdot} is @var{xdot_0}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function to\n\ ! call to compute the vector of residuals for the set of equations.\n\ ! It must have the form\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ --- 373,381 ---- row of the output @var{x} is @var{x_0} and the first row\n\ of the output @var{xdot} is @var{xdot_0}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function, a\n\ ! function handle, or an inline function, to call to compute the vector\n\ ! of residuals for the set of equations. It must have the form\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ *************** *** 216,224 **** in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ scalar.\n\ \n\ ! If @var{fcn} is a two-element string array, the first element names\n\ ! the function @math{f} described above, and the second element names\n\ ! a function to compute the modified Jacobian\n\ @tex\n\ $$\n\ J = {\\partial f \\over \\partial x}\n\ --- 385,392 ---- in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ scalar.\n\ \n\ ! If the function @var{fcn} can return a second argument representing a\n\ ! function to compute the modified Jacobian\n\ @tex\n\ $$\n\ J = {\\partial f \\over \\partial x}\n\ *************** *** 242,247 **** --- 410,422 ---- \n\ @end example\n\ \n\ + In this case the presence of the second return argument must be flagged\n\ + with the @code{jacobian} option of @code{daspk_options}.\n\ + \n\ + If @var{fcn} is a two-element string array, the first element names\n\ + the function @math{f} described above, and the second element names\n\ + the modified Jacobian\n\ + \n\ The second and third arguments to @code{daspk} specify the initial\n\ condition of the states and their derivatives, and the fourth argument\n\ specifies a vector of output times at which the solution is desired,\n\ *************** *** 285,330 **** if (nargin > 3 && nargin < 6) { daspk_fcn = 0; daspk_jac = 0; octave_value f_arg = args(0); ! switch (f_arg.rows ()) { ! case 1: ! daspk_fcn = extract_function ! (args(0), "daspk", "__daspk_fcn__", ! "function res = __daspk_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! break; ! ! case 2: ! { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! daspk_fcn = extract_function ! (tmp(0), "daspk", "__daspk_fcn__", ! "function res = __daspk_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! if (daspk_fcn) { ! daspk_jac = extract_function ! (tmp(1), "daspk", "__daspk_jac__", ! "function jac = __daspk_jac__ (x, xdot, t, cj) jac = ", ! "; endfunction"); ! ! if (! daspk_jac) ! daspk_fcn = 0; } } ! } } ! if (error_state || ! daspk_fcn) DASPK_ABORT (); ColumnVector state = ColumnVector (args(1).vector_value ()); --- 460,535 ---- if (nargin > 3 && nargin < 6) { + std::string fcn_name, fname, jac_name, jname; daspk_fcn = 0; daspk_jac = 0; + daspk_fcn_jac = 0; octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) { ! if (daspk_opts.jacobian ()) ! daspk_fcn_jac = f_arg.function_value (); ! else ! daspk_fcn = f_arg.function_value (); ! } ! else ! { ! switch (f_arg.rows ()) ! { ! case 1: ! if (daspk_opts.jacobian ()) ! { ! fcn_name = unique_symbol_name ("__daspk_fcn_jac__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t, cj) y = "); ! daspk_fcn_jac = extract_function ! (args(0), "daspk", fcn_name, fname, "; endfunction"); ! } ! else ! { ! fcn_name = unique_symbol_name ("__daspk_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! daspk_fcn = extract_function ! (args(0), "daspk", fcn_name, fname, "; endfunction"); ! } ! break; ! case 2: { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! fcn_name = unique_symbol_name ("__daspk_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! daspk_fcn = extract_function ! (tmp(0), "daspk", fcn_name, fname, "; endfunction"); ! ! if (daspk_fcn) ! { ! jac_name = unique_symbol_name ("__daspk_jac__"); ! jname = "function jac = "; ! jname.append(jac_name); ! jname.append (" (x, xdot, t, cj) jac = "); ! daspk_jac = extract_function ! (tmp(1), "daspk", jac_name, jname, "; endfunction"); ! ! if (! daspk_jac) ! daspk_fcn = 0; ! } } } ! } } ! if (error_state || ! (daspk_fcn || daspk_fcn_jac)) DASPK_ABORT (); ColumnVector state = ColumnVector (args(1).vector_value ()); *************** *** 359,367 **** double tzero = out_times (0); ! DAEFunc func (daspk_user_function); ! if (daspk_jac) ! func.set_jacobian_function (daspk_user_jacobian); DASPK dae (state, deriv, tzero, func); dae.set_options (daspk_opts); --- 564,581 ---- double tzero = out_times (0); ! DAEFunc func; ! if (daspk_fcn) ! { ! func.set_function (daspk_user_function); ! if (daspk_jac) ! func.set_jacobian_function (daspk_user_jacobian); ! } ! else ! { ! func.set_function_cached (daspk_user_function_cached); ! func.set_jacobian_function_cached (daspk_user_jacobian_cached); ! } DASPK dae (state, deriv, tzero, func); dae.set_options (daspk_opts); *************** *** 374,379 **** --- 588,598 ---- else output = dae.integrate (out_times, deriv_output); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (! error_state) { std::string msg = dae.error_message (); *************** *** 395,400 **** --- 614,623 ---- error ("daspk: %s", msg.c_str ()); } } + + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); + } else print_usage ("daspk"); *** ./src/DLD-FUNCTIONS/dasrt.cc.orig 2004-08-06 17:43:01.000000000 +0200 --- ./src/DLD-FUNCTIONS/dasrt.cc 2004-08-24 10:24:38.000000000 +0200 *************** *** 47,52 **** --- 47,54 ---- static octave_function *dasrt_f; static octave_function *dasrt_j; static octave_function *dasrt_cf; + static octave_function *dasrt_fj; + // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; *************** *** 186,191 **** --- 188,341 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_x; + static ColumnVector cached_xdot; + static double cached_t; + static double cached_cj; + + static ColumnVector + dasrt_user_f_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int&, double cj) + { + ColumnVector retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t) + retval = cached_fcn; + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (dasrt_fj) + { + octave_value_list tmp = dasrt_fj->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("dasrt"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + retval = cached_fcn; + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("dasrt"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + } + else + gripe_user_supplied_eval ("dasrt"); + } + } + + return retval; + } + + static Matrix + dasrt_user_j_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int&, double cj) + { + Matrix retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t && cached_cj == cj) + retval = cached_jac; + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (dasrt_fj) + { + octave_value_list tmp = dasrt_fj->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("dasrt"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + retval = cached_jac; + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("dasrt"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + } + else + gripe_user_supplied_eval ("dasrt"); + } + } + + return retval; + } + #define DASRT_ABORT \ do \ { \ *************** *** 248,256 **** @var{t_out} will be the point at which the stopping condition was met,\n\ and may not correspond to any element of the vector @var{t}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function to\n\ ! call to compute the vector of residuals for the set of equations.\n\ ! It must have the form\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ --- 398,406 ---- @var{t_out} will be the point at which the stopping condition was met,\n\ and may not correspond to any element of the vector @var{t}.\n\ \n\ ! The first argument, @var{fcn}, is either a string that names the function,\n\ ! a function handle, or an inline function, to call to compute the vector of\n\ ! residuals for the set of equations. It must have the form\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ *************** *** 260,268 **** in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ scalar.\n\ \n\ ! If @var{fcn} is a two-element string array, the first element names\n\ ! the function @math{f} described above, and the second element names\n\ ! a function to compute the modified Jacobian\n\ \n\ @tex\n\ $$\n\ --- 410,418 ---- in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ scalar.\n\ \n\ ! The function @var{fcn} can also return a second argument, in which can the\n\ ! second argument is the jacobian @code{j (@var{x}, @var{t})}. In this case\n\ ! the @code{jacobian} option of @code{dasrt_options} must be set.\n\ \n\ @tex\n\ $$\n\ *************** *** 288,293 **** --- 438,447 ---- \n\ @end example\n\ \n\ + If @var{fcn} is a two-element string array, the first element names\n\ + the function @math{f} described above, and the second element names\n\ + a function to compute the modified Jacobian\n\ + \n\ The optional second argument names a function that defines the\n\ constraint functions whose roots are desired during the integration.\n\ This function must have the form\n\ *************** *** 366,374 **** --- 520,530 ---- return retval; } + std::string fcn_name, fname, jac_name, jname; dasrt_f = 0; dasrt_j = 0; dasrt_cf = 0; + dasrt_fj = 0; // Check all the arguments. Are they the right animals? *************** *** 376,428 **** octave_value f_arg = args(0); ! switch (f_arg.rows ()) { ! case 1: ! dasrt_f = extract_function ! (args(0), "dasrt", "__dasrt_fcn__", ! "function res = __dasrt_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! break; ! ! case 2: ! { ! string_vector tmp = args(0).all_strings (); ! ! if (! error_state) { ! dasrt_f = extract_function ! (tmp(0), "dasrt", "__dasrt_fcn__", ! "function res = __dasrt_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! if (dasrt_f) { ! dasrt_j = extract_function ! (tmp(1), "dasrt", "__dasrt_jac__", ! "function jac = __dasrt_jac__ (x, xdot, t, cj) jac = ", ! "; endfunction"); ! ! if (! dasrt_j) ! dasrt_f = 0; } } ! } ! break; ! default: ! DASRT_ABORT1 ! ("first arg should be a string or 2-element string array"); } ! if (error_state || (! dasrt_f)) DASRT_ABORT; ! DAERTFunc func (dasrt_user_f); argp++; ! if (args(1).is_string ()) { dasrt_cf = is_valid_function (args(1), "dasrt", true); --- 532,627 ---- octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) { ! if (dasrt_opts.jacobian ()) ! dasrt_fj = f_arg.function_value (); ! else ! dasrt_f = f_arg.function_value (); ! } ! else ! { ! switch (f_arg.rows ()) ! { ! case 1: ! if (dasrt_opts.jacobian ()) ! { ! fcn_name = unique_symbol_name ("__dasrt_fcn_jac__"); ! fname = "function [y, jac] = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t, cj) [y, jac] = "); ! dasrt_fj = extract_function ! (f_arg, "dasrt", fcn_name, fname, "; endfunction"); ! } ! else ! { ! fcn_name = unique_symbol_name ("__dasrt_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! dasrt_f = extract_function ! (f_arg, "dasrt", fcn_name, fname, "; endfunction"); ! } ! break; ! ! case 2: { ! string_vector tmp = args(0).all_strings (); ! if (! error_state) { ! fcn_name = unique_symbol_name ("__dasrt_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! dasrt_f = extract_function ! (tmp(0), "dasrt", fcn_name, fname, "; endfunction"); ! ! if (dasrt_f) ! { ! jac_name = unique_symbol_name ("__dasrt_jac__"); ! jname = "function jac = "; ! jname.append(jac_name); ! jname.append (" (x, xdot, t, cj) jac = "); ! dasrt_j = extract_function ! (tmp(1), "dasrt", jac_name, jname, "; endfunction"); ! ! if (! dasrt_j) ! dasrt_f = 0; ! } } } ! break; ! default: ! DASRT_ABORT1 ! ("first arg should be a string or 2-element string array"); ! } } ! if (error_state || (! (dasrt_f || dasrt_fj))) DASRT_ABORT; ! DAERTFunc func; ! if (dasrt_f) ! func.set_function (dasrt_user_f); ! else ! func.set_function_cached (dasrt_user_f_cached); argp++; ! if (args(1).is_function_handle() || args(1).is_inline_function()) ! { ! dasrt_cf = args(1).function_value(); ! ! if (! dasrt_cf) ! DASRT_ABORT1 ("expecting function name as argument 2"); ! ! argp++; ! ! func.set_constraint_function (dasrt_user_cf); ! } ! else if (args(1).is_string ()) { dasrt_cf = is_valid_function (args(1), "dasrt", true); *************** *** 470,475 **** --- 669,676 ---- if (dasrt_j) func.set_jacobian_function (dasrt_user_j); + else if (dasrt_fj) + func.set_jacobian_function_cached (dasrt_user_j_cached); DASRT_result output; *************** *** 482,487 **** --- 683,693 ---- else output = dae.integrate (out_times); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (! error_state) { std::string msg = dae.error_message (); *************** *** 506,511 **** --- 712,720 ---- } } + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); + unwind_protect::run_frame ("Fdasrt"); return retval; *** ./src/DLD-FUNCTIONS/lsode.cc.orig 2004-08-06 18:43:55.000000000 +0200 --- ./src/DLD-FUNCTIONS/lsode.cc 2004-08-06 19:09:13.000000000 +0200 *************** *** 51,56 **** --- 51,59 ---- // Global pointer for optional user defined jacobian function used by lsode. static octave_function *lsode_jac; + // Global pointer to a user function that returns the jacobian as a second arg + static octave_function *lsode_fcn_jac; + // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; static bool warned_jac_imaginary = false; *************** *** 136,141 **** --- 139,277 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_x; + static double cached_t; + + + ColumnVector + lsode_user_function_cached (const ColumnVector& x, double t) + { + ColumnVector retval; + + if (cached_x == x && cached_t == t) + retval = cached_fcn; + else + { + octave_value_list args; + args(1) = t; + args(0) = x; + + if (lsode_fcn_jac) + { + octave_value_list tmp = lsode_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("lsode"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + retval = cached_fcn; + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("lsode"); + + cached_x = x; + cached_t = t; + } + else + gripe_user_supplied_eval ("lsode"); + } + } + + return retval; + } + + Matrix + lsode_user_jacobian_cached (const ColumnVector& x, double t) + { + Matrix retval; + + if (cached_x == x && cached_t == t) + retval = cached_jac; + else + { + octave_value_list args; + args(1) = t; + args(0) = x; + + if (lsode_fcn_jac) + { + octave_value_list tmp = lsode_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("lsode"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + retval = cached_jac; + } + + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("lsode"); + + cached_x = x; + cached_t = t; + } + else + gripe_user_supplied_eval ("lsode"); + } + } + + return retval; + } + #define LSODE_ABORT() \ do \ { \ *************** *** 190,198 **** state of the system @var{x_0}, so that the first row of the output\n\ is @var{x_0}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function to\n\ ! call to compute the vector of right hand sides for the set of equations.\n\ ! The function must have the form\n\ \n\ @example\n\ @var{xdot} = f (@var{x}, @var{t})\n\ --- 326,334 ---- state of the system @var{x_0}, so that the first row of the output\n\ is @var{x_0}.\n\ \n\ ! The first argument, @var{fcn}, is either a string that names the function,\n\ ! a function handle, or an inline function, to call to compute the vector of\n\ ! right hand sides for the set of equations. The function must have the form\n\ \n\ @example\n\ @var{xdot} = f (@var{x}, @var{t})\n\ *************** *** 201,209 **** @noindent\n\ in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\ \n\ ! If @var{fcn} is a two-element string array, the first element names the\n\ ! function @math{f} described above, and the second element names a function\n\ ! to compute the Jacobian of @math{f}. The Jacobian function must have the\n\ form\n\ \n\ @example\n\ --- 337,347 ---- @noindent\n\ in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\ \n\ ! The function @var{fcn} can also return a second argument, in which can the\n\ ! second argument is the jacobian @code{j (@var{x}, @var{t})}. In this case\n\ ! the @code{jacobian} option of @code{lsode_options} must be set.\n\ ! \n\ ! The Jacobian function must have the\n\ form\n\ \n\ @example\n\ *************** *** 249,254 **** --- 387,396 ---- \n\ @end ifinfo\n\ \n\ + If @var{fcn} is a two-element string array, the first element names the\n\ + function @math{f} described above, and the second element names a function\n\ + to compute the Jacobian of @math{f}.\n\ + \n\ The second and third arguments specify the intial state of the system,\n\ @math{x_0}, and the initial value of the independent variable @math{t_0}.\n\ \n\ *************** *** 285,335 **** if (nargin > 2 && nargin < 5 && nargout < 4) { lsode_fcn = 0; lsode_jac = 0; octave_value f_arg = args(0); ! switch (f_arg.rows ()) { - case 1: - lsode_fcn = extract_function - (f_arg, "lsode", "__lsode_fcn__", - "function xdot = __lsode_fcn__ (x, t) xdot = ", - "; endfunction"); - break; - - case 2: - { - string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! lsode_fcn = extract_function ! (tmp(0), "lsode", "__lsode_fcn__", ! "function xdot = __lsode_fcn__ (x, t) xdot = ", ! "; endfunction"); ! if (lsode_fcn) { ! lsode_jac = extract_function ! (tmp(1), "lsode", "__lsode_jac__", ! "function jac = __lsode_jac__ (x, t) jac = ", ! "; endfunction"); ! ! if (! lsode_jac) ! lsode_fcn = 0; } } ! } ! break; ! default: ! LSODE_ABORT1 ! ("first arg should be a string or 2-element string array"); } ! if (error_state || ! lsode_fcn) LSODE_ABORT (); ColumnVector state (args(1).vector_value ()); --- 427,508 ---- if (nargin > 2 && nargin < 5 && nargout < 4) { + std::string fcn_name, fname, jac_name, jname; lsode_fcn = 0; lsode_jac = 0; + lsode_fcn_jac = 0; octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) ! { ! if (lsode_opts.jacobian ()) ! lsode_fcn_jac = f_arg.function_value (); ! else ! lsode_fcn = f_arg.function_value (); ! } ! else { ! switch (f_arg.rows ()) ! { ! case 1: ! if (lsode_opts.jacobian ()) ! { ! fcn_name = unique_symbol_name ("__lsode_fcn_jac__"); ! fname = "function [y, jac] = "; ! fname.append (fcn_name); ! fname.append (" (x, t) [y, jac] = "); ! lsode_fcn_jac = extract_function ! (f_arg, "lsode", fcn_name, fname, "; endfunction"); ! } ! else ! { ! fcn_name = unique_symbol_name ("__lsode_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, t) y = "); ! lsode_fcn = extract_function ! (f_arg, "lsode", fcn_name, fname, "; endfunction"); ! } ! break; ! ! case 2: { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! fcn_name = unique_symbol_name ("__lsode_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, t) y = "); ! lsode_fcn = extract_function ! (tmp(0), "lsode", fcn_name, fname, "; endfunction"); ! ! if (lsode_fcn) ! { ! jac_name = unique_symbol_name ("__lsode_jac__"); ! jname = "function jac = "; ! jname.append(jac_name); ! jname.append (" (x, t) jac = "); ! lsode_jac = extract_function ! (tmp(1), "lsode", jac_name, jname, "; endfunction"); ! ! if (! lsode_jac) ! lsode_fcn = 0; ! } } } ! break; ! default: ! LSODE_ABORT1 ! ("first arg should be a string or 2-element string array"); ! } } ! if (error_state || ! (lsode_fcn || lsode_fcn_jac)) LSODE_ABORT (); ColumnVector state (args(1).vector_value ()); *************** *** 357,365 **** double tzero = out_times (0); ! ODEFunc func (lsode_user_function); ! if (lsode_jac) ! func.set_jacobian_function (lsode_user_jacobian); LSODE ode (state, tzero, func); --- 530,547 ---- double tzero = out_times (0); ! ODEFunc func; ! if (lsode_fcn) ! { ! func.set_function (lsode_user_function); ! if (lsode_jac) ! func.set_jacobian_function (lsode_user_jacobian); ! } ! else ! { ! func.set_function (lsode_user_function_cached); ! func.set_jacobian_function (lsode_user_jacobian_cached); ! } LSODE ode (state, tzero, func); *************** *** 371,376 **** --- 553,563 ---- else output = ode.integrate (out_times); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (! error_state) { std::string msg = ode.error_message (); *************** *** 388,393 **** --- 575,583 ---- error ("lsode: %s", msg.c_str ()); } } + + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); } else print_usage ("lsode"); *** ./src/DLD-FUNCTIONS/dassl.cc.orig 2004-08-08 18:38:05.000000000 +0200 --- ./src/DLD-FUNCTIONS/dassl.cc 2004-08-08 19:04:57.000000000 +0200 *************** *** 49,54 **** --- 49,57 ---- // Global pointer for optional user defined jacobian function. static octave_function *dassl_jac; + // Global pointer for user function that returns the jacobian as a second arg + static octave_function *dassl_fcn_jac; + // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; static bool warned_jac_imaginary = false; *************** *** 150,155 **** --- 153,319 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_x; + static ColumnVector cached_xdot; + static double cached_t; + static int cached_ires; + static double cached_cj; + + ColumnVector + dassl_user_function_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int& ires, double cj) + { + ColumnVector retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t) + { + retval = cached_fcn; + ires = cached_ires; + } + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (dassl_fcn_jac) + { + octave_value_list tmp = dassl_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("dassl"); + return retval; + } + + int tlen = tmp.length (); + if (tlen > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + retval = cached_fcn; + } + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function"); + warned_fcn_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + if (tlen > 2) + { + cached_ires = tmp(2).int_value (); + ires = cached_ires; + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("dassl"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + } + else + gripe_user_supplied_eval ("dassl"); + } + } + + return retval; + } + + Matrix + dassl_user_jacobian_cached (const ColumnVector& x, const ColumnVector& xdot, + double t, int&, double cj) + { + Matrix retval; + + if (cached_x == x && cached_xdot == xdot && cached_t == t && cached_cj == cj) + retval = cached_jac; + else + { + assert (x.capacity () == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + args(1) = xdot; + args(0) = x; + + if (dassl_fcn_jac) + { + octave_value_list tmp = dassl_fcn_jac->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("dassl"); + return retval; + } + + int tlen = tmp.length (); + if (tlen > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function"); + warned_fcn_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + retval = cached_jac; + } + + if (tlen > 2) + cached_ires = tmp(2).int_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("dassl"); + + cached_x = x; + cached_xdot = xdot; + cached_t = t; + cached_cj = cj; + } + else + gripe_user_supplied_eval ("dassl"); + } + } + + return retval; + } + #define DASSL_ABORT() \ do \ { \ *************** *** 207,215 **** row of the output @var{x} is @var{x_0} and the first row\n\ of the output @var{xdot} is @var{xdot_0}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function to\n\ ! call to compute the vector of residuals for the set of equations.\n\ ! It must have the form\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ --- 371,382 ---- row of the output @var{x} is @var{x_0} and the first row\n\ of the output @var{xdot} is @var{xdot_0}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function,\n\ ! a function handle, or an inline function, to call to compute the vector\n\ ! of residuals for the set of equations. It must have the form\n\ ! \n\ ! If the function @var{fcn} can return a second argument representing a\n\ ! function to compute the modified Jacobian\n\ \n\ @example\n\ @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ *************** *** 219,228 **** in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ scalar.\n\ \n\ - If @var{fcn} is a two-element string array, the first element names\n\ - the function @math{f} described above, and the second element names\n\ - a function to compute the modified Jacobian\n\ - \n\ @iftex\n\ @tex\n\ $$\n\ --- 386,391 ---- *************** *** 247,252 **** --- 410,422 ---- \n\ @end example\n\ \n\ + In this case the presence of the second return argument must be flagged\n\ + with the @code{jacobian} option of @code{dassl_options}.\n\ + \n\ + If @var{fcn} is a two-element string array, the first element names\n\ + the function @math{f} described above, and the second element names\n\ + a function to compute the modified Jacobian\n\ + \n\ The second and third arguments to @code{dassl} specify the initial\n\ condition of the states and their derivatives, and the fourth argument\n\ specifies a vector of output times at which the solution is desired,\n\ *************** *** 290,335 **** if (nargin > 3 && nargin < 6 && nargout < 5) { dassl_fcn = 0; dassl_jac = 0; octave_value f_arg = args(0); ! switch (f_arg.rows ()) { ! case 1: ! dassl_fcn = extract_function ! (f_arg, "dassl", "__dassl_fcn__", ! "function res = __dassl_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! break; ! ! case 2: ! { ! string_vector tmp = f_arg.all_strings (); ! if (! error_state) { ! dassl_fcn = extract_function ! (tmp(0), "dassl", "__dassl_fcn__", ! "function res = __dassl_fcn__ (x, xdot, t) res = ", ! "; endfunction"); ! ! if (dassl_fcn) { ! dassl_jac = extract_function ! (tmp(1), "dassl", "__dassl_jac__", ! "function jac = __dassl_jac__ (x, xdot, t, cj) jac = ", ! "; endfunction"); ! ! if (! dassl_jac) ! dassl_fcn = 0; } } ! } } ! if (error_state || ! dassl_fcn) DASSL_ABORT (); ColumnVector state = ColumnVector (args(1).vector_value ()); --- 460,535 ---- if (nargin > 3 && nargin < 6 && nargout < 5) { + std::string fcn_name, fname, jac_name, jname; dassl_fcn = 0; dassl_jac = 0; + dassl_fcn_jac = 0; octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) { ! if (dassl_opts.jacobian ()) ! dassl_fcn_jac = f_arg.function_value (); ! else ! dassl_fcn = f_arg.function_value (); ! } ! else ! { ! switch (f_arg.rows ()) ! { ! case 1: ! if (dassl_opts.jacobian ()) ! { ! fcn_name = unique_symbol_name ("__dassl_fcn_jac__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t, cj) y = "); ! dassl_fcn_jac = extract_function ! (args(0), "dassl", fcn_name, fname, "; endfunction"); ! } ! else ! { ! fcn_name = unique_symbol_name ("__dassl_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! dassl_fcn = extract_function ! (args(0), "dassl", fcn_name, fname, "; endfunction"); ! } ! break; ! case 2: { ! string_vector tmp = f_arg.all_strings (); ! ! if (! error_state) { ! fcn_name = unique_symbol_name ("__dassl_fcn__"); ! fname = "function y = "; ! fname.append (fcn_name); ! fname.append (" (x, xdot, t) y = "); ! dassl_fcn = extract_function ! (tmp(0), "dassl", fcn_name, fname, "; endfunction"); ! ! if (dassl_fcn) ! { ! jac_name = unique_symbol_name ("__dassl_jac__"); ! jname = "function jac = "; ! jname.append(jac_name); ! jname.append (" (x, xdot, t, cj) jac = "); ! dassl_jac = extract_function ! (tmp(1), "dassl", jac_name, jname, "; endfunction"); ! ! if (! dassl_jac) ! dassl_fcn = 0; ! } } } ! } } ! if (error_state || ! (dassl_fcn || dassl_fcn_jac)) DASSL_ABORT (); ColumnVector state = ColumnVector (args(1).vector_value ()); *************** *** 364,372 **** double tzero = out_times (0); ! DAEFunc func (dassl_user_function); ! if (dassl_jac) ! func.set_jacobian_function (dassl_user_jacobian); DASSL dae (state, deriv, tzero, func); --- 564,581 ---- double tzero = out_times (0); ! DAEFunc func; ! if (dassl_fcn) ! { ! func.set_function (dassl_user_function); ! if (dassl_jac) ! func.set_jacobian_function (dassl_user_jacobian); ! } ! else ! { ! func.set_function_cached (dassl_user_function_cached); ! func.set_jacobian_function_cached (dassl_user_jacobian_cached); ! } DASSL dae (state, deriv, tzero, func); *************** *** 380,385 **** --- 589,599 ---- else output = dae.integrate (out_times, deriv_output); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (! error_state) { std::string msg = dae.error_message (); *************** *** 401,406 **** --- 615,623 ---- error ("dassl: %s", msg.c_str ()); } } + + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); } else print_usage ("dassl"); *** ./src/DLD-FUNCTIONS/odessa.cc.orig 2004-08-08 19:15:43.000000000 +0200 --- ./src/DLD-FUNCTIONS/odessa.cc 2004-08-24 10:53:55.000000000 +0200 *************** *** 50,55 **** --- 50,56 ---- static octave_function *odessa_f; static octave_function *odessa_j; static octave_function *odessa_b; + static octave_function *odessa_fj; // Have we warned about imaginary values returned from user function? static bool warned_fcn_imaginary = false; *************** *** 230,235 **** --- 231,547 ---- return retval; } + // Storage for the cached values of the function and jacobian + static ColumnVector cached_fcn; + static Matrix cached_jac; + static ColumnVector cached_b; + static ColumnVector cached_x; + static double cached_t; + static ColumnVector cached_theta; + static int cached_column; + + static ColumnVector + odessa_user_f_cached (const ColumnVector& x, double t, + const ColumnVector& theta, int column) + { + ColumnVector retval; + + if (cached_x == x && cached_t == t && cached_theta == theta) + retval = cached_fcn; + else + { + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + if (column > 0) + args(3) = static_cast (column); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + args(1) = t; + + if (n > 1) + args(0) = x; + else if (n == 1) + args(0) = x(0); + else + args(0) = Matrix (); + + if (odessa_fj) + { + octave_value_list tmp; + if (column > 0) + tmp = odessa_fj->do_multi_index_op (3, args); + else + tmp = odessa_fj->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + retval = cached_fcn; + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + if (tmp(2).is_defined ()) + { + if (! warned_b_imaginary && tmp(2).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied inhomogeneity function"); + warned_b_imaginary = true; + } + + cached_b = ColumnVector (tmp(2).vector_value ()); + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + + cached_x = x; + cached_t = t; + cached_theta = theta; + cached_column = column; + } + else + gripe_user_supplied_eval ("odessa"); + } + } + + return retval; + } + + static Matrix + odessa_user_j_cached (const ColumnVector& x, double t, + const ColumnVector& theta, int column) + { + Matrix retval; + + if (cached_x == x && cached_t == t && cached_theta == theta) + retval = cached_jac; + else + { + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + if (column > 0) + args(3) = static_cast (column); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + args(1) = t; + + if (n > 1) + args(0) = x; + else if (n == 1) + args(0) = x(0); + else + args(0) = Matrix (); + + if (odessa_fj) + { + octave_value_list tmp; + if (column > 0) + tmp = odessa_fj->do_multi_index_op (3, args); + else + tmp = odessa_fj->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + retval = cached_jac; + } + + if (tmp(2).is_defined ()) + { + if (! warned_b_imaginary && tmp(2).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied inhomogeneity function"); + warned_b_imaginary = true; + } + + cached_b = ColumnVector (tmp(2).vector_value ()); + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + + cached_x = x; + cached_t = t; + cached_theta = theta; + cached_column = column; + } + else + gripe_user_supplied_eval ("odessa"); + } + } + + return retval; + } + + static ColumnVector + odessa_user_b_cached (const ColumnVector& x, double t, + const ColumnVector& theta, int column) + { + ColumnVector retval; + + if (cached_x == x && cached_t == t && cached_theta == theta && + cached_column == column) + retval = cached_b; + else if (column < 1) + gripe_user_supplied_eval ("odessa"); + else + { + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + args(3) = static_cast (column); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + args(1) = t; + + if (n > 1) + args(0) = x; + else if (n == 1) + args(0) = x(0); + else + args(0) = Matrix (); + + if (odessa_fj) + { + octave_value_list tmp; + if (column > 0) + tmp = odessa_fj->do_multi_index_op (3, args); + else + tmp = odessa_fj->do_multi_index_op (2, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 1) + { + if (tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + + cached_fcn = ColumnVector (tmp(0).vector_value ()); + } + + if (tmp(1).is_defined ()) + { + if (! warned_jac_imaginary && tmp(1).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + cached_jac = tmp(1).matrix_value (); + } + + if (tmp(2).is_defined ()) + { + if (! warned_b_imaginary && tmp(2).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied inhomogeneity function"); + warned_b_imaginary = true; + } + + cached_b = ColumnVector (tmp(2).vector_value ()); + retval = cached_b; + } + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + + cached_x = x; + cached_t = t; + cached_theta = theta; + cached_column = column; + } + else + gripe_user_supplied_eval ("odessa"); + } + } + + return retval; + } + static octave_value make_list (const Array& m_array) { *************** *** 319,327 **** with each element of the list corresponding to an element of the\n\ vector @var{t}.\n\ \n\ ! The first argument, @var{fcn}, is a string that names the function to\n\ ! call to compute the vector of right hand sides for the set of equations.\n\ ! The function must have the form\n\ \n\ @example\n\ @var{xdot} = f (@var{x}, @var{t}, @var{p})\n\ --- 631,639 ---- with each element of the list corresponding to an element of the\n\ vector @var{t}.\n\ \n\ ! The first argument, @var{fcn}, is either a string that names the function,\n\ ! a function handle, or an inline function, to call to compute the vector of\n\ ! right hand sides for the set of equations. The function must have the form\n\ \n\ @example\n\ @var{xdot} = f (@var{x}, @var{t}, @var{p})\n\ *************** *** 332,345 **** \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\ ! [\"f\"; \"j\"; \"b\"]\n\ ! @end example\n\ ! \n\ ! in which the first element names the function @math{f} described\n\ ! above, the second element names a function to compute the Jacobian\n\ of @math{f}, and the third element names a function to compute the\n\ inhomogeneity matrix.\n\ \n\ --- 644,652 ---- \n\ The variable @var{p} is a vector of parameters.\n\ \n\ ! The function @var{fcn} can also return a second argument or even a third\n\ ! argument, in which case the first element names the function @math{f}\n\ ! described above, the second element names a function to compute the Jacobian\n\ of @math{f}, and the third element names a function to compute the\n\ inhomogeneity matrix.\n\ \n\ *************** *** 382,387 **** --- 689,700 ---- dp_2\n\ @end example\n\ \n\ + The @var{fcn} argument may also be an array of strings\n\ + \n\ + @example\n\ + [\"f\"; \"j\"; \"b\"]\n\ + @end example\n\ + \n\ The second argument, @var{x_0}, specifies the intial state of the system.\n\ \n\ The third argument, @var{p}, specifies the set of parameters.\n\ *************** *** 428,480 **** return retval; } odessa_f = 0; odessa_j = 0; odessa_b = 0; octave_value f_arg = args(0); ! int nr = f_arg.rows (); ! ! if (nr == 1) { ! 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) { ! odessa_f = extract_function ! (tmp(0), "odessa", "__odessa_fcn__", ! "function xdot = __odessa_fcn__ (x, t, p) xdot = ", ! "; endfunction"); ! ! if (odessa_f) { ! 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(2), "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; } } } --- 741,843 ---- return retval; } + std::string fcn_name, fname, jac_name, jname, inhom_name, iname; odessa_f = 0; odessa_j = 0; odessa_b = 0; + odessa_fj = 0; octave_value f_arg = args(0); ! if (f_arg.is_function_handle () || f_arg.is_inline_function ()) { ! if (odessa_opts.jacobian ()) ! odessa_fj = f_arg.function_value (); ! else ! odessa_f = f_arg.function_value (); } ! else { ! switch (f_arg.rows ()) { ! case 1: ! if (odessa_opts.jacobian ()) { ! fcn_name = unique_symbol_name ("__odessa_fcn_jac__"); ! if (odessa_opts.inhomogeneity ()) { ! fname = "function [y, jac, dxdp] = "; ! fname.append (fcn_name); ! fname.append (" (x, t, p, c) [y, jac] = "); } ! else ! { ! fname = "function [y, jac] = "; ! fname.append (fcn_name); ! fname.append (" (x, t, p) [y, jac] = "); ! } ! odessa_fj = extract_function ! (f_arg, "odessa", fcn_name, fname, "; endfunction"); } + else + { + fcn_name = unique_symbol_name ("__odessa_fcn__"); + fname = "function y = "; + fname.append (fcn_name); + fname.append (" (x, t, p) y = "); + odessa_f = extract_function + (f_arg, "odessa", fcn_name, fname, "; endfunction"); + } + break; + + case 2: + case 3: + { + string_vector tmp = args(0).all_strings (); + + if (! error_state) + { + fcn_name = unique_symbol_name ("__odessa_fcn__"); + fname = "function y = "; + fname.append (fcn_name); + fname.append (" (x, t, p) y = "); + odessa_f = extract_function + (tmp(0), "odessa", fcn_name, fname, "; endfunction"); + + if (odessa_f) + { + jac_name = unique_symbol_name ("__odessa_jac__"); + jname = "function jac = "; + jname.append(jac_name); + jname.append (" (x, t, p) jac = "); + odessa_j = extract_function + (tmp(1), "odessa", jac_name, jname, "; endfunction"); + + if (odessa_j && f_arg.rows () == 3) + { + inhom_name = + unique_symbol_name ("__odessa_inhom__"); + iname = "function dxdp = "; + iname.append(inhom_name); + iname.append (" (x, t, p, c) dxdp = "); + odessa_b = extract_function + (tmp(2), "odessa", inhom_name, iname, + "; endfunction"); + + if (! odessa_b) + odessa_j = 0; + } + + if (! odessa_j) + odessa_f = 0; + } + } + } + break; + + default: + ODESSA_ABORT1 + ("first arg should be a string or 2-element string array"); } } *************** *** 543,555 **** crit_times_set = true; } ! ODESFunc func (odessa_user_f); ! ! if (odessa_j) ! func.set_jsub_function (odessa_user_j); ! if (odessa_b) ! func.set_bsub_function (odessa_user_b); double tzero = out_times (0); --- 906,928 ---- crit_times_set = true; } ! ODESFunc func; ! if (odessa_f) ! { ! func.set_fsub_function (odessa_user_f); ! if (odessa_j) ! func.set_jsub_function (odessa_user_j); ! if (odessa_b) ! func.set_bsub_function (odessa_user_b); ! } ! else ! { ! func.set_fsub_function_cached (odessa_user_f_cached); ! func.set_jsub_function_cached (odessa_user_j_cached); ! if (odessa_opts.inhomogeneity ()) ! func.set_bsub_function_cached (odessa_user_b_cached); ! } double tzero = out_times (0); *************** *** 566,571 **** --- 939,951 ---- else output = ode.integrate (out_times); + if (fcn_name.length()) + clear_function (fcn_name); + if (jac_name.length()) + clear_function (jac_name); + if (inhom_name.length()) + clear_function (inhom_name); + if (! error_state) { int k = have_parameters ? 3 : 2; *************** *** 594,599 **** --- 974,982 ---- } } + // Destory cached values to prevent them being reused + cached_x = ColumnVector (); + unwind_protect::run_frame ("Fodessa"); return retval; *** ./src/ov-fcn-inline.h.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/ov-fcn-inline.h 2004-08-06 12:11:13.000000000 +0200 *************** *** 57,62 **** --- 57,64 ---- ~octave_fcn_inline (void) { } + bool is_inline_function (void) const { return true; } + octave_fcn_inline *fcn_inline_value (bool = false) { return this; } std::string fcn_text (void) const { return iftext; } *** ./src/ov-base.h.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/ov-base.h 2004-08-06 12:11:13.000000000 +0200 *************** *** 171,176 **** --- 171,178 ---- bool is_function_handle (void) const { return false; } + bool is_inline_function (void) const { return false; } + bool is_function (void) const { return false; } bool is_builtin_function (void) const { return false; } *** ./src/ov.h.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/ov.h 2004-08-06 12:11:13.000000000 +0200 *************** *** 484,489 **** --- 484,492 ---- virtual bool is_function_handle (void) const { return rep->is_function_handle (); } + virtual bool is_inline_function (void) const + { return rep->is_function_handle (); } + virtual bool is_function (void) const { return rep->is_function (); } *** ./src/variables.h.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/variables.h 2004-08-06 12:11:13.000000000 +0200 *************** *** 77,82 **** --- 77,85 ---- extern int symbol_exist (const std::string& name, const std::string& type = "any"); + extern std::string + unique_symbol_name (const std::string& basename); + extern bool lookup (symbol_record *s, bool exec_script = true); extern symbol_record * *************** *** 116,121 **** --- 119,126 ---- extern void munlock (const std::string&); extern bool mislocked (const std::string&); + extern bool clear_function (const std::string& nm); + // Symbol table for symbols at the top level. extern symbol_table *top_level_sym_tab; *** ./src/variables.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/variables.cc 2004-08-06 12:11:13.000000000 +0200 *************** *** 662,667 **** --- 662,677 ---- return retval; } + std::string + unique_symbol_name (const std::string& basename) + { + // XXX FIXME XXX Can we be smarter than just adding characters? + std::string name = basename; + while (symbol_exist (name, "any")) + name.append ("X"); + return name; + } + DEFUN (exist, args, , "-*- texinfo -*-\n\ @deftypefn {Built-in Function} {} exist (@var{name}, @var{type})\n\ *************** *** 1897,1902 **** --- 1907,1919 ---- } \ while (0) + + bool + clear_function (const std::string& nm) + { + return do_clear_function (nm); + } + DEFCMD (clear, args, , "-*- texinfo -*-\n\ @deffn {Command} clear [-x] pattern @dots{}\n\ *** ./src/ov-fcn-inline.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/ov-fcn-inline.cc 2004-08-06 14:17:04.000000000 +0200 *************** *** 55,65 **** : octave_fcn_handle (0, n), iftext (f), ifargs (a) { // Find a function name that isn't already in the symbol table. ! ! std::string fname = "__inline__"; ! ! while (symbol_exist (fname)) ! fname.append ("X"); // Form a string representing the function. --- 55,61 ---- : octave_fcn_handle (0, n), iftext (f), ifargs (a) { // Find a function name that isn't already in the symbol table. ! std::string fname = unique_symbol_name ("__inline__"); // Form a string representing the function. *************** *** 91,100 **** { fcn = tmp; ! // XXX FIXME XXX -- probably shouldn't be directly altering the ! // symbol table here. ! ! fbi_sym_tab->clear_function (fname); } else error ("inline: unable to define function"); --- 87,93 ---- { fcn = tmp; ! clear_function (fname); } else error ("inline: unable to define function"); *** ./src/dynamic-ld.cc.orig 2004-08-06 12:09:55.000000000 +0200 --- ./src/dynamic-ld.cc 2004-08-06 12:11:13.000000000 +0200 *************** *** 195,201 **** } static ! void clear_function (const std::string& fcn_name) { if (Vwarn_reload_forces_clear) warning (" %s", fcn_name.c_str ()); --- 195,201 ---- } static ! void do_clear_function (const std::string& fcn_name) { if (Vwarn_reload_forces_clear) warning (" %s", fcn_name.c_str ()); *************** *** 235,241 **** warning ("reloading %s clears the following functions:", oct_file.file_name().c_str ()); ! oct_file.close (clear_function); function = 0; } --- 235,241 ---- warning ("reloading %s clears the following functions:", oct_file.file_name().c_str ()); ! oct_file.close (do_clear_function); function = 0; }