Mercurial > octave-nkf
changeset 4049:a35a3c5d4740
[project @ 2002-08-16 08:54:31 by jwe]
| author | jwe |
|---|---|
| date | Fri, 16 Aug 2002 08:54:31 +0000 |
| parents | c9991c59cf5c |
| children | 6481f41a79f3 |
| files | liboctave/ChangeLog liboctave/DAEFunc.h liboctave/DAERTFunc.h liboctave/DASPK-opts.in liboctave/DASPK.cc liboctave/DASPK.h liboctave/DASRT-opts.in liboctave/DASRT.cc liboctave/DASRT.h liboctave/DASSL-opts.in liboctave/DASSL.cc liboctave/DASSL.h liboctave/LSODE-opts.in liboctave/LSODE.cc liboctave/LSODE.h liboctave/ODEFunc.h liboctave/ODESSA-opts.in liboctave/ODESSA.h liboctave/base-de.h mk-opts.pl |
| diffstat | 20 files changed, 766 insertions(+), 747 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/ChangeLog Fri Aug 16 08:54:31 2002 +0000 @@ -1,5 +1,29 @@ +2002-08-16 John W. Eaton <jwe@bevo.che.wisc.edu> + + * LSODE.h (rel_tol, abs_tol, px, pabs_tol, piwork, prwork, itol): + New data members. + (LSODE::sanity_checked): Delete unused data member. + + * DASPKL.h (initialized, abs_tol, rel_tol, px, pxdot, pabs_tol, + prel_tol, pinfo, piwork, prwork): New data members. + * DASSL.h (DASSL): Likewise. + + * DASRT.h (DASRT::sanity_checked): Delete unused data member. + + * DASRT.cc (DASRT::integrate (double)): Better handling of + initialization, changes in options, etc. + * DASPK.cc (DASPK::do_integrate): Likewise. + * DASSL.cc (DASSL::do_integrate): Likewise. + * LSODE.cc (LSODE::do_integrate): Likewise. + 2002-08-15 John W. Eaton <jwe@bevo.che.wisc.edu> + * DAEFunc.h (DAEFunc::reset): New data member. + * DAERTFunc.h (DAERTFunc::reset): Likewise. + + * base-de.h (base_diff_eqn::set_stop_time): Force restart here. + (base_diff_eqn::clear_stop_time): Likewise. + * DASSL.cc (DASSL::do_integrate (double)): Handle more optoins. * DASPK.cc (DASPK::do_integrate (double)): Likewise.
--- a/liboctave/DAEFunc.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DAEFunc.h Fri Aug 16 08:54:31 2002 +0000 @@ -44,16 +44,16 @@ double t, double cj); DAEFunc (void) - : fun (0), jac (0) { } + : fun (0), jac (0), reset (true) { } DAEFunc (DAERHSFunc f) - : fun (f), jac (0) { } + : fun (f), jac (0), reset (true) { } DAEFunc (DAERHSFunc f, DAEJacFunc j) - : fun (f), jac (j) { } + : fun (f), jac (j), reset (true) { } DAEFunc (const DAEFunc& a) - : fun (a.fun), jac (a.jac) { } + : fun (a.fun), jac (a.jac), reset (a.reset) { } DAEFunc& operator = (const DAEFunc& a) { @@ -61,6 +61,7 @@ { fun = a.fun; jac = a.jac; + reset = a.reset; } return *this; } @@ -72,6 +73,7 @@ DAEFunc& set_function (DAERHSFunc f) { fun = f; + reset = true; return *this; } @@ -80,6 +82,7 @@ DAEFunc& set_jacobian_function (DAEJacFunc j) { jac = j; + reset = true; return *this; } @@ -87,6 +90,13 @@ DAERHSFunc fun; DAEJacFunc jac; + + // This variable is TRUE when this object is constructed, and also + // after any internal data has changed. Derived classes may use + // this information (and change it) to know when to (re)initialize + // their own internal data related to this object. + + bool reset; }; #endif
--- a/liboctave/DAERTFunc.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DAERTFunc.h Fri Aug 16 08:54:31 2002 +0000 @@ -33,22 +33,22 @@ typedef ColumnVector (*DAERTConstrFunc) (const ColumnVector& x, double t); DAERTFunc (void) - : DAEFunc (), constr (0) { } + : DAEFunc (), constr (0), reset (true) { } DAERTFunc (DAERHSFunc f) - : DAEFunc (f), constr (0) { } + : DAEFunc (f), constr (0), reset (true) { } DAERTFunc (DAERHSFunc f, DAEJacFunc j) - : DAEFunc (f, j), constr (0) { } + : DAEFunc (f, j), constr (0), reset (true) { } DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf) - : DAEFunc (f), constr (cf) { } + : DAEFunc (f), constr (cf), reset (true) { } DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf, DAEJacFunc j) - : DAEFunc (f, j), constr (cf) { } + : DAEFunc (f, j), constr (cf), reset (true) { } DAERTFunc (const DAERTFunc& a) - : DAEFunc (a), constr (a.constr) { } + : DAEFunc (a), constr (a.constr), reset (a.reset) { } DAERTFunc& operator = (const DAERTFunc& a) { @@ -56,6 +56,7 @@ { DAEFunc::operator = (a); constr = a.constr; + reset = a.reset; } return *this; } @@ -67,12 +68,20 @@ DAERTFunc& set_constraint_function (DAERTConstrFunc cf) { constr = cf; + reset = true; return *this; } protected: DAERTConstrFunc constr; + + // This variable is TRUE when this object is constructed, and also + // after any internal data has changed. Derived classes may use + // this information (and change it) to know when to (re)initialize + // their own internal data related to this object. + + bool reset; }; #endif
--- a/liboctave/DASPK-opts.in Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASPK-opts.in Fri Aug 16 08:54:31 2002 +0000 @@ -15,10 +15,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -35,10 +36,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -62,10 +64,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = val; + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -89,10 +92,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = val; + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -153,4 +157,3 @@ INIT_VALUE = "0" SET_EXPR = "val" END_OPTION -
--- a/liboctave/DASPK.cc Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASPK.cc Fri Aug 16 08:54:31 2002 +0000 @@ -64,47 +64,7 @@ static DAEFunc::DAEJacFunc user_jac; static int nn; -DASPK::DASPK (void) : DAE () -{ - sanity_checked = 0; - - info.resize (20); - - for (int i = 0; i < 20; i++) - info.elem (i) = 0; -} - -DASPK::DASPK (const ColumnVector& state, double time, DAEFunc& f) - : DAE (state, time, f) -{ - n = size (); - - sanity_checked = 0; - - info.resize (20); - - for (int i = 0; i < 20; i++) - info.elem (i) = 0; -} - -DASPK::DASPK (const ColumnVector& state, const ColumnVector& deriv, - double time, DAEFunc& f) - : DAE (state, deriv, time, f) -{ - n = size (); - - DAEFunc::set_function (f.function ()); - DAEFunc::set_jacobian_function (f.jacobian_function ()); - - sanity_checked = 0; - - info.resize (20); - - for (int i = 0; i < 20; i++) - info.elem (i) = 0; -} - -int +static int ddaspk_f (const double& time, const double *state, const double *deriv, const double&, double *delta, int& ires, double *, int *) { @@ -137,7 +97,7 @@ //NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT, //C WP, IWP, B, EPLIN, IER, RPAR, IPAR) -int +static int ddaspk_psol (const int& neq, const double& time, const double *state, const double *deriv, const double *savr, const double& cj, const double *wght, double *wp, @@ -149,7 +109,7 @@ } -int +static int ddaspk_j (const double& time, const double *state, const double *deriv, double *pd, const double& cj, double *, int *) { @@ -181,178 +141,220 @@ ColumnVector retval; - if (restart) + if (! initialized || restart || DAEFunc::reset|| DASPK_options::reset) { - restart = false; - info.elem (0) = 0; - } + integration_error = false; + + initialized = true; + + info.resize (20); - liw = 40 + n; - if (info(9) == 1 || info(9) == 3) - liw += n; - if (info (10) == 1 || info(15) == 1) - liw += n; + for (int i = 0; i < 20; i++) + info(i) = 0; + + pinfo = info.fortran_vec (); - lrw = 50 + 9*n; - if (info(5) == 0) - lrw += n*n; - if (info(15) == 1) - lrw += n; + int n = size (); - if (iwork.length () != liw) - iwork.resize (liw); + nn = n; + + info(0) = 0; - if (rwork.length () != lrw) - rwork.resize (lrw); - - integration_error = false; + if (stop_time_set) + { + rwork(0) = stop_time; + info(3) = 1; + } + else + info(3) = 0; - if (DAEFunc::jacobian_function ()) - info.elem (4) = 1; - else - info.elem (4) = 0; + px = x.fortran_vec (); + pxdot = xdot.fortran_vec (); + + // DAEFunc + + user_fun = DAEFunc::function (); + user_jac = DAEFunc::jacobian_function (); - double *px = x.fortran_vec (); - double *pxdot = xdot.fortran_vec (); + if (user_fun) + { + int ires = 0; - nn = n; - user_fun = DAEFunc::fun; - user_jac = DAEFunc::jac; + ColumnVector res = (*user_fun) (x, xdot, t, ires); - if (! sanity_checked) - { - int ires = 0; + if (res.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("daspk: inconsistent sizes for state and residual vectors"); - ColumnVector res = (*user_fun) (x, xdot, t, ires); - - if (res.length () != x.length ()) + integration_error = true; + return retval; + } + } + else { (*current_liboctave_error_handler) - ("daspk: inconsistent sizes for state and residual vectors"); + ("daspk: no user supplied RHS subroutine!"); + + integration_error = true; + return retval; + } + + info(4) = user_jac ? 1 : 0; + + DAEFunc::reset = false; + + // DASPK_options + + Array<double> abs_tol = absolute_tolerance (); + Array<double> rel_tol = relative_tolerance (); + + int abs_tol_len = abs_tol.length (); + int rel_tol_len = rel_tol.length (); + + if (abs_tol_len == 1 && rel_tol_len == 1) + { + info(1) = 0; + } + else if (abs_tol_len == n && rel_tol_len == n) + { + info(1) = 1; + } + else + { + (*current_liboctave_error_handler) + ("daspk: inconsistent sizes for tolerance arrays"); integration_error = true; return retval; } - sanity_checked = 1; - } - - if (stop_time_set) - { - rwork.elem (0) = stop_time; - info.elem (3) = 1; - } - else - info.elem (3) = 0; - - Array<double> abs_tol = absolute_tolerance (); - Array<double> rel_tol = relative_tolerance (); - - int abs_tol_len = abs_tol.length (); - int rel_tol_len = rel_tol.length (); - - if (abs_tol_len == 1 && rel_tol_len == 1) - { - info.elem (1) = 0; - } - else if (abs_tol_len == n && rel_tol_len == n) - { - info.elem (1) = 1; - } - else - { - (*current_liboctave_error_handler) - ("daspk: inconsistent sizes for tolerance arrays"); + double hmax = maximum_step_size (); + if (hmax >= 0.0) + { + rwork(1) = hmax; + info(6) = 1; + } + else + info(6) = 0; - integration_error = true; - return retval; - } - - double hmax = maximum_step_size (); - if (hmax >= 0.0) - { - rwork.elem (1) = hmax; - info.elem (6) = 1; - } - else - info.elem (6) = 0; - - double h0 = initial_step_size (); - if (h0 >= 0.0) - { - rwork.elem (2) = h0; - info.elem (7) = 1; - } - else - info.elem (7) = 0; - - int maxord = maximum_order (); - if (maxord >= 0) - { - if (maxord > 0 && maxord < 6) + double h0 = initial_step_size (); + if (h0 >= 0.0) { - info(8) = 1; - iwork(2) = maxord; + rwork(2) = h0; + info(7) = 1; } else + info(7) = 0; + + int maxord = maximum_order (); + if (maxord >= 0) { + if (maxord > 0 && maxord < 6) + { + info(8) = 1; + iwork(2) = maxord; + } + else + { + (*current_liboctave_error_handler) + ("daspk: invalid value for maximum order"); + integration_error = true; + return retval; + } + } + + int eiq = enforce_inequality_constraints (); + switch (eiq) + { + case 1: + case 3: + { + Array<int> ict = inequality_constraint_types (); + + if (ict.length () == n) + { + for (int i = 0; i < n; i++) + { + int val = ict(i); + if (val < -2 || val > 2) + { + (*current_liboctave_error_handler) + ("daspk: invalid value for inequality constraint type"); + integration_error = true; + return retval; + } + iwork(40+i) = val; + } + } + else + { + (*current_liboctave_error_handler) + ("daspk: inequality constraint types size mismatch"); + integration_error = true; + return retval; + } + } + // Fall through... + + case 2: + info(9) = eiq; + break; + + default: (*current_liboctave_error_handler) - ("daspk: invalid value for maximum order"); + ("daspk: invalid value for enforce inequality constraints option"); integration_error = true; return retval; } - } + + int ccic = compute_consistent_initial_condition (); + if (ccic) + { + if (ccic == 1) + { + // XXX FIXME XXX -- this code is duplicated below. + + Array<int> av = algebraic_variables (); - int eiq = enforce_inequality_constraints (); - switch (eiq) - { - case 1: - case 3: - { - Array<int> ict = inequality_constraint_types (); + if (av.length () == n) + { + int lid; + if (eiq == 0 || eiq == 2) + lid = 40; + else if (eiq == 1 || eiq == 3) + lid = 40 + n; + else + abort (); - if (ict.length () == n) - { - for (int i = 0; i < n; i++) - { - int val = ict(i); - if (val < -2 || val > 2) - { - (*current_liboctave_error_handler) - ("daspk: invalid value for inequality constraint type"); - integration_error = true; - return retval; - } - iwork(40+i) = val; - } - } - else - { - (*current_liboctave_error_handler) - ("daspk: inequality constraint types size mismatch"); - integration_error = true; - return retval; - } - } - // Fall through... + for (int i = 0; i < n; i++) + iwork(lid+i) = av(i) ? -1 : 1; + } + else + { + (*current_liboctave_error_handler) + ("daspk: algebraic variables size mismatch"); + integration_error = true; + return retval; + } + } + else if (ccic != 2) + { + (*current_liboctave_error_handler) + ("daspk: invalid value for compute consistent initial condition option"); + integration_error = true; + return retval; + } - case 2: - info(9) = eiq; - break; + info(10) = ccic; + } - default: - (*current_liboctave_error_handler) - ("daspk: invalid value for enforce inequality constraints option"); - integration_error = true; - return retval; - } + int eavfet = exclude_algebraic_variables_from_error_test (); + if (eavfet) + { + info(15) = 1; - int ccic = compute_consistent_initial_condition (); - if (ccic) - { - if (ccic == 1) - { - // XXX FIXME XXX -- this code is duplicated below. + // XXX FIXME XXX -- this code is duplicated above. Array<int> av = algebraic_variables (); @@ -369,102 +371,77 @@ for (int i = 0; i < n; i++) iwork(lid+i) = av(i) ? -1 : 1; } + } + + if (use_initial_condition_heuristics ()) + { + Array<double> ich = initial_condition_heuristics (); + + if (ich.length () == 6) + { + iwork(31) = NINT (ich(0)); + iwork(32) = NINT (ich(1)); + iwork(33) = NINT (ich(2)); + iwork(34) = NINT (ich(3)); + + rwork(13) = ich(4); + rwork(14) = ich(5); + } else { (*current_liboctave_error_handler) - ("daspk: algebraic variables size mismatch"); + ("daspk: invalid initial condition heuristics option"); integration_error = true; return retval; } + + info(16) = 1; } - else if (ccic != 2) + + int pici = print_initial_condition_info (); + switch (pici) { + case 0: + case 1: + case 2: + info(17) = pici; + break; + + default: (*current_liboctave_error_handler) - ("daspk: invalid value for compute consistent initial condition option"); + ("daspk: invalid value for print initial condition info option"); integration_error = true; return retval; + break; } - info(10) = ccic; - } + DASPK_options::reset = false; - if (exclude_algebraic_variables_from_error_test ()) - { - info(15) = 1; - - // XXX FIXME XXX -- this code is duplicated above. - - Array<int> av = algebraic_variables (); + liw = 40 + n; + if (eiq == 1 || eiq == 3) + liw += n; + if (ccic == 1 || eavfet == 1) + liw += n; - if (av.length () == n) - { - int lid; - if (eiq == 0 || eiq == 2) - lid = 40; - else if (eiq == 1 || eiq == 3) - lid = 40 + n; - else - abort (); + lrw = 50 + 9*n; + if (! user_jac) + lrw += n*n; + if (eavfet == 1) + lrw += n; - for (int i = 0; i < n; i++) - iwork(lid+i) = av(i) ? -1 : 1; - } + iwork.resize (liw); + rwork.resize (lrw); + + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + + restart = false; } - if (use_initial_condition_heuristics ()) - { - Array<double> ich = initial_condition_heuristics (); - - if (ich.length () == 6) - { - iwork(31) = NINT (ich(0)); - iwork(32) = NINT (ich(1)); - iwork(33) = NINT (ich(2)); - iwork(34) = NINT (ich(3)); - - rwork(13) = ich(4); - rwork(14) = ich(5); - } - else - { - (*current_liboctave_error_handler) - ("daspk: invalid initial condition heuristics option"); - integration_error = true; - return retval; - } - - info(16) = 1; - } + static double *dummy = 0; + static int *idummy = 0; - int pici = print_initial_condition_info (); - switch (pici) - { - case 0: - case 1: - case 2: - info(17) = pici; - break; - - default: - (*current_liboctave_error_handler) - ("daspk: invalid value for print initial condition info option"); - integration_error = true; - return retval; - break; - } - - double *dummy = 0; - int *idummy = 0; - - int *pinfo = info.fortran_vec (); - int *piwork = iwork.fortran_vec (); - double *prwork = rwork.fortran_vec (); - double *pabs_tol = abs_tol.fortran_vec (); - double *prel_tol = rel_tol.fortran_vec (); - -// again: - - F77_XFCN (ddaspk, DDASPK, (ddaspk_f, n, t, px, pxdot, tout, pinfo, + F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo, prel_tol, pabs_tol, istate, prwork, lrw, piwork, liw, dummy, idummy, ddaspk_j, ddaspk_psol)); @@ -546,7 +523,9 @@ DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) { @@ -589,7 +568,9 @@ const ColumnVector& tcrit) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) {
--- a/liboctave/DASPK.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASPK.h Fri Aug 16 08:54:31 2002 +0000 @@ -37,12 +37,14 @@ { public: - DASPK (void); + DASPK (void) : DAE (), DASPK_options (), initialized (false) { } - DASPK (const ColumnVector& x, double time, DAEFunc& f); + DASPK (const ColumnVector& state, double time, DAEFunc& f) + : DAE (state, time, f), DASPK_options (), initialized (false) { } - DASPK (const ColumnVector& x, const ColumnVector& xdot, - double time, DAEFunc& f); + DASPK (const ColumnVector& state, const ColumnVector& deriv, + double time, DAEFunc& f) + : DAE (state, deriv, time, f), DASPK_options (), initialized (false) { } ~DASPK (void) { } @@ -61,20 +63,26 @@ private: - int n; + bool initialized; + int liw; int lrw; - int sanity_checked; + Array<int> info; Array<int> iwork; + Array<double> rwork; - friend int ddaspk_j (double *time, double *state, double *deriv, - double *pd, double *cj, double *rpar, int *ipar); + Array<double> abs_tol; + Array<double> rel_tol; - friend int ddaspk_f (double *time, double *state, double *deriv, - double *delta, int *ires, double *rpar, int *ipar); - + double *px; + double *pxdot; + double *pabs_tol; + double *prel_tol; + int *pinfo; + int *piwork; + double *prwork; }; #endif
--- a/liboctave/DASRT-opts.in Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASRT-opts.in Fri Aug 16 08:54:31 2002 +0000 @@ -15,10 +15,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -35,10 +36,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -50,6 +52,13 @@ END_OPTION OPTION + NAME = "maximum order" + TYPE = "int" + INIT_VALUE = "-1" + SET_EXPR = "val" +END_OPTION + +OPTION NAME = "maximum step size" TYPE = "double" INIT_VALUE = "-1.0" @@ -57,13 +66,6 @@ END_OPTION OPTION - NAME = "minimum step size" - TYPE = "double" - INIT_VALUE = "0.0" - SET_EXPR = "(val >= 0.0) ? val : 0.0" -END_OPTION - -OPTION NAME = "step limit" TYPE = "int" INIT_VALUE = "-1"
--- a/liboctave/DASRT.cc Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASRT.cc Fri Aug 16 08:54:31 2002 +0000 @@ -48,10 +48,6 @@ #include "f77-fcn.h" #include "lo-error.h" -#ifndef F77_FUNC -#define F77_FUNC(x, X) F77_FCN (x, X) -#endif - typedef int (*dasrt_fcn_ptr) (const double&, const double*, const double*, double*, int&, double*, int*); @@ -142,102 +138,63 @@ return 0; } -DASRT::DASRT (void) - : DAERT () -{ - initialized = false; - - sanity_checked = false; - - info.resize (30, 0); - - ng = 0; - - liw = 0; - lrw = 0; -} - -DASRT::DASRT (const ColumnVector& state, double time, DAERTFunc& f) - : DAERT (state, time, f) -{ - n = size (); - - initialized = false; - - liw = 20 + n; - lrw = 50 + 9*n + n*n; - - sanity_checked = false; - - info.resize (15, 0); - - DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); - - if (tmp_csub) - { - ColumnVector tmp = tmp_csub (state, time); - ng = tmp.length (); - } - else - ng = 0; - - jroot.resize (ng, 1); -} - -DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv, - double time, DAERTFunc& f) - : DAERT (state, deriv, time, f) -{ - n = size (); - - initialized = false; - - sanity_checked = false; - - info.resize (30, 0); - - DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); - - if (tmp_csub) - { - ColumnVector tmp = tmp_csub (state, time); - ng = tmp.length (); - } - else - ng = 0; - - liw = 20 + n + 1000; - lrw = 50 + 9*n + n*n + 1000; - - jroot.resize (ng, 1); -} - void DASRT::integrate (double tout) { DASRT_result retval; - if (! initialized) + // I suppose this is the safe thing to do. If this is the first + // call, or if anything about the problem has changed, we should + // start completely fresh. + + if (! initialized || restart + || DAEFunc::reset || DAERTFunc::reset || DASRT_options::reset) { - info(0) = 0; - integration_error = false; - user_fsub = DAEFunc::function (); - user_jsub = DAEFunc::jacobian_function (); - user_csub = DAERTFunc::constraint_function (); + initialized = true; + + info.resize (15); + + for (int i = 0; i < 15; i++) + info(i) = 0; + + pinfo = info.fortran_vec (); + + int n = size (); + + nn = n; - if (user_jsub) - info(4) = 1; + liw = 20 + n; + lrw = 50 + 9*n + n*n; + + iwork.resize (liw); + rwork.resize (lrw); + + info(0) = 0; + + if (stop_time_set) + { + info(3) = 1; + rwork(0) = stop_time; + } else - info(4) = 0; + info(3) = 0; px = x.fortran_vec (); pxdot = xdot.fortran_vec (); - nn = n; + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + + restart = false; - if (! sanity_checked) + // DAEFunc + + user_fsub = DAEFunc::function (); + user_jsub = DAEFunc::jacobian_function (); + + if (user_fsub) { int ires = 0; @@ -246,20 +203,85 @@ if (fval.length () != x.length ()) { (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for state and residual vectors"); + ("dasrt: inconsistent sizes for state and residual vectors"); integration_error = true; return; } + } + else + { + (*current_liboctave_error_handler) + ("dasrt: no user supplied RHS subroutine!"); - sanity_checked = true; + integration_error = true; + return; + } + + info(4) = user_jsub ? 1 : 0; + + DAEFunc::reset = false; + + // DAERTFunc + + user_csub = DAERTFunc::constraint_function (); + + if (user_csub) + { + ColumnVector tmp = (*user_csub) (x, t); + ng = tmp.length (); } - - if (iwork.length () != liw) - iwork.resize (liw); + else + ng = 0; + + jroot.resize (ng, 1); + + pjroot = jroot.fortran_vec (); + + DAERTFunc::reset = false; + + // DASRT_options + + if (maximum_step_size () >= 0.0) + { + rwork(1) = maximum_step_size (); + info(6) = 1; + } + else + info(6) = 0; - if (rwork.length () != lrw) - rwork.resize (lrw); + if (initial_step_size () >= 0.0) + { + rwork(2) = initial_step_size (); + info(7) = 1; + } + else + info(7) = 0; + + int maxord = maximum_order (); + if (maxord >= 0) + { + if (maxord > 0 && maxord < 6) + { + info(8) = 1; + iwork(2) = maxord; + } + else + { + (*current_liboctave_error_handler) + ("dassl: invalid value for maximum order"); + integration_error = true; + return; + } + } + + if (step_limit () >= 0) + { + info(11) = 1; + iwork(18) = step_limit (); + } + else + info(11) = 0; abs_tol = absolute_tolerance (); rel_tol = relative_tolerance (); @@ -278,65 +300,22 @@ else { (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for tolerance arrays"); + ("dasrt: inconsistent sizes for tolerance arrays"); integration_error = true; return; } - if (initial_step_size () >= 0.0) - { - rwork(2) = initial_step_size (); - info(7) = 1; - } - else - info(7) = 0; - - if (step_limit () >= 0) - { - info(11) = 1; - iwork(18) = step_limit (); - } - else - info(11) = 0; - - if (maximum_step_size () >= 0.0) - { - rwork(1) = maximum_step_size (); - info(6) = 1; - } - else - info(6) = 0; - - pinfo = info.fortran_vec (); - piwork = iwork.fortran_vec (); pabs_tol = abs_tol.fortran_vec (); prel_tol = rel_tol.fortran_vec (); - prwork = rwork.fortran_vec (); - pjroot = jroot.fortran_vec (); - info(5) = 0; - info(8) = 0; - initialized = true; + DASRT_options::reset = false; } - if (restart) - { - info(0) = 0; + static double *dummy = 0; + static int *idummy = 0; - if (stop_time_set) - { - info(3) = 1; - rwork(0) = stop_time; - } - else - info(3) = 0; - } - - double *dummy = 0; - int *idummy = 0; - - F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, px, pxdot, tout, pinfo, + F77_XFCN (ddasrt, DASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo, prel_tol, pabs_tol, istate, prwork, lrw, piwork, liw, dummy, idummy, ddasrt_j, ddasrt_g, ng, pjroot)); @@ -344,7 +323,7 @@ if (f77_exception_encountered) { integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in dassl"); + (*current_liboctave_error_handler) ("unrecoverable error in dasrt"); } else { @@ -409,6 +388,7 @@ ColumnVector t_out = tout; int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) { @@ -467,6 +447,7 @@ ColumnVector t_outs = tout; int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) {
--- a/liboctave/DASRT.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASRT.h Fri Aug 16 08:54:31 2002 +0000 @@ -74,12 +74,14 @@ { public: - DASRT (void); + DASRT (void) : DAERT (), DASRT_options (), initialized (false) { } - DASRT (const ColumnVector& state, double time, DAERTFunc& f); + DASRT (const ColumnVector& state, double time, DAERTFunc& f) + : DAERT (state, time, f), DASRT_options (), initialized (false) { } DASRT (const ColumnVector& state, const ColumnVector& deriv, - double time, DAERTFunc& f); + double time, DAERTFunc& f) + : DAERT (state, deriv, time, f), DASRT_options (), initialized (false) { } ~DASRT (void) { } @@ -94,12 +96,9 @@ bool initialized; - bool sanity_checked; - int liw; int lrw; - int n; int ng; Array<int> info;
--- a/liboctave/DASSL-opts.in Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASSL-opts.in Fri Aug 16 08:54:31 2002 +0000 @@ -15,10 +15,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION @@ -35,10 +36,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION
--- a/liboctave/DASSL.cc Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASSL.cc Fri Aug 16 08:54:31 2002 +0000 @@ -56,53 +56,7 @@ static int nn; -DASSL::DASSL (void) : DAE () -{ - liw = 0; - lrw = 0; - - sanity_checked = false; - - info.resize (15); - - for (int i = 0; i < 15; i++) - info.elem (i) = 0; -} - -DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f) - : DAE (state, time, f) -{ - n = size (); - - liw = 20 + n; - lrw = 40 + 9*n + n*n; - - sanity_checked = false; - - info.resize (15, 0); -} - -DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv, - double time, DAEFunc& f) - : DAE (state, deriv, time, f) -{ - n = size (); - - DAEFunc::set_function (f.function ()); - DAEFunc::set_jacobian_function (f.jacobian_function ()); - - liw = 20 + n; - lrw = 40 + 9*n + n*n; - - sanity_checked = false; - - info.resize (15); - - for (int i = 0; i < 15; i++) - info.elem (i) = 0; -} - -int +static int ddassl_f (const double& time, const double *state, const double *deriv, double *delta, int& ires, double *, int *) { @@ -134,7 +88,7 @@ return 0; } -int +static int ddassl_j (const double& time, const double *state, const double *deriv, double *pd, const double& cj, double *, int *) { @@ -161,140 +115,158 @@ ColumnVector DASSL::do_integrate (double tout) { - // XXX FIXME XXX -- should handle all this option stuff just once - // for each new problem. - ColumnVector retval; - if (restart) + if (! initialized || restart || DAEFunc::reset|| DASSL_options::reset) { - restart = false; - info.elem (0) = 0; - } + integration_error = false; + + initialized = true; + + info.resize (15); + + for (int i = 0; i < 15; i++) + info(i) = 0; - if (iwork.length () != liw) - iwork.resize (liw); + pinfo = info.fortran_vec (); + + int n = size (); + + liw = 20 + n; + lrw = 40 + 9*n + n*n; - if (rwork.length () != lrw) - rwork.resize (lrw); + nn = n; - integration_error = false; + iwork.resize (liw); + rwork.resize (lrw); - user_fun = DAEFunc::fun; - user_jac = DAEFunc::jac; + info(0) = 0; - if (user_jac) - info.elem (4) = 1; - else - info.elem (4) = 0; + if (stop_time_set) + { + rwork(0) = stop_time; + info(3) = 1; + } + else + info(3) = 0; - double *px = x.fortran_vec (); - double *pxdot = xdot.fortran_vec (); + px = x.fortran_vec (); + pxdot = xdot.fortran_vec (); + + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + + restart = false; + + // DAEFunc - nn = n; + user_fun = DAEFunc::function (); + user_jac = DAEFunc::jacobian_function (); + + if (user_fun) + { + int ires = 0; + + ColumnVector res = (*user_fun) (x, xdot, t, ires); - if (! sanity_checked) - { - int ires = 0; + if (res.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("dassl: inconsistent sizes for state and residual vectors"); - ColumnVector res = (*user_fun) (x, xdot, t, ires); - - if (res.length () != x.length ()) + integration_error = true; + return retval; + } + } + else { (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for state and residual vectors"); + ("dassl: no user supplied RHS subroutine!"); integration_error = true; return retval; } - sanity_checked = true; - } + info(4) = user_jac ? 1 : 0; - if (stop_time_set) - { - rwork.elem (0) = stop_time; - info.elem (3) = 1; - } - else - info.elem (3) = 0; + DAEFunc::reset = false; - Array<double> abs_tol = absolute_tolerance (); - Array<double> rel_tol = relative_tolerance (); + // DASSL_options - int abs_tol_len = abs_tol.length (); - int rel_tol_len = rel_tol.length (); + double hmax = maximum_step_size (); + if (hmax >= 0.0) + { + rwork(1) = hmax; + info(6) = 1; + } + else + info(6) = 0; - if (abs_tol_len == 1 && rel_tol_len == 1) - { - info.elem (1) = 0; - } - else if (abs_tol_len == n && rel_tol_len == n) - { - info.elem (1) = 1; - } - else - { - (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for tolerance arrays"); + double h0 = initial_step_size (); + if (h0 >= 0.0) + { + rwork(2) = h0; + info(7) = 1; + } + else + info(7) = 0; - integration_error = true; - return retval; - } - - double hmax = maximum_step_size (); - if (hmax >= 0.0) - { - rwork.elem (1) = hmax; - info.elem (6) = 1; - } - else - info.elem (6) = 0; + int maxord = maximum_order (); + if (maxord >= 0) + { + if (maxord > 0 && maxord < 6) + { + info(8) = 1; + iwork(2) = maxord; + } + else + { + (*current_liboctave_error_handler) + ("dassl: invalid value for maximum order"); + integration_error = true; + return retval; + } + } - double h0 = initial_step_size (); - if (h0 >= 0.0) - { - rwork.elem (2) = h0; - info.elem (7) = 1; - } - else - info.elem (7) = 0; + int enc = enforce_nonnegativity_constraints (); + info(9) = enc ? 1 : 0; + + int ccic = compute_consistent_initial_condition (); + info(10) = ccic ? 1 : 0; + + abs_tol = absolute_tolerance (); + rel_tol = relative_tolerance (); - int maxord = maximum_order (); - if (maxord >= 0) - { - if (maxord > 0 && maxord < 6) + int abs_tol_len = abs_tol.length (); + int rel_tol_len = rel_tol.length (); + + if (abs_tol_len == 1 && rel_tol_len == 1) { - info(8) = 1; - iwork(2) = maxord; + info(1) = 0; + } + else if (abs_tol_len == n && rel_tol_len == n) + { + info(1) = 1; } else { (*current_liboctave_error_handler) - ("dassl: invalid value for maximum order"); + ("dassl: inconsistent sizes for tolerance arrays"); + integration_error = true; return retval; } + + pabs_tol = abs_tol.fortran_vec (); + prel_tol = rel_tol.fortran_vec (); + + DASSL_options::reset = false; } - int enc = enforce_nonnegativity_constraints (); - info (9) = enc ? 1 : 0; - - int ccic = compute_consistent_initial_condition (); - info(10) = ccic ? 1 : 0; - - double *dummy = 0; - int *idummy = 0; + static double *dummy = 0; + static int *idummy = 0; - int *pinfo = info.fortran_vec (); - int *piwork = iwork.fortran_vec (); - double *prwork = rwork.fortran_vec (); - double *pabs_tol = abs_tol.fortran_vec (); - double *prel_tol = rel_tol.fortran_vec (); - -// again: - - F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo, + F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo, prel_tol, pabs_tol, istate, prwork, lrw, piwork, liw, dummy, idummy, ddassl_j)); @@ -366,7 +338,9 @@ DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) { @@ -409,7 +383,9 @@ const ColumnVector& tcrit) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) {
--- a/liboctave/DASSL.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/DASSL.h Fri Aug 16 08:54:31 2002 +0000 @@ -37,12 +37,14 @@ { public: - DASSL (void); + DASSL (void) : DAE (), DASSL_options (), initialized (false) { } - DASSL (const ColumnVector& state, double time, DAEFunc& f); + DASSL (const ColumnVector& state, double time, DAEFunc& f) + : DAE (state, time, f), DASSL_options (), initialized (false) { } - DASSL (const ColumnVector& state, const ColumnVector& xdot, - double time, DAEFunc& f); + DASSL (const ColumnVector& state, const ColumnVector& deriv, + double time, DAEFunc& f) + : DAE (state, deriv, time, f), DASSL_options (), initialized (false) { } ~DASSL (void) { } @@ -61,20 +63,26 @@ private: - int n; + bool initialized; + int liw; int lrw; - bool sanity_checked; + Array<int> info; Array<int> iwork; + Array<double> rwork; - friend int ddassl_j (double *time, double *state, double *deriv, - double *pd, double *cj, double *rpar, int *ipar); + Array<double> abs_tol; + Array<double> rel_tol; - friend int ddassl_f (double *time, double *state, double *deriv, - double *delta, int *ires, double *rpar, int *ipar); - + double *px; + double *pxdot; + double *pabs_tol; + double *prel_tol; + int *pinfo; + int *piwork; + double *prwork; }; #endif
--- a/liboctave/LSODE-opts.in Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/LSODE-opts.in Fri Aug 16 08:54:31 2002 +0000 @@ -15,10 +15,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION
--- a/liboctave/LSODE.cc Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/LSODE.cc Fri Aug 16 08:54:31 2002 +0000 @@ -54,28 +54,7 @@ static ODEFunc::ODEJacFunc user_jac; static ColumnVector *tmp_x; -LSODE::LSODE (void) : ODE (), LSODE_options () -{ - n = size (); - - itask = 1; - iopt = 0; - - sanity_checked = false; -} - -LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f) - : ODE (state, time, f), LSODE_options () -{ - n = size (); - - itask = 1; - iopt = 0; - - sanity_checked = false; -} - -int +static int lsode_f (const int& neq, const double& time, double *, double *deriv, int& ierr) { @@ -98,7 +77,7 @@ return 0; } -int +static int lsode_j (const int& neq, const double& time, double *, const int&, const int&, double *pd, const int& nrowpd) { @@ -122,60 +101,77 @@ { ColumnVector retval; - if (restart) + static int nn = 0; + + if (! initialized || restart || ODEFunc::reset || LSODE_options::reset) { - restart = false; - istate = 1; - } + integration_error = false; - if (integration_method () == "stiff") - { - if (jac) - method_flag = 21; - else - method_flag = 22; + initialized = true; + + istate = 1; + + int n = size (); + + nn = n; - liw = 20 + n; - lrw = 22 + n * (9 + n); - } - else - { - method_flag = 10; + if (integration_method () == "stiff") + { + if (jac) + method_flag = 21; + else + method_flag = 22; - liw = 20; - lrw = 22 + 16 * n; - } + liw = 20 + n; + lrw = 22 + n * (9 + n); + } + else + { + method_flag = 10; - if (iwork.length () != liw) - { + liw = 20; + lrw = 22 + 16 * n; + } + iwork.resize (liw); for (int i = 4; i < 9; i++) - iwork.elem (i) = 0; - } + iwork(i) = 0; - if (rwork.length () != lrw) - { rwork.resize (lrw); for (int i = 4; i < 9; i++) - rwork.elem (i) = 0; - } + rwork(i) = 0; - integration_error = false; + if (stop_time_set) + { + itask = 4; + rwork(0) = stop_time; + iopt = 1; + } + else + { + itask = 1; + } - double *xp = x.fortran_vec (); + px = x.fortran_vec (); - // NOTE: this won't work if LSODE passes copies of the state vector. - // In that case we have to create a temporary vector object - // and copy. + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + + restart = false; + + // ODEFunc - tmp_x = &x; - user_fun = function (); - user_jac = jacobian_function (); + // NOTE: this won't work if LSODE passes copies of the state vector. + // In that case we have to create a temporary vector object + // and copy. - if (! sanity_checked) - { + tmp_x = &x; + + user_fun = function (); + user_jac = jacobian_function (); + ColumnVector xdot = (*user_fun) (x, t); if (x.length () != xdot.length ()) @@ -187,70 +183,62 @@ return retval; } - sanity_checked = true; - } + ODEFunc::reset = false; + + // LSODE_options + + rel_tol = relative_tolerance (); + abs_tol = absolute_tolerance (); + + int abs_tol_len = abs_tol.length (); + + if (abs_tol_len == 1) + itol = 1; + else if (abs_tol_len == n) + itol = 2; + else + { + (*current_liboctave_error_handler) + ("lsode: inconsistent sizes for state and absolute tolerance vectors"); + + integration_error = true; + return retval; + } - if (stop_time_set) - { - itask = 4; - rwork.elem (0) = stop_time; - iopt = 1; - } - else - { - itask = 1; + double iss = initial_step_size (); + if (iss >= 0.0) + { + rwork(4) = iss; + iopt = 1; + } + + double maxss = maximum_step_size (); + if (maxss >= 0.0) + { + rwork(5) = maxss; + iopt = 1; + } + + double minss = minimum_step_size (); + if (minss >= 0.0) + { + rwork(6) = minss; + iopt = 1; + } + + int sl = step_limit (); + if (sl > 0) + { + iwork(5) = sl; + iopt = 1; + } + + pabs_tol = abs_tol.fortran_vec (); + + LSODE_options::reset = false; } - double rel_tol = relative_tolerance (); - - const Array<double> abs_tol = absolute_tolerance (); - - int abs_tol_len = abs_tol.length (); - - int itol; - - if (abs_tol_len == 1) - itol = 1; - else if (abs_tol_len == n) - itol = 2; - else - { - (*current_liboctave_error_handler) - ("lsode: inconsistent sizes for state and absolute tolerance vectors"); - - integration_error = true; - return retval; - } - - if (initial_step_size () >= 0.0) - { - rwork.elem (4) = initial_step_size (); - iopt = 1; - } - - if (maximum_step_size () >= 0.0) - { - rwork.elem (5) = maximum_step_size (); - iopt = 1; - } - - if (minimum_step_size () >= 0.0) - { - rwork.elem (6) = minimum_step_size (); - iopt = 1; - } - - if (step_limit () > 0) - { - iwork.elem (5) = step_limit (); - iopt = 1; - } - - const double *pabs_tol = abs_tol.fortran_vec (); - int *piwork = iwork.fortran_vec (); - double *prwork = rwork.fortran_vec (); - - F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol, + F77_XFCN (lsode, LSODE, (lsode_f, nn, px, t, tout, itol, rel_tol, pabs_tol, itask, istate, iopt, prwork, lrw, piwork, liw, lsode_j, method_flag)); @@ -365,7 +353,9 @@ LSODE::do_integrate (const ColumnVector& tout) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) { @@ -393,7 +383,9 @@ LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) { Matrix retval; + int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) {
--- a/liboctave/LSODE.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/LSODE.h Fri Aug 16 08:54:31 2002 +0000 @@ -37,11 +37,10 @@ { public: - LSODE (void); + LSODE (void) : ODE (), LSODE_options (), initialized (false) { } - LSODE (int n); - - LSODE (const ColumnVector& state, double time, const ODEFunc& f); + LSODE (const ColumnVector& state, double time, const ODEFunc& f) + : ODE (state, time, f), LSODE_options (), initialized (false) { } ~LSODE (void) { } @@ -55,20 +54,27 @@ private: - int n; + bool initialized; + int method_flag; - Array<int> iwork; - Array<double> rwork; int itask; int iopt; + int itol; + int liw; int lrw; - bool sanity_checked; + + Array<int> iwork; + Array<double> rwork; + + double rel_tol; - friend int lsode_f (int *neq, double *t, double *y, double *ydot); + Array<double> abs_tol; - friend int lsode_j (int *neq, double *t, double *y, int *ml, int *mu, - double *pd, int *nrowpd); + double *px; + double *pabs_tol; + int *piwork; + double *prwork; }; #endif
--- a/liboctave/ODEFunc.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/ODEFunc.h Fri Aug 16 08:54:31 2002 +0000 @@ -35,16 +35,16 @@ typedef Matrix (*ODEJacFunc) (const ColumnVector&, double); ODEFunc (void) - : fun (0), jac (0) { } + : fun (0), jac (0), reset (true) { } ODEFunc (ODERHSFunc f) - : fun (f), jac (0) { } + : fun (f), jac (0), reset (true) { } ODEFunc (ODERHSFunc f, ODEJacFunc j) - : fun (f), jac (j) { } + : fun (f), jac (j), reset (true) { } ODEFunc (const ODEFunc& a) - : fun (a.fun), jac (a.jac) { } + : fun (a.fun), jac (a.jac), reset (true) { } ODEFunc& operator = (const ODEFunc& a) { @@ -52,6 +52,7 @@ { fun = a.fun; jac = a.jac; + reset = a.reset; } return *this; } @@ -63,6 +64,7 @@ ODEFunc& set_function (ODERHSFunc f) { fun = f; + reset = true; return *this; } @@ -71,6 +73,7 @@ ODEFunc& set_jacobian_function (ODEJacFunc j) { jac = j; + reset = true; return *this; } @@ -78,6 +81,13 @@ ODERHSFunc fun; ODEJacFunc jac; + + // This variable is TRUE when this object is constructed, and also + // after any internal data has changed. Derived classes may use + // this information (and change it) to know when to (re)initialize + // their own internal data related to this object. + + bool reset; }; #endif
--- a/liboctave/ODESSA-opts.in Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/ODESSA-opts.in Fri Aug 16 08:54:31 2002 +0000 @@ -15,10 +15,11 @@ { $OPTVAR.resize (1); $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + reset = true; } void set_$OPT (const $TYPE& val) - { $OPTVAR = val; } + { $OPTVAR = val; reset = true; } END_SET_CODE END_OPTION
--- a/liboctave/ODESSA.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/ODESSA.h Fri Aug 16 08:54:31 2002 +0000 @@ -109,14 +109,13 @@ int n; int npar; - // Hey, check out this crap: ZZZZ + // XXX FIXME XXX -- ??? Array<double> par; Matrix sx0; Matrix y; - double *py; double *ppar; int *piwork;
--- a/liboctave/base-de.h Fri Aug 16 02:43:14 2002 +0000 +++ b/liboctave/base-de.h Fri Aug 16 08:54:31 2002 +0000 @@ -82,9 +82,14 @@ { stop_time_set = true; stop_time = t; + force_restart (); } - void clear_stop_time (void) { stop_time_set = false; } + void clear_stop_time (void) + { + stop_time_set = false; + force_restart (); + } virtual void force_restart (void) { restart = true; }
--- a/mk-opts.pl Fri Aug 16 02:43:14 2002 +0000 +++ b/mk-opts.pl Fri Aug 16 08:54:31 2002 +0000 @@ -374,7 +374,8 @@ } } - print " }\n"; + print " reset = true; + }\n"; print "\n void copy (const ${class_name}& opt)\n {\n"; @@ -383,7 +384,8 @@ print " $optvar[$i] = opt.$optvar[$i];\n"; } - print " }\n"; + print " reset = opt.reset; + }\n"; print "\n void set_default_options (void) { init (); }\n"; @@ -393,7 +395,7 @@ { &emit_set_decl ($i); - print "\n { $optvar[$i] = $set_expr[$i]; }\n"; + print "\n { $optvar[$i] = $set_expr[$i]; reset = true; }\n"; } elsif ($set_body[$i]) { @@ -403,7 +405,7 @@ chop ($s); $s =~ s/^/ /g; $s =~ s/\n/\n /g; - print "\n {\n$s\n }\n"; + print "\n {\n$s\n reset = true;\n }\n"; } elsif ($set_code[$i]) { @@ -427,7 +429,7 @@ print " $type[$i] $optvar[$i];\n"; } - print "};\n\n#endif\n"; + print "\nprotected:\n\n bool reset;\n};\n\n#endif\n"; } sub emit_set_decl