Mercurial > hg > octave-lyh
changeset 3984:addebffd4961
[project @ 2002-07-11 03:39:33 by jwe]
author | jwe |
---|---|
date | Thu, 11 Jul 2002 03:39:34 +0000 |
parents | 7a37caf6ed43 |
children | fa0ae9105656 |
files | liboctave/ChangeLog liboctave/Makefile.in liboctave/ODE.h liboctave/ODES.cc liboctave/ODES.h liboctave/ODESFunc.h liboctave/ODESSA.cc liboctave/ODESSA.h liboctave/base-de.h src/ChangeLog src/DLD-FUNCTIONS/odessa.cc src/Makefile.in |
diffstat | 12 files changed, 2138 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,13 @@ +2002-07-10 John W. Eaton <jwe@bevo.che.wisc.edu> + + * ODE.h: Move integrate and do_integrate method declarations and + definitions here. + * base-de.h: From here. + + * ODES.h, ODES.cc, ODESFunc.h, ODESSA.h, ODESSA.cc: New files. + * Makefile.in: Add them to the appropriate lists. + (LIBOCTAVE_CXX_SOURCES): + 2002-07-02 John W. Eaton <jwe@bevo.che.wisc.edu> * NLEqn.cc (NLEqn::error_message): New function.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -46,7 +46,8 @@ INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DASPK.h DASSL.h FEGrid.h \ LinConst.h LP.h LPsolve.h LSODE.h NLConst.h NLEqn.h NLFunc.h \ - NLP.h ODE.h ODEFunc.h Objective.h QP.h Quad.h Range.h base-de.h \ + NLP.h ODE.h ODEFunc.h ODES.h ODESFunc.h ODESSA.h \ + Objective.h QP.h Quad.h Range.h base-de.h \ base-min.h byte-swap.h cmd-edit.h cmd-hist.h data-conv.h \ dir-ops.h file-ops.h file-stat.h getopt.h glob-match.h \ idx-vector.h lo-ieee.h lo-mappers.h lo-specfun.h lo-sysdep.h \ @@ -84,9 +85,10 @@ vx-rv-cs.cc vx-s-ccv.cc vx-s-crv.cc \ vx-rv-crv.cc vx-cv-ccv.cc vx-crv-rv.cc vx-ccv-cv.cc -LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc DASPK.cc \ - DASSL.cc FEGrid.cc LinConst.cc LPsolve.cc LSODE.cc \ - NLEqn.cc Quad.cc Range.cc data-conv.cc dir-ops.cc \ +LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc \ + DASPK.cc DASSL.cc FEGrid.cc LinConst.cc LPsolve.cc \ + LSODE.cc NLEqn.cc ODES.cc ODESSA.cc Quad.cc \ + Range.cc data-conv.cc dir-ops.cc \ file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \ lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \ lo-utils.cc mach-info.cc oct-alloc.cc oct-env.cc \
--- a/liboctave/ODE.h +++ b/liboctave/ODE.h @@ -51,6 +51,64 @@ } ~ODE (void) { } + + // Derived classes must provide functions to actually do the + // integration. + + // Return the vector of states at output time t. + virtual ColumnVector do_integrate (double tt) = 0; + + // Return a matrix of states at each output time specified by t. + // The rows of the result matrix should each correspond to a new + // output time. + virtual Matrix do_integrate (const ColumnVector& tt) = 0; + + virtual Matrix do_integrate (const ColumnVector& tt, + const ColumnVector& ttcrit) = 0; + + // Lots of ways to call the single function and optionally set and + // get additional information. + + // Integrate to t from current point. + virtual ColumnVector integrate (double tt) + { return do_integrate (tt); } + + // Set new x0, t0 and integrate to t. + virtual ColumnVector integrate (const ColumnVector& x0, double t0, double tt) + { + initialize (x0, t0); + return do_integrate (tt); + } + + // Integrate from current point and return output at all points + // specified by t. + virtual Matrix integrate (const ColumnVector& tt) + { return do_integrate (tt); } + + // Set new x0, t0 and integrate to return output at all points + // specified by t. + virtual Matrix integrate (const ColumnVector& x0, double t0, + const ColumnVector& tt) + { + initialize (x0, t0); + return do_integrate (tt); + } + + // Integrate from current point and return output at all points + // specified by t. + virtual Matrix integrate (const ColumnVector& tt, + const ColumnVector& ttcrit) + { return do_integrate (tt, ttcrit); } + + // Set new x0, t0 and integrate to return output at all points + // specified by t. + virtual Matrix integrate (const ColumnVector& x0, double t0, + const ColumnVector& tt, + const ColumnVector& ttcrit) + { + initialize (x0, t0); + return do_integrate (tt, ttcrit); + } }; #endif
new file mode 100644 --- /dev/null +++ b/liboctave/ODES.cc @@ -0,0 +1,54 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "ODES.h" +#include "lo-error.h" + +void +ODES::initialize (const ColumnVector& xx, double t) +{ + base_diff_eqn::initialize (xx, t); + xdot = ColumnVector (xx.length (), 0.0); +} + +void +ODES::initialize (const ColumnVector& xx, double t, + const ColumnVector& xtheta) +{ + base_diff_eqn::initialize (xx, t); + xdot = ColumnVector (xx.length (), 0.0); + theta = xtheta; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/ODES.h @@ -0,0 +1,89 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if !defined (octave_ODES_h) +#define octave_ODES_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "ODESFunc.h" +#include "base-de.h" + +class +ODES : public base_diff_eqn, public ODESFunc +{ +public: + + ODES (void) + : base_diff_eqn (), ODESFunc (), theta () { } + + ODES (const ColumnVector& x, double t, ODESFunc& f) + : base_diff_eqn (x, t), ODESFunc (f), xdot (x.length (), 0.0), theta () { } + + ODES (const ColumnVector& x, const ColumnVector& xtheta, double t, + ODESFunc& f) + : base_diff_eqn (x, t), ODESFunc (f), xdot (x.length (), 0.0), + theta (xtheta) { } + + ODES (const ODES& a) + : base_diff_eqn (a), ODESFunc (a), theta (a.theta) { } + + ODES& operator = (const ODES& a) + { + if (this != &a) + { + base_diff_eqn::operator = (a); + ODESFunc::operator = (a); + + xdot = a.xdot; + theta = a.theta; + } + return *this; + } + + ~ODES (void) { } + + ColumnVector parameter_vector (void) { return theta; } + + void initialize (const ColumnVector& x, double t); + + void initialize (const ColumnVector& x, double t, + const ColumnVector& theta); + +protected: + + // State vector time derivatives. + ColumnVector xdot; + + // Parameter vector. + ColumnVector theta; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/ODESFunc.h @@ -0,0 +1,117 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if !defined (octave_ODESFunc_h) +#define octave_ODESFunc_h 1 + +#include "dMatrix.h" + +class +ODESFunc +{ +public: + + struct DAEJac + { + Matrix *dfdxdot; + Matrix *dfdx; + }; + + typedef ColumnVector (*ODES_fsub) (double, const ColumnVector& x, + const ColumnVector& theta); + + typedef ColumnVector (*ODES_bsub) (double, const ColumnVector& x, + const ColumnVector& theta, int column); + + typedef Matrix (*ODES_jsub) (double, const ColumnVector& x, + 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) + { + if (this != &a) + { + fsub = a.fsub; + bsub = a.bsub; + jsub = a.jsub; + } + return *this; + } + + ~ODESFunc (void) { } + + ODES_fsub fsub_function (void) const { return fsub; } + + ODESFunc& set_fsub_function (ODES_fsub f) + { + fsub = f; + return *this; + } + + ODES_bsub bsub_function (void) const { return bsub; } + + ODESFunc& set_bsub_function (ODES_bsub b) + { + bsub = b; + return *this; + } + + ODES_jsub jsub_function (void) const { return jsub; } + + ODESFunc& set_jsub_function (ODES_jsub j) + { + jsub = j; + return *this; + } + +protected: + + ODES_fsub fsub; + ODES_bsub bsub; + ODES_jsub jsub; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + + + +
new file mode 100644 --- /dev/null +++ b/liboctave/ODESSA.cc @@ -0,0 +1,722 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream.h> + +#include <cstdlib> +#include <cfloat> +#include <cmath> + +// For instantiating the Array<Matrix> object. +#include "Array.h" +#include "Array.cc" + +#include "ODESSA.h" +#include "f77-fcn.h" +#include "lo-error.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "ov-fcn.h" +#include "pager.h" +#include "parse.h" +#include "unwind-prot.h" +#include "utils.h" +#include "variables.h" + +#ifndef F77_FUNC +#define F77_FUNC(x, X) F77_FCN (x, X) +#endif +void +ODESSA_options::set_integration_method (const std::string& val) +{ + if (val == "stiff" || val == "bdf") + x_integration_method = "stiff"; + else if (val == "non-stiff" || val == "adams") + x_integration_method = "non-stiff"; + else + (*current_liboctave_error_handler) + ("odessa_options: method must be \"stiff\", \"bdf\", \"non-stiff\", or \"adams\""); +} + +typedef int (*odessa_fcn_ptr) (int*, const double&, double*, + double*, double*); + +typedef int (*odessa_jac_ptr) (int*, const double&, double*, + double*, const int&, const int&, + double*, const int&); + +typedef int (*odessa_dfdp_ptr) (int*, const double&, double*, + double*, double*, const int&); + + +extern "C" +int F77_FUNC (odessa, ODESSA) (odessa_fcn_ptr, odessa_dfdp_ptr, int*, + double*, double*, double&, double&, + int&, double&, const double*, int&, + int&, int*, double*, int&, int*, int&, + odessa_jac_ptr, int&); + +template class Array<Matrix>; + +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, + double *par, double *fval) +{ + int n = neq[0]; + int n_par = neq[1]; + + // Load the state and parameter arrays as Octave objects + + ColumnVector tmp_state (n); + for (int i = 0; i < n; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_param (n_par); + for (int i = 0; i < n_par; i++) + tmp_param(i) = par[i]; + + ColumnVector tmp_fval = user_fsub (t, tmp_state, tmp_param); + + // Return the function value as a C array object + + for (int i = 0; i < n; i++) + { + fval[i] = tmp_fval(i); + } + + return 0; +} + +static int +odessa_j (int* neq, const double& t, double *state, + double *par, const int& ml, const int& mu, double *pd, const + int& nrowpd) +{ + int n = neq[0]; + int n_par = neq[1]; + + // Load the state and parameter arrays as Octave objects + ColumnVector tmp_state (n); + for (int i = 0; i < n; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_param (n_par); + for (int i = 0; i < n_par; i++) + tmp_param(i) = par[i]; + + Matrix tmp_fval = user_jsub (t, tmp_state, tmp_param); + + for (int j = 0; j < n; j++) + for (int i = 0; i < nrowpd; i++) + pd[nrowpd*j+i] = tmp_fval(i,j); + + return 0; +} + +static int +odessa_b (int* neq, const double& t, double *state, + double *par, double *dfdp, const int& jpar) + +{ + int n = neq[0]; + int n_par = neq[1]; + + // Load the state and parameter arrays as Octave objects + ColumnVector tmp_state (n); + for (int i = 0; i < n; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_param (n_par); + for (int i = 0; i < n_par; i++) + tmp_param(i) = par[i]; + + ColumnVector tmp_fval = user_bsub (t, tmp_state, tmp_param, jpar); + + for (int i = 0; i < n; i++) + dfdp[i] = tmp_fval(i); + + return 0; +} + +ODESSA::ODESSA (void) : ODES (), ODESSA_options () +{ + neq.resize(2); + n = size (); + + stop_time_set = false; + stop_time = 0.0; + + integration_error = false; + restart = true; + iopt.resize(4); + + istate = 1; + itask = 1; + iopt(0) = 0; + isopt = 0; + iopt(1) = isopt; + npar = 0; + neq(0) = n; + neq(1) = npar; + + sanity_checked = false; +} + +ODESSA::ODESSA (const ColumnVector& state, double time, ODESFunc& f) + : ODES (state, time, f), ODESSA_options () +{ + neq.resize(2); + n = size (); + + stop_time_set = false; + stop_time = 0.0; + + ODESFunc::operator = (f); + + integration_error = false; + restart = true; + + iopt.resize(4); + istate = 1; + itask = 1; + iopt(0) = 0; + isopt = 0; + iopt(1) = isopt; + + sanity_checked = false; + + npar = 0; + neq(0) = n; + neq(1) = npar; + + y.resize (n, 1, 0.0); +} + +ODESSA::ODESSA (const ColumnVector& state, const ColumnVector& theta, + const Matrix& sensitivity_guess, double time, ODESFunc& f) + : ODES (state, theta, time, f) +{ + initialized = false; + restart = false; + + neq.resize(2); + n = state.length(); + npar = theta.length(); + + neq(0) = n; + neq(1) = npar; + + stop_time_set = false; + stop_time = 0.0; + + ODESFunc::operator = (f); + + sx0 = sensitivity_guess; + par.resize (npar); + + for (int i = 0; i < npar; i++) + { + par(i) = theta(i); + } + + sanity_checked = false; + + npar = theta.length (); + + iopt.resize(4); + istate = 1; + itask = 1; + iopt(0) = 0; + isopt = 1; + iopt(1) = isopt; + + y.resize (n, npar+1, 0.0); +} + +void +ODESSA::force_restart (void) +{ + restart = true; +} + +void +ODESSA::set_stop_time (double time) +{ + stop_time_set = true; + stop_time = time; +} + +void +ODESSA::clear_stop_time (void) +{ + stop_time_set = 0; +} + +void +ODESSA::integrate (double tout) +{ + ODESSA_result retval; + if (! initialized) + { + + for (int i = 0; i < n; i++) + y(i,0) = x(i); + + if (npar > 0) + { + for (int j = 0; j < npar; j++) + for (int i = 0; i < n; i++) + y(i,j+1) = sx0(i,j); + } + + 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; + + iopt(2) = idf; + + + if (restart) + { + restart = false; + istate = 1; + } + + if (integration_method () == "stiff") + { + if (user_jsub) + method_flag = 21; + else + method_flag = 22; + if (isopt) + { + liw = 21 + n + npar; + lrw = 22 + 8*(npar+1)*n + n*n + n; + } + else + { + liw = 20 + n; + lrw = 22 + 9*n + n*n; + } + } + else + { + if (isopt) + { + if (user_jsub) + method_flag = 11; + else + method_flag = 12; + liw = 21 + n + npar; + lrw = 22 + 15*(npar+1)*n + n*n + n; + } + else + { + method_flag = 10; + liw = 20 + n; + lrw = 22 + 16 * n; + } + } + + if (iwork.length () != liw) + { + iwork.resize (liw); + + for (int i = 4; i < 10; i++) + iwork.elem (i) = 0; + } + + if (rwork.length () != lrw) + { + rwork.resize (lrw); + + for (int i = 4; i < 10; i++) + rwork.elem (i) = 0.0; + } + + initialized = true; + } + integration_error = false; + + // 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) + { + ColumnVector fval = user_fsub (t, x, theta); + + if (fval.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("odessa: inconsistent sizes for state and residual vectors"); + + integration_error = true; + return; + } + + sanity_checked = true; + } + + if (stop_time_set) + { + itask = 4; + rwork.elem (0) = stop_time; + iopt(0) = 1; + } + else + { + itask = 1; + } + + double rel_tol = relative_tolerance (); + + Array<double> abs_tol = absolute_tolerance (); //note; this should + // be a matrix, not a vector + + 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 = 1; + return; + } + + if (initial_step_size () >= 0.0) + { + rwork.elem (4) = initial_step_size (); + iopt(0) = 1; + } + + if (maximum_step_size () >= 0.0) + { + rwork.elem (5) = maximum_step_size (); + iopt(0) = 1; + } + + if (minimum_step_size () >= 0.0) + { + rwork.elem (6) = minimum_step_size (); + iopt(0) = 1; + } + + if (step_limit () > 0) + { + iwork.elem (5) = step_limit (); + iopt(0) = 1; + } + + py = y.fortran_vec (); + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + ppar = par.fortran_vec (); + piopt = iopt.fortran_vec (); + pneq = neq.fortran_vec (); + + const double *pabs_tol = abs_tol.fortran_vec (); + + F77_XFCN (odessa, ODESSA, (odessa_f,odessa_b,pneq,py,ppar,t, tout, + itol,rel_tol,pabs_tol,itask,istate, + piopt,prwork,lrw,piwork,liw,odessa_j, + method_flag)); + + if (f77_exception_encountered) + { + integration_error = true; + (*current_liboctave_error_handler) ("unrecoverable error in odessa"); + } + else + { + switch (istate) + { + case -13: // Return requested in user-supplied function. + case -6: // error weight became zero during problem. (solution + // component i vanished, and atol or atol(i) = 0.) + case -5: // repeated convergence failures (perhaps bad jacobian + // supplied or wrong choice of mf or tolerances). + case -4: // repeated error test failures (check all inputs). + case -3: // illegal input detected (see printed message). + case -2: // excess accuracy requested (tolerances too small). + case -1: // excess work done on this call (perhaps wrong mf). + integration_error = 1; + break; + + case 1: // tout is same as t + case 2: // lsode was successful + t = tout; + break; + + default: // Error? + (*current_liboctave_error_handler) + ("unrecognized value of istate returned from lsode"); + break; + } + } +} + +std::string +ODESSA::error_message (void) const +{ + std::string retval; + + switch (istate) + { + case -13: + retval = "return requested in user-supplied function"; + break; + + case -6: + retval = "Error weight became zero during problem.\ + (solution component i vanished, and atol or atol(i) == 0)"; + break; + + case -5: + retval = "repeated convergence failures (perhaps bad jacobian\ + supplied or wrong choice of integration method or tolerances)"; + break; + + case -4: + retval = "repeated error test failures (check all inputs)"; + break; + + case -3: + retval = "invalid input detected (see printed message)"; + break; + + case -2: + retval = "excess accuracy requested (tolerances too small)"; + break; + + case -1: + retval = "excess work on this call (perhaps wrong integration method)"; + break; + + case 1: + retval = "prior to initial call"; + break; + + case 2: + retval = "successful exit"; + break; + + case 3: + retval = "prior to continuation call with modified parameters"; + break; + + default: + retval = "unknown error state"; + break; + } + + return retval; +} + + +ODESSA_result +ODESSA::integrate (const ColumnVector& tout) +{ + ODESSA_result retval; + + Matrix x_out; + + Array<Matrix> x_s_out; + + int n_out = tout.capacity (); + + if (n_out > 0 && n > 0) + { + x_out.resize (n_out, n); + + x_s_out.resize (npar, Matrix (n_out, n, 0.0)); + + for (int j = 0; j < n_out; j++) + { + integrate (tout(j)); + + if (integration_error) + { + retval = ODESSA_result (x_out, x_s_out); + return retval; + } + + for (int i = 0; i < n; i++) + { + x_out(j,i) = y(i,0); + + for (int k = 0; k < npar; k++) + { + x_s_out(k)(j,i) = y(i,k+1); + } + } + } + } + + retval = ODESSA_result (x_out, x_s_out); + + return retval; +} + +ODESSA_result +ODESSA::integrate (const ColumnVector& tout, const ColumnVector& tcrit) +{ + ODESSA_result retval; + + Matrix x_out; + + Array<Matrix> x_s_out; + + int n_out = tout.capacity (); + + if (n_out > 0 && n > 0) + { + x_out.resize (n_out, n); + + x_s_out.resize (npar, Matrix (n_out, n, 0.0)); + + int n_crit = tcrit.capacity (); + + if (n_crit > 0) + { + int i_crit = 0; + int i_out = 0; + double next_crit = tcrit(0); + double next_out; + while (i_out < n_out) + { + bool do_restart = false; + + next_out = tout(i_out); + if (i_crit < n_crit) + next_crit = tcrit(i_crit); + + int save_output; + double t_out; + + if (next_crit == next_out) + { + set_stop_time (next_crit); + t_out = next_out; + save_output = true; + i_out++; + i_crit++; + do_restart = true; + } + else if (next_crit < next_out) + { + if (i_crit < n_crit) + { + set_stop_time (next_crit); + t_out = next_crit; + save_output = false; + i_crit++; + do_restart = true; + } + else + { + clear_stop_time (); + t_out = next_out; + save_output = true; + i_out++; + } + } + else + { + set_stop_time (next_crit); + t_out = next_out; + save_output = true; + i_out++; + } + integrate (t_out); + if (integration_error) + { + retval = ODESSA_result (x_out, x_s_out); + return retval; + } + if (save_output) + { + for (int i = 0; i < n; i++) + { + x_out(i_out-1,i) = y(i,0); + + for (int k = 0; k < npar; k++) + { + x_s_out(k)(i_out-1,i) = y(i,k+1); + } + } + } + + if (do_restart) + force_restart (); + } + + retval = ODESSA_result (x_out, x_s_out); + } + else + { + retval = integrate (tout); + + if (integration_error) + return retval; + } + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/ODESSA.h @@ -0,0 +1,265 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if !defined (octave_ODESSA_h) +#define octave_ODESSA_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include <cfloat> +#include <cmath> + +#include "ODES.h" + +class +ODESSA_options +{ +public: + + ODESSA_options (void) { init (); } + + ODESSA_options (const ODESSA_options& opt) { copy (opt); } + + ODESSA_options& operator = (const ODESSA_options& opt) + { + if (this != &opt) + copy (opt); + + return *this; + } + + ~ODESSA_options (void) { } + + void init (void) + { + double sqrt_eps = ::sqrt (DBL_EPSILON); + x_absolute_tolerance.resize (1); + x_absolute_tolerance(0) = sqrt_eps; + x_initial_step_size = -1.0; + x_integration_method = "stiff"; + x_maximum_step_size = -1.0; + x_minimum_step_size = 0.0; + x_relative_tolerance = sqrt_eps; + + // This is consistent with earlier versions of Octave, and is + // much larger than the default of 500 specified in the LSODE + // sources. + x_step_limit = 100000; + } + + void copy (const ODESSA_options& opt) + { + x_absolute_tolerance = opt.x_absolute_tolerance; + x_initial_step_size = opt.x_initial_step_size; + x_integration_method = opt.x_integration_method; + 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; + x_step_limit = opt.x_step_limit; + } + + void set_default_options (void) { init (); } + + void set_absolute_tolerance (double val) + { + x_absolute_tolerance.resize (1); + x_absolute_tolerance(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); + } + + void set_absolute_tolerance (const Array<double>& val) + { x_absolute_tolerance = val; } + + void set_initial_step_size (double val) + { x_initial_step_size = (val >= 0.0) ? val : -1.0; } + + void set_integration_method (const std::string& val); + + + void set_maximum_step_size (double val) + { x_maximum_step_size = (val >= 0.0) ? val : -1.0; } + + void set_minimum_step_size (double val) + { x_minimum_step_size = (val >= 0.0) ? val : 0.0; } + + void set_relative_tolerance (double val) + { x_relative_tolerance = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); } + + void set_step_limit (int val) + { x_step_limit = val; } + + Array<double> absolute_tolerance (void) const + { return x_absolute_tolerance; } + + double initial_step_size (void) const + { return x_initial_step_size; } + + std::string integration_method (void) const + { return x_integration_method; } + + double maximum_step_size (void) const + { return x_maximum_step_size; } + + double minimum_step_size (void) const + { return x_minimum_step_size; } + + double relative_tolerance (void) const + { return x_relative_tolerance; } + + int step_limit (void) const + { return x_step_limit; } + +private: + + Array<double> x_absolute_tolerance; + double x_initial_step_size; + std::string x_integration_method; + double x_maximum_step_size; + double x_minimum_step_size; + double x_relative_tolerance; + + int x_step_limit; +}; + + +class +ODESSA_result +{ +public: + + ODESSA_result (void) { } + + ODESSA_result (const Matrix& xx, + const Array<Matrix>& xx_s) + + : x (xx), x_s (xx_s) { } + + ODESSA_result (const ODESSA_result& r) + : x (r.x), x_s (r.x_s) { } + + ODESSA_result& operator = (const ODESSA_result& r) + { + if (this != &r) + { + x = r.x; + x_s = r.x_s; + } + return *this; + } + + ~ODESSA_result (void) { } + + Matrix state (void) const { return x; } + Array<Matrix> state_sensitivity (void) const { return x_s; } + +private: + + Matrix x; + Array<Matrix> x_s; +}; + +class +ODESSA : public ODES, public ODESSA_options +{ +public: + + ODESSA (void); + + ODESSA (const ColumnVector& x, double time, ODESFunc& f); + + ODESSA (const ColumnVector& x, const ColumnVector& theta, + const Matrix& sensitivity_guess, double time, ODESFunc& f); + + ~ODESSA (void) { } + + void force_restart (void); + + void set_stop_time (double t); + void clear_stop_time (void); + + ODESSA_result integrate (const ColumnVector& tout); + + ODESSA_result integrate (const ColumnVector& tout, + const ColumnVector& tcrit); + + int integration_state (void) const { return istate; } + + bool integration_ok (void) const { return ! integration_error; } + + std::string error_message (void) const; + +private: + + bool initialized; + + bool sanity_checked; + + bool stop_time_set; + double stop_time; + + bool restart; + + bool integration_error; + + int liw; + int lrw; + int method_flag; + Array<int> iwork; + Array<double> rwork; + int istate; + int itask; + Array<int> iopt; + int isopt; + + Array<int> neq; + + int n; + int npar; + + // Hey, check out this crap: ZZZZ + Array<double> par; + + Matrix sx0; + + Matrix y; + + + double *py; + double *ppar; + int *piwork; + int *piopt; + int *pneq; + double *prwork; + + void init_work_size (int); + + void integrate (double t); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/base-de.h +++ b/liboctave/base-de.h @@ -49,69 +49,11 @@ return *this; } - // Derived classes must provide functions to actually do the - // integration. - - // Return the vector of states at output time t. - virtual ColumnVector do_integrate (double tt) = 0; + // There must be a way for us to force the integration to restart. - // Return a matrix of states at each output time specified by t. - // The rows of the result matrix should each correspond to a new - // output time. - virtual Matrix do_integrate (const ColumnVector& tt) = 0; - - virtual Matrix do_integrate (const ColumnVector& tt, - const ColumnVector& ttcrit) = 0; - - // There must also be a way for us to force the integration to - // restart. virtual void force_restart (void) = 0; - // Lots of ways to call the single function and optionally set and - // get additional information. - - // Integrate to t from current point. - virtual ColumnVector integrate (double tt) - { return do_integrate (tt); } - - // Set new x0, t0 and integrate to t. - virtual ColumnVector integrate (const ColumnVector& x0, double t0, double tt) - { - initialize (x0, t0); - return do_integrate (tt); - } - - // Integrate from current point and return output at all points - // specified by t. - virtual Matrix integrate (const ColumnVector& tt) - { return do_integrate (tt); } - - // Set new x0, t0 and integrate to return output at all points - // specified by t. - virtual Matrix integrate (const ColumnVector& x0, double t0, - const ColumnVector& tt) - { - initialize (x0, t0); - return do_integrate (tt); - } - - // Integrate from current point and return output at all points - // specified by t. - virtual Matrix integrate (const ColumnVector& tt, - const ColumnVector& ttcrit) - { return do_integrate (tt, ttcrit); } - - // Set new x0, t0 and integrate to return output at all points - // specified by t. - virtual Matrix integrate (const ColumnVector& x0, double t0, - const ColumnVector& tt, - const ColumnVector& ttcrit) - { - initialize (x0, t0); - return do_integrate (tt, ttcrit); - } - - virtual void initialize (const ColumnVector& x0, double t0) + void initialize (const ColumnVector& x0, double t0) { x = x0; t = t0;
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2002-07-10 John W. Eaton <jwe@bevo.che.wisc.edu> + + * DLD-FUNCTIONS/odessa.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. + 2002-07-05 John W. Eaton <jwe@bevo.che.wisc.edu> * pt-assign.cc (tree_multi_assignment::rvalue): Call
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/odessa.cc @@ -0,0 +1,803 @@ +/* + +Copyright (C) 2002 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iomanip> +#include <iostream> + +#include "ODESSA.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "ov-fcn.h" +#include "pager.h" +#include "pr-output.h" +#include "unwind-prot.h" +#include "utils.h" +#include "variables.h" +#include "parse.h" + +// Global pointer for user defined function required by odessa. +static octave_function *odessa_f; +static octave_function *odessa_j; +static octave_function *odessa_b; + +static ODESSA_options odessa_opts; + +// Is this a recursive call? +static int call_depth = 0; + +static ColumnVector +odessa_user_f (double t, const ColumnVector& x, const ColumnVector& theta) +{ + ColumnVector retval; + + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + if (n > 1) + args(1) = x; + else if (n == 1) + args(1) = x(0); + else + args(1) = Matrix (); + + args(0) = t; + + if (odessa_f) + { + octave_value_list tmp = odessa_f->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = ColumnVector (tmp(0).vector_value ()); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + } + else + gripe_user_supplied_eval ("odessa"); + } + + return retval; +} + +static Matrix +odessa_user_mf (double t, const ColumnVector& x, const ColumnVector& theta, + octave_function *mf) +{ + Matrix retval; + + if (mf) + { + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + if (n > 1) + args(1) = x; + else if (n == 1) + args(1) = x(0); + else + args(1) = Matrix (); + + args(0) = t; + + octave_value_list tmp = mf->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).matrix_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + } + else + gripe_user_supplied_eval ("odessa"); + } + + return retval; +} + +static Matrix +odessa_user_j (double t, const ColumnVector& x, const ColumnVector& theta) +{ + return odessa_user_mf (t, x, theta, odessa_j); +} + +static ColumnVector +odessa_user_b (double t, const ColumnVector& x, + const ColumnVector& theta, int column) +{ + ColumnVector retval; + + if (odessa_b) + { + octave_value_list args; + + int n = x.length (); + int npar = theta.length (); + + args(3) = static_cast<double> (column); + + if (npar > 1) + args(2) = theta; + else if (npar == 1) + args(2) = theta(0); + else + args(2) = Matrix (); + + if (n > 1) + args(1) = x; + else if (n == 1) + args(1) = x(0); + else + args(1) = Matrix (); + + args(0) = t; + + octave_value_list tmp = odessa_b->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("odessa"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = ColumnVector (tmp(0).vector_value ()); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("odessa"); + } + else + gripe_user_supplied_eval ("odessa"); + } + + return retval; +} + +static octave_value_list +make_list (const Array<Matrix>& m_array) +{ + octave_value_list retval; + + int len = m_array.length (); + + retval.resize (len); + + for (int i = 0; i < len; i++) + retval(i) = m_array(i); + + return retval; +} + +#define ODESSA_ABORT() \ + do \ + { \ + unwind_protect::run_frame ("Fodessa"); \ + return retval; \ + } \ + while (0) + +#define ODESSA_ABORT1(msg) \ + do \ + { \ + ::error ("odessa: " msg); \ + ODESSA_ABORT (); \ + } \ + while (0) + +#define ODESSA_ABORT2(fmt, arg) \ + do \ + { \ + ::error ("odessa: " fmt, arg); \ + ODESSA_ABORT (); \ + } \ + while (0) + +// -------------------------------- +// Everthing is so great above here +// -------------------------------- + +DEFUN_DLD (odessa, args, , + "odessa (\"f\", x_0, theta, sx_0, t_out, t_crit)\n\ +\n\ +The string \"f\" may be substituted for the vector of strings\n\ +\n\ + [\"f\"; \"j\"; \"b\"] \n\ +\n\ +You can use the function @code{odessa_options} to set optional\n\ +parameters for @code{odessa}.") +{ + octave_value_list retval; + + unwind_protect::begin_frame ("Fodessa"); + + unwind_protect_int (call_depth); + call_depth++; + + if (call_depth > 1) + ODESSA_ABORT1 ("invalid recursive call"); + + int nargin = args.length (); + + if (nargin < 5 || nargin > 6) + { + print_usage ("odessa"); + unwind_protect::run_frame ("Fodessa"); + return retval; + } + + odessa_f = 0; + odessa_j = 0; + odessa_b = 0; + + octave_value f_arg = args(0); + + if (f_arg.is_string ()) + { + string_vector f_str_arg = f_arg.all_strings (); + + int len = f_str_arg.length (); + + if (len > 2) + { + string t = f_str_arg(2); + + if (t.length () > 0) + { + odessa_b = is_valid_function (t, "odessa", true); + + if (! odessa_b) + ODESSA_ABORT1 + ("expecting b function name as argument 1"); + } + } + + if (len > 1) + { + string t = f_str_arg(1); + + if (t.length () > 0) + { + odessa_j = is_valid_function (t, "odessa", true); + + if (! odessa_j) + ODESSA_ABORT1 + ("expecting function name as argument 1"); + } + } + + if (len > 0) + odessa_f = is_valid_function (f_str_arg(0), "odessa", true); + + if (! odessa_f) + ODESSA_ABORT1 ("expecting function name as argument 1"); + } + + ColumnVector state (args(1).vector_value ()); + + if (error_state) + ODESSA_ABORT1 ("expecting state vector as argument 2"); + + bool have_parameters = false; + + ColumnVector theta; + Matrix sensitivity_guess; + + if (nargin == 5 || nargin == 6) + { + octave_value theta_arg = args(2); + + if (! theta_arg.is_empty ()) + { + theta = ColumnVector (theta_arg.vector_value ()); + + if (error_state) + ODESSA_ABORT1 + ("expecting parameter vector as argument 3"); + } + + have_parameters = (theta.length () > 0); + + if (have_parameters) + { + sensitivity_guess = args(3).matrix_value (); + + if (error_state) + ODESSA_ABORT1 + ("expecting sensitivity guess as argument 4"); + + if (sensitivity_guess.rows () != state.length () + || sensitivity_guess.columns () != theta.length ()) + ODESSA_ABORT1 + ("incorrect dimension for sensitivity guess"); + } + } + + ColumnVector out_times (args(4).vector_value ()); + + if (error_state) + ODESSA_ABORT1 + ("expecting output time vector as %s argument 5"); + + ColumnVector crit_times; + + bool crit_times_set = false; + + if (nargin==6) + { + crit_times = ColumnVector (args(5).vector_value ()); + + if (error_state) + ODESSA_ABORT1 + ("expecting critical time vector as argument 6"); + + 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); + + ODESSA_result output; + + if (have_parameters) + { + ODESSA dae = ODESSA (state, theta, sensitivity_guess, tzero, func); + + dae.copy (odessa_opts); + + if (crit_times_set) + output = dae.integrate (out_times, crit_times); + else + output = dae.integrate (out_times); + } + else + { + ODESSA dae = ODESSA (state, tzero, func); + + dae.copy (odessa_opts); + + if (crit_times_set) + output = dae.integrate (out_times, crit_times); + else + output = dae.integrate (out_times); + } + + if (! error_state) + { + if (have_parameters) + { + retval(1) = make_list (output.state_sensitivity ()); + } + + retval(0) = output.state (); + } + + unwind_protect::run_frame ("Fodessa"); + + return retval; +} +// ----------------------------- +// EVERYTHING SWELL BELOW HERE +// ----------------------------- + + +typedef void (ODESSA_options::*da_set_opt_mf) (const Array<double>&); +typedef void (ODESSA_options::*d_set_opt_mf) (double); +typedef void (ODESSA_options::*i_set_opt_mf) (int); +typedef void (ODESSA_options::*s_set_opt_mf) (const std::string&); + +typedef Array<double> (ODESSA_options::*da_get_opt_mf) (void) const; +typedef double (ODESSA_options::*d_get_opt_mf) (void) const; +typedef int (ODESSA_options::*i_get_opt_mf) (void) const; +typedef std::string (ODESSA_options::*s_get_opt_mf) (void) const; + +#define MAX_TOKENS 3 + +struct ODESSA_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + da_set_opt_mf da_set_fcn; + d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; + s_set_opt_mf s_set_fcn; + da_get_opt_mf da_get_fcn; + d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; + s_get_opt_mf s_get_fcn; +}; + +static ODESSA_OPTIONS odessa_option_table [] = +{ + { "absolute tolerance", + { "absolute", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + &ODESSA_options::set_absolute_tolerance, 0, 0, 0, + &ODESSA_options::absolute_tolerance, 0, 0, 0, }, + + { "initial step size", + { "initial", "step", "size", 0, }, + { 3, 0, 0, 0, }, 1, + 0, &ODESSA_options::set_initial_step_size, 0, 0, + 0, &ODESSA_options::initial_step_size, 0, 0, }, + + { "integration method", + { "integration", "method", 0, 0, }, + { 3, 0, 0, 0, }, 1, + 0, 0, 0, &ODESSA_options::set_integration_method, + 0, 0, 0, &ODESSA_options::integration_method, }, + + { "maximum step size", + { "maximum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + 0, &ODESSA_options::set_maximum_step_size, 0, 0, + 0, &ODESSA_options::maximum_step_size, 0, 0, }, + + { "minimum step size", + { "minimum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + 0, &ODESSA_options::set_minimum_step_size, 0, 0, + 0, &ODESSA_options::minimum_step_size, 0, 0, }, + + { "relative tolerance", + { "relative", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + 0, &ODESSA_options::set_relative_tolerance, 0, 0, + 0, &ODESSA_options::relative_tolerance, 0, 0, }, + + { "step limit", + { "step", "limit", 0, 0, }, + { 1, 0, 0, 0, }, 1, + 0, 0, &ODESSA_options::set_step_limit, 0, + 0, 0, &ODESSA_options::step_limit, 0, }, + + { 0, + { 0, 0, 0, 0, }, + { 0, 0, 0, 0, }, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, }, +}; + +static void +print_odessa_option_list (std::ostream& os) +{ + print_usage ("odessa_options", 1); + + os << "\n" + << "Options for odessa include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + ODESSA_OPTIONS *list = odessa_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os << " " + << std::setiosflags (std::ios::left) << std::setw (40) + << keyword + << std::resetiosflags (std::ios::left) + << " "; + + if (list->da_get_fcn) + { + Array<double> val = (odessa_opts.*list->da_get_fcn) (); + if (val.length () == 1) + { + if (val(0) < 0.0) + os << "computed automatically"; + else + os << val(0); + } + else + { + os << "\n\n"; + // XXX FIXME XXX + Matrix tmp = Matrix (ColumnVector (val)); + octave_print_internal (os, tmp, false, 2); + os << "\n"; + } + } + else if (list->d_get_fcn) + { + double val = (odessa_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } + else if (list->i_get_fcn) + { + int val = (odessa_opts.*list->i_get_fcn) (); + if (val < 0) + os << "infinite"; + else + os << val; + } + else if (list->s_get_fcn) + { + os << (odessa_opts.*list->s_get_fcn) (); + } + else + panic_impossible (); + + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_odessa_option (const std::string& keyword, double val) +{ + ODESSA_OPTIONS *list = odessa_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->da_set_fcn) + { + Array<double> tmp (1, val); + (odessa_opts.*list->da_set_fcn) (tmp); + } + else if (list->d_set_fcn) + { + (odessa_opts.*list->d_set_fcn) (val); + } + else + { + if (xisnan (val)) + { + error ("odessa_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (odessa_opts.*list->i_set_fcn) (NINT (val)); + } + return; + } + list++; + } + + warning ("odessa_options: no match for `%s'", keyword.c_str ()); +} + +static void +set_odessa_option (const std::string& keyword, const Array<double>& val) +{ + ODESSA_OPTIONS *list = odessa_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->da_set_fcn) + (odessa_opts.*list->da_set_fcn) (val); + else + error ("odessa_options: no function to handle vector option"); + + return; + } + list++; + } + + warning ("odessa_options: no match for `%s'", keyword.c_str ()); +} + +static void +set_odessa_option (const std::string& keyword, const std::string& val) +{ + ODESSA_OPTIONS *list = odessa_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->s_set_fcn) + (odessa_opts.*list->s_set_fcn) (val); + else + error ("odessa_options: no function to handle string option"); + + return; + } + list++; + } + + warning ("odessa_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_odessa_option (const std::string& keyword) +{ + octave_value retval; + + ODESSA_OPTIONS *list = odessa_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->da_get_fcn) + { + Array<double> val = (odessa_opts.*list->da_get_fcn) (); + if (val.length () == 1) + { + if (val(0) < 0.0) + retval = "computed automatically"; + else + retval = val(0); + } + else + retval = ColumnVector (val); + } + else if (list->d_get_fcn) + { + double val = (odessa_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } + else if (list->i_get_fcn) + { + int val = (odessa_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "infinite"; + else + retval = static_cast<double> (val); + } + else if (list->s_get_fcn) + { + retval = (odessa_opts.*list->s_get_fcn) (); + } + else + panic_impossible (); + + return retval; + } + list++; + } + + warning ("odessa_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (odessa_options, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} odessa_options (@var{opt}, @var{val})\n\ +When called with two arguments, this function allows you set options\n\ +parameters for the function @code{odessa}. Given one argument,\n\ +@code{odessa_options} returns the value of the corresponding option. If\n\ +no arguments are supplied, the names of all the available options and\n\ +their current values are displayed.\n\ +@end deftypefn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_odessa_option_list (octave_stdout); + } + else if (nargin == 1 || nargin == 2) + { + std::string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + retval = show_odessa_option (keyword); + else + { + if (args(1).is_string ()) + { + std::string val = args(1).string_value (); + + if (! error_state) + set_odessa_option (keyword, val); + } + else if (args(1).is_scalar_type ()) + { + double val = args(1).double_value (); + + if (! error_state) + set_odessa_option (keyword, val); + } + else + { + Array<double> val = args(1).vector_value (); + + if (! error_state) + set_odessa_option (keyword, val); + } + } + } + } + else + print_usage ("odessa_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -39,12 +39,13 @@ endif endif -DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc daspk.cc \ - dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \ - fsolve.cc gammainc.cc getgrent.cc getpwent.cc getrusage.cc \ - givens.cc hess.cc ifft.cc ifft2.cc inv.cc kron.cc log.cc lpsolve.cc \ - lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qz.cc rand.cc \ - schur.cc sort.cc svd.cc syl.cc time.cc +DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc \ + daspk.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \ + filter.cc find.cc fsolve.cc gammainc.cc getgrent.cc \ + getpwent.cc getrusage.cc givens.cc hess.cc ifft.cc ifft2.cc \ + inv.cc kron.cc log.cc lpsolve.cc lsode.cc lu.cc minmax.cc \ + odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \ + sort.cc svd.cc syl.cc time.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))