Mercurial > hg > octave-terminal
changeset 289:c23f50e61c58
[project @ 1994-01-13 06:25:58 by jwe]
author | jwe |
---|---|
date | Thu, 13 Jan 1994 06:26:54 +0000 |
parents | f8ae4f4dc9fd |
children | 5325afdb0ef5 |
files | liboctave/DASSL.cc liboctave/LSODE.cc liboctave/NLEqn.cc liboctave/NLEqn.h liboctave/ODE.h liboctave/QPSOL.cc liboctave/QPSOL.h liboctave/Quad.cc liboctave/Quad.h src/dassl.cc src/fsolve.cc src/g-builtins.cc src/lsode.cc src/qpsol.cc src/quad.cc |
diffstat | 15 files changed, 1173 insertions(+), 115 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -75,9 +75,6 @@ n = size; t = 0.0; - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -105,9 +102,6 @@ x = state; xdot.resize (n, 0.0); - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -144,9 +138,6 @@ xdot = deriv; x = state; - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -279,6 +270,25 @@ else info [3] = 0; + double abs_tol = absolute_tolerance (); + double rel_tol = relative_tolerance (); + + if (initial_step_size () >= 0.0) + { + rwork[2] = initial_step_size (); + info[7] = 1; + } + else + info[7] = 0; + + if (maximum_step_size () >= 0.0) + { + rwork[2] = maximum_step_size (); + info[6] = 1; + } + else + info[6] = 0; + double dummy; int idummy; @@ -291,9 +301,8 @@ again: F77_FCN (ddassl) (ddassl_f, &n, &t, px, pxdot, &tout, info, - &relative_tolerance, &absolute_tolerance, &idid, - rwork, &lrw, iwork, &liw, &dummy, &idummy, - ddassl_j); + &rel_tol, &abs_tol, &idid, rwork, &lrw, iwork, + &liw, &dummy, &idummy, ddassl_j); switch (idid) { @@ -475,3 +484,10 @@ return retval; } + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
--- a/liboctave/LSODE.cc +++ b/liboctave/LSODE.cc @@ -25,6 +25,9 @@ #include "config.h" #endif +#include <math.h> +#include <float.h> + #include "ODE.h" #include "f77-uscore.h" #include "lo-error.h" @@ -48,9 +51,6 @@ n = 0; t = 0.0; - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -82,9 +82,6 @@ n = size; t = 0.0; - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -117,9 +114,6 @@ t = time; x = state; - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; - stop_time_set = 0; stop_time = 0.0; @@ -129,7 +123,7 @@ istate = 1; itol = 1; itask = 1; - iopt = 0; + iopt = 1; liw = 20 + n; lrw = 22 + n * (9 + n); @@ -224,16 +218,21 @@ if (stop_time_set) { - iopt = 1; itask = 4; rwork [0] = stop_time; } else { - iopt = 0; itask = 1; } + double abs_tol = absolute_tolerance (); + double rel_tol = relative_tolerance (); + + rwork[4] = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0; + rwork[5] = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0; + rwork[6] = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0; + if (restart) { restart = 0; @@ -243,9 +242,9 @@ again: (void) F77_FCN (lsode) (lsode_f, &n, xp, &t, &tout, &itol, - &relative_tolerance, &absolute_tolerance, - &itask, &istate, &iopt, rwork, &lrw, iwork, - &liw, lsode_j, &method_flag); + &rel_tol, &abs_tol, &itask, &istate, &iopt, + rwork, &lrw, iwork, &liw, lsode_j, + &method_flag); switch (istate) { @@ -480,6 +479,116 @@ stop_time_set = 0; } +ODE_options::ODE_options (void) +{ + init (); +} + +ODE_options::ODE_options (const ODE_options& opt) +{ + copy (opt); +} + +ODE_options& +ODE_options::operator = (const ODE_options& opt) +{ + if (this != &opt) + copy (opt); + + return *this; +} + +ODE_options::~ODE_options (void) +{ +} + +void +ODE_options::init (void) +{ + double sqrt_eps = sqrt (DBL_EPSILON); + x_absolute_tolerance = sqrt_eps; + x_initial_step_size = -1.0; + x_maximum_step_size = -1.0; + x_minimum_step_size = 0.0; + x_relative_tolerance = sqrt_eps; +} + +void +ODE_options::copy (const ODE_options& opt) +{ + x_absolute_tolerance = opt.x_absolute_tolerance; + x_initial_step_size = opt.x_initial_step_size; + x_maximum_step_size = opt.x_maximum_step_size; + x_minimum_step_size = opt.x_minimum_step_size; + x_relative_tolerance = opt.x_relative_tolerance; +} + +void +ODE_options::set_default_options (void) +{ + init (); +} + +void +ODE_options::set_absolute_tolerance (double val) +{ + x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +void +ODE_options::set_initial_step_size (double val) +{ + x_initial_step_size = (val >= 0.0) ? val : -1.0; +} + +void +ODE_options::set_maximum_step_size (double val) +{ + x_maximum_step_size = (val >= 0.0) ? val : -1.0; +} + +void +ODE_options::set_minimum_step_size (double val) +{ + x_minimum_step_size = (val >= 0.0) ? val : 0.0; +} + +void +ODE_options::set_relative_tolerance (double val) +{ + x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +double +ODE_options::absolute_tolerance (void) +{ + return x_absolute_tolerance; +} + +double +ODE_options::initial_step_size (void) +{ + return x_initial_step_size; +} + +double +ODE_options::maximum_step_size (void) +{ + return x_maximum_step_size; +} + +double +ODE_options::minimum_step_size (void) +{ + return x_minimum_step_size; +} + +double +ODE_options::relative_tolerance (void) +{ + return x_relative_tolerance; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/NLEqn.cc +++ b/liboctave/NLEqn.cc @@ -211,7 +211,7 @@ return Vector (); } - double tol = sqrt (DBL_EPSILON); + double tol = tolerance (); double *fvec = new double [n]; double *px = new double [n]; @@ -258,6 +258,60 @@ return retval; } +NLEqn_options::NLEqn_options (void) +{ + init (); +} + +NLEqn_options::NLEqn_options (const NLEqn_options& opt) +{ + copy (opt); +} + +NLEqn_options& +NLEqn_options::operator = (const NLEqn_options& opt) +{ + if (this != &opt) + copy (opt); + + return *this; +} + +NLEqn_options::~NLEqn_options (void) +{ +} + +void +NLEqn_options::init (void) +{ + double sqrt_eps = sqrt (DBL_EPSILON); + x_tolerance = sqrt_eps; +} + +void +NLEqn_options::copy (const NLEqn_options& opt) +{ + x_tolerance = opt.x_tolerance; +} + +void +NLEqn_options::set_default_options (void) +{ + init (); +} + +void +NLEqn_options::set_tolerance (double val) +{ + x_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +double +NLEqn_options::tolerance (void) +{ + return x_tolerance; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/NLEqn.h +++ b/liboctave/NLEqn.h @@ -31,7 +31,32 @@ #define Vector ColumnVector #endif -class NLEqn : public NLFunc +class NLEqn_options +{ + public: + + NLEqn_options (void); + NLEqn_options (const NLEqn_options& opt); + + NLEqn_options& operator = (const NLEqn_options& opt); + + ~NLEqn_options (void); + + void init (void); + void copy (const NLEqn_options& opt); + + void set_default_options (void); + + void set_tolerance (double); + + double tolerance (void); + + private: + + double x_tolerance; +}; + +class NLEqn : public NLFunc, public NLEqn_options { public:
--- a/liboctave/ODE.h +++ b/liboctave/ODE.h @@ -29,7 +29,44 @@ #include "Matrix.h" #include "ODEFunc.h" -class ODE : public ODEFunc +class ODE_options +{ + public: + + ODE_options (void); + ODE_options (const ODE_options& opt); + + ODE_options& operator = (const ODE_options& opt); + + ~ODE_options (void); + + void init (void); + void copy (const ODE_options& opt); + + void set_default_options (void); + + void set_absolute_tolerance (double); + void set_initial_step_size (double); + void set_maximum_step_size (double); + void set_minimum_step_size (double); + void set_relative_tolerance (double); + + double absolute_tolerance (void); + double initial_step_size (void); + double maximum_step_size (void); + double minimum_step_size (void); + double relative_tolerance (void); + + private: + + double x_absolute_tolerance; + double x_initial_step_size; + double x_maximum_step_size; + double x_minimum_step_size; + double x_relative_tolerance; +}; + +class ODE : public ODEFunc, public ODE_options { public: @@ -68,9 +105,6 @@ double t; ColumnVector x; - double absolute_tolerance; - double relative_tolerance; - double stop_time; int stop_time_set;
--- a/liboctave/QPSOL.cc +++ b/liboctave/QPSOL.cc @@ -26,6 +26,7 @@ #endif #include <math.h> +#include <float.h> #ifndef QPSOL_MISSING @@ -36,9 +37,11 @@ { int F77_FCN (qpsol) (int*, int*, int*, int*, int*, int*, int*, int*, double*, double*, double*, double*, double*, - double*, double*, int (*)(), int*, int*, int*, - int*, double*, int*, int*, double*, double*, - int*, int*, double*, int*); + double*, double*, + int (*)(int*, int*, int*, int*, double*, + double*, double*), + int*, int*, int*, int*, double*, int*, int*, + double*, double*, int*, int*, double*, int*); int F77_FCN (dgemv) (const char*, const int*, const int*, const double*, const double*, const int*, @@ -80,13 +83,13 @@ int i; int n = x.capacity (); - - int itmax = 50 * n; - int msglvl = 0; + + int itmax = (iteration_limit () < 0) ? 50 * n : iteration_limit (); + int msglvl = print_level (); int nclin = lc.size (); int nctotl = nclin + n; - double bigbnd = 1e30; + double bigbnd = infinite_bound (); double dummy; double *pa = &dummy; @@ -125,10 +128,10 @@ double *pc = c.fortran_vec (); - double sqrt_eps = sqrt (DBL_EPSILON); double *featol = new double [nctotl]; + double tmp = feasibility_tolerance (); for (i = 0; i < nctotl; i++) - featol[i] = sqrt_eps; + featol[i] = tmp; double *ph = H.fortran_vec (); @@ -169,10 +172,99 @@ return x; } +QPSOL_options::QPSOL_options (void) +{ + init (); +} + +QPSOL_options::QPSOL_options (const QPSOL_options& opt) +{ + copy (opt); +} + +QPSOL_options& +QPSOL_options::operator = (const QPSOL_options& opt) +{ + if (this != &opt) + copy (opt); + + return *this; +} + +QPSOL_options::~QPSOL_options (void) +{ +} + void -QPSOL::set_default_options (void) +QPSOL_options::init (void) +{ + x_feasibility_tolerance = sqrt (DBL_EPSILON); + x_infinite_bound = 1.0e+30; + x_iteration_limit = -1; + x_print_level = 0; +} + +void +QPSOL_options::copy (const QPSOL_options& opt) +{ + x_feasibility_tolerance = opt.x_feasibility_tolerance; + x_infinite_bound = opt.x_infinite_bound; + x_iteration_limit = opt.x_iteration_limit; + x_print_level = opt.x_print_level; +} + +void +QPSOL_options::set_default_options (void) +{ + init (); +} + +void +QPSOL_options::set_feasibility_tolerance (double val) +{ + x_feasibility_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +void +QPSOL_options::set_infinite_bound (double val) { - iprint = 0; + x_infinite_bound = (val > 0.0) ? val : 1.0e+30; +} + +void +QPSOL_options::set_iteration_limit (int val) +{ + x_iteration_limit = (val > 0) ? val : -1; +} + +void +QPSOL_options::set_print_level (int val) +{ + x_print_level = (val >= 0) ? val : 0; +} + +double +QPSOL_options::feasibility_tolerance (void) +{ + return x_feasibility_tolerance; +} + +double +QPSOL_options::infinite_bound (void) +{ + return x_infinite_bound; +} + +int +QPSOL_options::iteration_limit (void) +{ + return x_iteration_limit; +} + +int +QPSOL_options::print_level (void) +{ + return x_print_level; } #endif /* QPSOL_MISSING */
--- a/liboctave/QPSOL.h +++ b/liboctave/QPSOL.h @@ -33,51 +33,81 @@ #define Vector ColumnVector #endif -class QPSOL : public QP +class QPSOL_options +{ + public: + + QPSOL_options (void); + QPSOL_options (const QPSOL_options& opt); + + QPSOL_options& operator = (const QPSOL_options& opt); + + ~QPSOL_options (void); + + void init (void); + void copy (const QPSOL_options& opt); + + void set_default_options (void); + + void set_feasibility_tolerance (double); + void set_infinite_bound (double); + void set_iteration_limit (int); + void set_print_level (int); + + double feasibility_tolerance (void); + double infinite_bound (void); + int iteration_limit (void); + int print_level (void); + + private: + + double x_feasibility_tolerance; + double x_infinite_bound; + int x_iteration_limit; + int x_print_level; +}; + +class QPSOL : public QP, public QPSOL_options { public: QPSOL (void) : QP () - { set_default_options (); } + { } QPSOL (const Vector& x, const Matrix& H) : QP (x, H) - { set_default_options (); } + { } QPSOL (const Vector& x, const Matrix& H, const Vector& c) : QP (x, H, c) - { set_default_options (); } + { } QPSOL (const Vector& x, const Matrix& H, const Bounds& b) : QP (x, H, b) - { set_default_options (); } + { } QPSOL (const Vector& x, const Matrix& H, const LinConst& lc) : QP (x, H, lc) - { set_default_options (); } + { } QPSOL (const Vector& x, const Matrix& H, const Vector& c, const Bounds& b) - : QP (x, H, c, b) { set_default_options (); } + : QP (x, H, c, b) { } QPSOL (const Vector& x, const Matrix& H, const Vector& c, const LinConst& lc) - : QP (x, H, c, lc) { set_default_options (); } + : QP (x, H, c, lc) { } QPSOL (const Vector& x, const Matrix& H, const Bounds& b, const LinConst& lc) - : QP (x, H, b, lc) { set_default_options (); } + : QP (x, H, b, lc) { } QPSOL (const Vector& x, const Matrix& H, const Vector& c, const Bounds& b, const LinConst& lc) - : QP (x, H, c, b, lc) { set_default_options (); } + : QP (x, H, c, b, lc) { } QPSOL (const QPSOL& a); QPSOL& operator = (const QPSOL& a); Vector minimize (double& objf, int& inform, Vector& lambda); - -private: - void set_default_options (void); - int iprint; }; inline QPSOL::QPSOL (const QPSOL& a) : QP (a.x, a.H, a.c, a.bnds, a.lc) - { set_default_options (); } + { } inline QPSOL& QPSOL::operator = (const QPSOL& a) @@ -87,7 +117,6 @@ c = a.c; bnds = a.bnds; lc = a.lc; - iprint = a.iprint; return *this; }
--- a/liboctave/Quad.cc +++ b/liboctave/Quad.cc @@ -25,6 +25,9 @@ #include "config.h" #endif +#include <math.h> +#include <float.h> + #include "Quad.h" #include "f77-uscore.h" #include "sun-utils.h" @@ -53,15 +56,11 @@ Quad::Quad (integrand_fcn fcn) { - absolute_tolerance = 1.0e-6; - relative_tolerance = 1.0e-6; f = fcn; } Quad::Quad (integrand_fcn fcn, double abs, double rel) { - absolute_tolerance = abs; - relative_tolerance = rel; f = fcn; } @@ -173,10 +172,12 @@ user_fcn = f; int last; + double abs_tol = absolute_tolerance (); + double rel_tol = relative_tolerance (); + F77_FCN (dqagp) (user_function, &lower_limit, &upper_limit, &npts, - points, &absolute_tolerance, &relative_tolerance, - &result, &abserr, &neval, &ier, &leniw, &lenw, - &last, iwork, work); + points, &abs_tol, &rel_tol, &result, &abserr, + &neval, &ier, &leniw, &lenw, &last, iwork, work); delete [] iwork; delete [] work; @@ -239,9 +240,12 @@ break; } - F77_FCN (dqagi) (user_function, &bound, &inf, &absolute_tolerance, - &relative_tolerance, &result, &abserr, &neval, - &ier, &leniw, &lenw, &last, iwork, work); + double abs_tol = absolute_tolerance (); + double rel_tol = relative_tolerance (); + + F77_FCN (dqagi) (user_function, &bound, &inf, &abs_tol, &rel_tol, + &result, &abserr, &neval, &ier, &leniw, &lenw, + &last, iwork, work); delete [] iwork; delete [] work; @@ -249,6 +253,74 @@ return result; } +Quad_options::Quad_options (void) +{ + init (); +} + +Quad_options::Quad_options (const Quad_options& opt) +{ + copy (opt); +} + +Quad_options& +Quad_options::operator = (const Quad_options& opt) +{ + if (this != &opt) + copy (opt); + + return *this; +} + +Quad_options::~Quad_options (void) +{ +} + +void +Quad_options::init (void) +{ + double sqrt_eps = sqrt (DBL_EPSILON); + x_absolute_tolerance = sqrt_eps; + x_relative_tolerance = sqrt_eps; +} + +void +Quad_options::copy (const Quad_options& opt) +{ + x_absolute_tolerance = opt.x_absolute_tolerance; + x_relative_tolerance = opt.x_relative_tolerance; +} + +void +Quad_options::set_default_options (void) +{ + init (); +} + +void +Quad_options::set_absolute_tolerance (double val) +{ + x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +void +Quad_options::set_relative_tolerance (double val) +{ + x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +double +Quad_options::absolute_tolerance (void) +{ + return x_absolute_tolerance; +} + +double +Quad_options::relative_tolerance (void) +{ + return x_relative_tolerance; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/Quad.h +++ b/liboctave/Quad.h @@ -43,7 +43,35 @@ // function, and the user wants us to quit. extern int quad_integration_error; -class Quad +class Quad_options +{ + public: + + Quad_options (void); + Quad_options (const Quad_options& opt); + + Quad_options& operator = (const Quad_options& opt); + + ~Quad_options (void); + + void init (void); + void copy (const Quad_options& opt); + + void set_default_options (void); + + void set_absolute_tolerance (double); + void set_relative_tolerance (double); + + double absolute_tolerance (void); + double relative_tolerance (void); + + private: + + double x_absolute_tolerance; + double x_relative_tolerance; +}; + +class Quad : public Quad_options { public: @@ -57,9 +85,6 @@ protected: - double absolute_tolerance; - double relative_tolerance; - integrand_fcn f; };
--- a/src/dassl.cc +++ b/src/dassl.cc @@ -25,6 +25,8 @@ #include "config.h" #endif +#include <strstream.h> + #include "DAE.h" #include "tree-const.h" @@ -32,6 +34,7 @@ #include "gripes.h" #include "error.h" #include "utils.h" +#include "pager.h" #include "f-dassl.h" // Global pointer for user defined function required by dassl. @@ -51,6 +54,8 @@ } #endif +static ODE_options dassl_opts; + ColumnVector dassl_user_function (const ColumnVector& x, const ColumnVector& xdot, double t) { @@ -153,6 +158,7 @@ DAEFunc func (dassl_user_function); DAE dae (state, deriv, tzero, func); + dae.copy (dassl_opts); Matrix output; Matrix deriv_output; @@ -168,13 +174,131 @@ return retval; } +typedef void (ODE_options::*d_set_opt_mf) (double); +typedef double (ODE_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 3 + +struct ODE_OPTIONS +{ + char *keyword; + char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static ODE_OPTIONS dassl_option_table[] = +{ + { "absolute tolerance", + { "absolute", "tolerance", NULL, NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_absolute_tolerance, + ODE_options::absolute_tolerance, }, + + { "initial step size", + { "initial", "step", "size", NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_initial_step_size, + ODE_options::initial_step_size, }, + + { "maximum step size", + { "maximum", "step", "size", NULL, }, + { 2, 0, 0, 0, }, 1, + ODE_options::set_maximum_step_size, + ODE_options::maximum_step_size, }, + + { "relative tolerance", + { "relative", "tolerance", NULL, NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_relative_tolerance, + ODE_options::relative_tolerance, }, + + { NULL, + { NULL, NULL, NULL, NULL, }, + { 0, 0, 0, 0, }, 0, + NULL, NULL, }, +}; + +static void +print_dassl_option_list (void) +{ + ostrstream output_buf; + + print_usage ("dassl_options", 1); + + output_buf << "\n" + << "Options for dassl include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + ODE_OPTIONS *list = dassl_option_table; + + char *keyword; + while ((keyword = list->keyword) != (char *) NULL) + { + output_buf.form (" %-40s ", keyword); + + double val = (dassl_opts.*list->d_get_fcn) (); + if (val < 0.0) + output_buf << "computed automatically"; + else + output_buf << val; + + output_buf << "\n"; + list++; + } + + output_buf << "\n" << ends; + maybe_page_output (output_buf); +} + +static void +do_dassl_option (char *keyword, double val) +{ + ODE_OPTIONS *list = dassl_option_table; + + while (list->keyword != (char *) NULL) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (dassl_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("dassl_options: no match for `%s'", keyword); +} + tree_constant * dassl_options (const tree_constant *args, int nargin, int nargout) { -// Assumes that we have been given the correct number of arguments. + tree_constant *retval = NULL_TREE_CONST; - tree_constant *retval = NULL_TREE_CONST; - error ("dassl_options: not implemented yet"); + if (nargin == 1) + { + print_dassl_option_list (); + } + else if (nargin == 3) + { + if (args[1].is_string_type ()) + { + char *keyword = args[1].string_value (); + double val = args[2].double_value (); + do_dassl_option (keyword, val); + } + else + print_usage ("dassl_options"); + } + else + { + print_usage ("dassl_options"); + } + return retval; }
--- a/src/fsolve.cc +++ b/src/fsolve.cc @@ -25,6 +25,8 @@ #include "config.h" #endif +#include <strstream.h> + #include "NLEqn.h" #include "tree-const.h" @@ -32,6 +34,7 @@ #include "gripes.h" #include "error.h" #include "utils.h" +#include "pager.h" #include "f-fsolve.h" // Global pointer for user defined function required by hybrd1. @@ -51,6 +54,8 @@ } #endif +static NLEqn_options fsolve_opts; + int hybrd_info_to_fsolve_info (int info) { @@ -143,13 +148,14 @@ ColumnVector x = args[2].to_vector (); if (nargin > 3) - warning ("fsolve: ignoring optional arguments"); + warning ("fsolve: ignoring extra arguments"); if (nargout > 2) warning ("fsolve: can't compute path output yet"); NLFunc foo_fcn (fsolve_user_function); NLEqn foo (x, foo_fcn); + foo.copy (fsolve_opts); int info; ColumnVector soln = foo.solve (info); @@ -168,13 +174,113 @@ return retval; } +typedef void (NLEqn_options::*d_set_opt_mf) (double); +typedef double (NLEqn_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 1 + +struct NLEQN_OPTIONS +{ + char *keyword; + char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static NLEQN_OPTIONS fsolve_option_table[] = +{ + { "tolerance", + { "tolerance", NULL, }, + { 1, 0, }, 1, + NLEqn_options::set_tolerance, + NLEqn_options::tolerance, }, + + { NULL, + { NULL, NULL, }, + { 0, 0, }, 0, + NULL, NULL, }, +}; + +static void +print_fsolve_option_list (void) +{ + ostrstream output_buf; + + print_usage ("fsolve_options", 1); + + output_buf << "\n" + << "Options for fsolve include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + NLEQN_OPTIONS *list = fsolve_option_table; + + char *keyword; + while ((keyword = list->keyword) != (char *) NULL) + { + output_buf.form (" %-40s ", keyword); + + double val = (fsolve_opts.*list->d_get_fcn) (); + if (val < 0.0) + output_buf << "computed automatically"; + else + output_buf << val; + + output_buf << "\n"; + list++; + } + + output_buf << "\n" << ends; + maybe_page_output (output_buf); +} + +static void +do_fsolve_option (char *keyword, double val) +{ + NLEQN_OPTIONS *list = fsolve_option_table; + + while (list->keyword != (char *) NULL) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (fsolve_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("fsolve_options: no match for `%s'", keyword); +} + tree_constant * fsolve_options (const tree_constant *args, int nargin, int nargout) { -// Assumes that we have been given the correct number of arguments. + tree_constant *retval = NULL_TREE_CONST; - tree_constant *retval = NULL_TREE_CONST; - error ("fsolve_options: not implemented yet"); + if (nargin == 1) + { + print_fsolve_option_list (); + } + else if (nargin == 3) + { + if (args[1].is_string_type ()) + { + char *keyword = args[1].string_value (); + double val = args[2].double_value (); + do_fsolve_option (keyword, val); + } + else + print_usage ("fsolve_options"); + } + else + { + print_usage ("fsolve_options"); + } + return retval; }
--- a/src/g-builtins.cc +++ b/src/g-builtins.cc @@ -335,11 +335,8 @@ { tree_constant *retval = NULL_TREE_CONST; - if (nargin >= 3 && nargin <= 7 && nargout >= 1 && nargout <= 3) - DLD_BUILTIN (args, nargin, nargout, dassl_options, - retval = dassl_options (args, nargin, nargout);) - else - print_usage ("dassl_options"); + DLD_BUILTIN (args, nargin, nargout, dassl_options, + retval = dassl_options (args, nargin, nargout);) return retval; } @@ -797,11 +794,8 @@ { tree_constant *retval = NULL_TREE_CONST; - if (nargin >= 3 && nargin <= 7 && nargout >= 1 && nargout <= 3) - DLD_BUILTIN (args, nargin, nargout, fsolve_options, - retval = fsolve_options (args, nargin, nargout);) - else - print_usage ("fsolve_options"); + DLD_BUILTIN (args, nargin, nargout, fsolve_options, + retval = fsolve_options (args, nargin, nargout);) return retval; } @@ -1049,11 +1043,8 @@ { tree_constant *retval = NULL_TREE_CONST; - if (nargin >= 3 && nargin <= 7 && nargout >= 1 && nargout <= 3) - DLD_BUILTIN (args, nargin, nargout, lpsolve_options, - retval = lpsolve_options (args, nargin, nargout);) - else - print_usage ("lpsolve_options"); + DLD_BUILTIN (args, nargin, nargout, lpsolve_options, + retval = lpsolve_options (args, nargin, nargout);) return retval; } @@ -1080,11 +1071,8 @@ { tree_constant *retval = NULL_TREE_CONST; - if (nargin >= 3 && nargin <= 7 && nargout >= 1 && nargout <= 3) - DLD_BUILTIN (args, nargin, nargout, lsode_options, - retval = lsode_options (args, nargin, nargout);) - else - print_usage ("lsode_options"); + DLD_BUILTIN (args, nargin, nargout, lsode_options, + retval = lsode_options (args, nargin, nargout);) return retval; } @@ -1417,11 +1405,8 @@ { tree_constant *retval = NULL_TREE_CONST; - if (nargin >= 3 && nargin <= 7 && nargout >= 1 && nargout <= 3) - DLD_BUILTIN (args, nargin, nargout, quad_options, - retval = quad_options (args, nargin, nargout);) - else - print_usage ("quad_options"); + DLD_BUILTIN (args, nargin, nargout, quad_options, + retval = quad_options (args, nargin, nargout);) return retval; }
--- a/src/lsode.cc +++ b/src/lsode.cc @@ -25,6 +25,8 @@ #include "config.h" #endif +#include <strstream.h> + #include "ODE.h" #include "tree-const.h" @@ -32,6 +34,7 @@ #include "gripes.h" #include "error.h" #include "utils.h" +#include "pager.h" #include "f-lsode.h" // Global pointer for user defined function required by lsode. @@ -51,6 +54,8 @@ } #endif +static ODE_options lsode_opts; + ColumnVector lsode_user_function (const ColumnVector& x, double t) { @@ -136,6 +141,7 @@ ODEFunc func (lsode_user_function); ODE ode (state, tzero, func); + ode.copy (lsode_opts); int nstates = state.capacity (); Matrix output (nsteps, nstates + 1); @@ -150,13 +156,137 @@ return retval; } +typedef void (ODE_options::*d_set_opt_mf) (double); +typedef double (ODE_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 3 + +struct ODE_OPTIONS +{ + char *keyword; + char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static ODE_OPTIONS lsode_option_table[] = +{ + { "absolute tolerance", + { "absolute", "tolerance", NULL, NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_absolute_tolerance, + ODE_options::absolute_tolerance, }, + + { "initial step size", + { "initial", "step", "size", NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_initial_step_size, + ODE_options::initial_step_size, }, + + { "maximum step size", + { "maximum", "step", "size", NULL, }, + { 2, 0, 0, 0, }, 1, + ODE_options::set_maximum_step_size, + ODE_options::maximum_step_size, }, + + { "minimum step size", + { "minimum", "step", "size", NULL, }, + { 2, 0, 0, 0, }, 1, + ODE_options::set_minimum_step_size, + ODE_options::minimum_step_size, }, + + { "relative tolerance", + { "relative", "tolerance", NULL, NULL, }, + { 1, 0, 0, 0, }, 1, + ODE_options::set_relative_tolerance, + ODE_options::relative_tolerance, }, + + { NULL, + { NULL, NULL, NULL, NULL, }, + { 0, 0, 0, 0, }, 0, + NULL, NULL, }, +}; + +static void +print_lsode_option_list (void) +{ + ostrstream output_buf; + + print_usage ("lsode_options", 1); + + output_buf << "\n" + << "Options for lsode include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + ODE_OPTIONS *list = lsode_option_table; + + char *keyword; + while ((keyword = list->keyword) != (char *) NULL) + { + output_buf.form (" %-40s ", keyword); + + double val = (lsode_opts.*list->d_get_fcn) (); + if (val < 0.0) + output_buf << "computed automatically"; + else + output_buf << val; + + output_buf << "\n"; + list++; + } + + output_buf << "\n" << ends; + maybe_page_output (output_buf); +} + +static void +do_lsode_option (char *keyword, double val) +{ + ODE_OPTIONS *list = lsode_option_table; + + while (list->keyword != (char *) NULL) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (lsode_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("lsode_options: no match for `%s'", keyword); +} + tree_constant * lsode_options (const tree_constant *args, int nargin, int nargout) { -// Assumes that we have been given the correct number of arguments. + tree_constant *retval = NULL_TREE_CONST; - tree_constant *retval = NULL_TREE_CONST; - error ("lsode_options: not implemented yet"); + if (nargin == 1) + { + print_lsode_option_list (); + } + else if (nargin == 3) + { + if (args[1].is_string_type ()) + { + char *keyword = args[1].string_value (); + double val = args[2].double_value (); + do_lsode_option (keyword, val); + } + else + print_usage ("lsode_options"); + } + else + { + print_usage ("lsode_options"); + } + return retval; }
--- a/src/qpsol.cc +++ b/src/qpsol.cc @@ -27,6 +27,8 @@ #ifndef QPSOL_MISSING +#include <strstream.h> + #include "QPSOL.h" #include "tree-const.h" @@ -34,6 +36,7 @@ #include "gripes.h" #include "error.h" #include "utils.h" +#include "pager.h" #include "f-qpsol.h" // This should probably be defined in some shared file and declared in @@ -57,6 +60,8 @@ } #endif +static QPSOL_options qpsol_opts; + tree_constant * qpsol (const tree_constant *args, int nargin, int nargout) { @@ -126,6 +131,7 @@ // 1. qpsol (x, H, c) QPSOL qp (x, H, c); + qp.copy (qpsol_opts); soln = qp.minimize (objf, inform, lambda); goto solved; @@ -136,6 +142,7 @@ // 2. qpsol (x, H, c, lb, ub) QPSOL qp (x, H, c, bounds); + qp.copy (qpsol_opts); soln = qp.minimize (objf, inform, lambda); goto solved; @@ -163,6 +170,7 @@ // 3. qpsol (x, H, c, lb, ub, llb, A, lub) QPSOL qp (x, H, c, bounds, linear_constraints); + qp.copy (qpsol_opts); soln = qp.minimize (objf, inform, lambda); } else @@ -170,6 +178,7 @@ // 4. qpsol (x, H, c, llb, A, lub) QPSOL qp (x, H, c, linear_constraints); + qp.copy (qpsol_opts); soln = qp.minimize (objf, inform, lambda); } goto solved; @@ -191,13 +200,147 @@ return retval; } +typedef void (QPSOL_options::*d_set_opt_mf) (double); +typedef void (QPSOL_options::*i_set_opt_mf) (int); +typedef double (QPSOL_options::*d_get_opt_mf) (void); +typedef int (QPSOL_options::*i_get_opt_mf) (void); + +#define MAX_TOKENS 2 + +struct QPSOL_OPTIONS +{ + char *keyword; + char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; + d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; +}; + +static QPSOL_OPTIONS qpsol_option_table[] = +{ + { "feasibility tolerance", + { "feasibility", "tolerance", NULL, }, + { 1, 0, 0, }, 1, + QPSOL_options::set_feasibility_tolerance, NULL, + QPSOL_options::feasibility_tolerance, NULL, }, + + { "infinite bound", + { "infinite", "bound", NULL, }, + { 2, 0, 0, }, 1, + QPSOL_options::set_infinite_bound, NULL, + QPSOL_options::infinite_bound, NULL, }, + + { "iteration limit", + { "iteration", "limit", NULL, }, + { 2, 0, 0, }, 1, + NULL, QPSOL_options::set_iteration_limit, + NULL, QPSOL_options::iteration_limit, }, + + { "print level", + { "print", "level", NULL, }, + { 1, 0, 0, }, 1, + NULL, QPSOL_options::set_print_level, + NULL, QPSOL_options::print_level, }, + + { NULL, + { NULL, NULL, NULL, }, + { 0, 0, 0, }, 0, + NULL, NULL, NULL, NULL, }, +}; + +static void +print_qpsol_option_list (void) +{ + ostrstream output_buf; + + print_usage ("qpsol_options", 1); + + output_buf << "\n" + << "Options for qpsol include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + QPSOL_OPTIONS *list = qpsol_option_table; + + char *keyword; + while ((keyword = list->keyword) != (char *) NULL) + { + output_buf.form (" %-40s ", keyword); + if (list->d_get_fcn) + { + double val = (qpsol_opts.*list->d_get_fcn) (); + if (val < 0.0) + output_buf << "computed automatically"; + else + output_buf << val; + } + else + { + int val = (qpsol_opts.*list->i_get_fcn) (); + if (val < 0) + output_buf << "depends on problem size"; + else + output_buf << val; + } + output_buf << "\n"; + list++; + } + + output_buf << "\n" << ends; + maybe_page_output (output_buf); +} + +static void +do_qpsol_option (char *keyword, double val) +{ + QPSOL_OPTIONS *list = qpsol_option_table; + + while (list->keyword != (char *) NULL) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_set_fcn) + (qpsol_opts.*list->d_set_fcn) (val); + else + (qpsol_opts.*list->i_set_fcn) (NINT (val)); + + return; + } + list++; + } + + warning ("qpsol_options: no match for `%s'", keyword); +} + tree_constant * qpsol_options (const tree_constant *args, int nargin, int nargout) { -// Assumes that we have been given the correct number of arguments. + tree_constant *retval = NULL_TREE_CONST; - tree_constant *retval = NULL_TREE_CONST; - error ("qpsol_options: not implemented yet"); + if (nargin == 1) + { + print_qpsol_option_list (); + } + else if (nargin == 3) + { + if (args[1].is_string_type ()) + { + char *keyword = args[1].string_value (); + double val = args[2].double_value (); + do_qpsol_option (keyword, val); + } + else + print_usage ("qpsol_options"); + } + else + { + print_usage ("qpsol_options"); + } + return retval; }
--- a/src/quad.cc +++ b/src/quad.cc @@ -25,6 +25,8 @@ #include "config.h" #endif +#include <strstream.h> + #include "Quad.h" #include "tree-const.h" @@ -33,6 +35,7 @@ #include "gripes.h" #include "error.h" #include "utils.h" +#include "pager.h" #include "f-quad.h" // Global pointer for user defined function required by quadrature functions. @@ -52,6 +55,8 @@ } #endif +static Quad_options quad_opts; + double quad_user_function (double x) { @@ -164,6 +169,7 @@ if (indefinite) { IndefQuad iq (quad_user_function, bound, indef_type, abstol, reltol); + iq.copy (quad_opts); val = iq.integrate (ier, nfun, abserr); } else @@ -171,11 +177,13 @@ if (have_sing) { DefQuad dq (quad_user_function, a, b, sing, abstol, reltol); + dq.copy (quad_opts); val = dq.integrate (ier, nfun, abserr); } else { DefQuad dq (quad_user_function, a, b, abstol, reltol); + dq.copy (quad_opts); val = dq.integrate (ier, nfun, abserr); } } @@ -195,13 +203,119 @@ return retval; } +typedef void (Quad_options::*d_set_opt_mf) (double); +typedef double (Quad_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 2 + +struct QUAD_OPTIONS +{ + char *keyword; + char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static QUAD_OPTIONS quad_option_table[] = +{ + { "absolute tolerance", + { "absolute", "tolerance", NULL, }, + { 1, 0, 0, }, 1, + Quad_options::set_absolute_tolerance, + Quad_options::absolute_tolerance, }, + + { "relative tolerance", + { "relative", "tolerance", NULL, }, + { 1, 0, 0, }, 1, + Quad_options::set_relative_tolerance, + Quad_options::relative_tolerance, }, + + { NULL, + { NULL, NULL, NULL, }, + { 0, 0, 0, }, 0, + NULL, NULL, }, +}; + +static void +print_quad_option_list (void) +{ + ostrstream output_buf; + + print_usage ("quad_options", 1); + + output_buf << "\n" + << "Options for quad include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + QUAD_OPTIONS *list = quad_option_table; + + char *keyword; + while ((keyword = list->keyword) != (char *) NULL) + { + output_buf.form (" %-40s ", keyword); + + double val = (quad_opts.*list->d_get_fcn) (); + if (val < 0.0) + output_buf << "computed automatically"; + else + output_buf << val; + + output_buf << "\n"; + list++; + } + + output_buf << "\n" << ends; + maybe_page_output (output_buf); +} + +static void +do_quad_option (char *keyword, double val) +{ + QUAD_OPTIONS *list = quad_option_table; + + while (list->keyword != (char *) NULL) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (quad_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("quad_options: no match for `%s'", keyword); +} + tree_constant * quad_options (const tree_constant *args, int nargin, int nargout) { -// Assumes that we have been given the correct number of arguments. + tree_constant *retval = NULL_TREE_CONST; - tree_constant *retval = NULL_TREE_CONST; - error ("quad_options: not implemented yet"); + if (nargin == 1) + { + print_quad_option_list (); + } + else if (nargin == 3) + { + if (args[1].is_string_type ()) + { + char *keyword = args[1].string_value (); + double val = args[2].double_value (); + do_quad_option (keyword, val); + } + else + print_usage ("quad_options"); + } + else + { + print_usage ("quad_options"); + } + return retval; }