Mercurial > hg > octave-avbm
changeset 3990:46388d6a4e44
[project @ 2002-07-16 06:20:39 by jwe]
author | jwe |
---|---|
date | Tue, 16 Jul 2002 06:20:40 +0000 |
parents | bdde4f33221e |
children | 48d2bc4a3729 |
files | liboctave/ChangeLog liboctave/DAE.cc liboctave/DAE.h liboctave/DAEFunc.h liboctave/DAERT.h liboctave/DAERTFunc.h liboctave/DASPK.cc liboctave/DASRT.cc liboctave/DASRT.h liboctave/DASSL.cc liboctave/Makefile.in liboctave/base-dae.h src/ChangeLog src/DLD-FUNCTIONS/dasrt.cc src/Makefile.in |
diffstat | 15 files changed, 1977 insertions(+), 113 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,15 @@ +2002-07-16 John W. Eaton <jwe@bevo.che.wisc.edu> + + * DAE.cc: Delete. + + * DAERT.h, DAERTFunc.h, DASRT.h, DASRT.cc: New files for DAE + solving with root finding. + * Makefile.in: Add them to the appropriate lists. + + * base-dae.h: New file. + * Makefile.in (INCLUDES): Add it to the list. + * DAE.h (DAE): Derive from base_diff_alg_eqn, not base_diff_eqn. + 2002-07-10 John W. Eaton <jwe@bevo.che.wisc.edu> * ODE.h: Move integrate and do_integrate method declarations and
deleted file mode 100644 --- a/liboctave/DAE.cc +++ /dev/null @@ -1,68 +0,0 @@ -/* - -Copyright (C) 1996, 1997 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 "DAE.h" -#include "lo-error.h" - -DAE::DAE (const ColumnVector& xx, const ColumnVector& xxdot, - double tt, DAEFunc& f) - : base_diff_eqn (xx, tt), DAEFunc (f), xdot (xxdot) -{ - if (x.length () != xdot.length ()) - ; // XXX FIXME XXX -- exception! -} - -void -DAE::initialize (const ColumnVector& xx, double tt) -{ - if (xx.length () != xdot.length ()) - ; // XXX FIXME XXX -- exception! - else - base_diff_eqn::initialize (xx, tt); -} - -void -DAE::initialize (const ColumnVector& xx, const ColumnVector& xxdot, - double tt) -{ - if (xx.length () != xxdot.length ()) - ; // XXX FIXME XXX -- exception! - else - { - base_diff_eqn::initialize (xx, tt); - xdot = xxdot; - } -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/liboctave/DAE.h +++ b/liboctave/DAE.h @@ -23,54 +23,38 @@ #if !defined (octave_DAE_h) #define octave_DAE_h 1 -#if defined (__GNUG__) -#pragma interface -#endif - #include "DAEFunc.h" -#include "base-de.h" +#include "base-dae.h" class -DAE : public base_diff_eqn, public DAEFunc +DAE : public base_diff_alg_eqn, public DAEFunc { public: DAE (void) - : base_diff_eqn (), DAEFunc (), xdot () { } + : base_diff_alg_eqn (), DAEFunc () { } DAE (const ColumnVector& xx, double tt, DAEFunc& f) - : base_diff_eqn (xx, tt), DAEFunc (f), xdot (x.capacity (), 0.0) { } + : base_diff_alg_eqn (xx, tt), DAEFunc (f) { } DAE (const ColumnVector& xx, const ColumnVector& xxdot, - double tt, DAEFunc& f); + double tt, DAEFunc& f) + : base_diff_alg_eqn (xx, xxdot, tt), DAEFunc (f) { } DAE (const DAE& a) - : base_diff_eqn (a), DAEFunc (a), xdot (a.xdot) { } + : base_diff_alg_eqn (a), DAEFunc (a){ } DAE& operator = (const DAE& a) { if (this != &a) { - base_diff_eqn::operator = (a); + base_diff_alg_eqn::operator = (a); DAEFunc::operator = (a); - - xdot = a.xdot; } return *this; } ~DAE (void) { } - - ColumnVector state_derivative (void) { return xdot; } - - void initialize (const ColumnVector& xx, double tt); - - void initialize (const ColumnVector& xx, const ColumnVector& xxdot, - double tt); - -protected: - - ColumnVector xdot; }; #endif
--- a/liboctave/DAEFunc.h +++ b/liboctave/DAEFunc.h @@ -39,10 +39,10 @@ typedef ColumnVector (*DAERHSFunc) (const ColumnVector& x, const ColumnVector& xdot, - double, int&); + double t, int& ires); typedef DAEJac (*DAEJacFunc) (const ColumnVector& x, - const ColumnVector& xdot, double); + const ColumnVector& xdot, double t); DAEFunc (void) : fun (0), jac (0) { }
new file mode 100644 --- /dev/null +++ b/liboctave/DAERT.h @@ -0,0 +1,70 @@ +/* + +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_DAERT_h) +#define octave_DAERT_h 1 + +#include "DAE.h" +#include "DAERTFunc.h" +#include "base-dae.h" + +class +DAERT : public base_diff_alg_eqn, public DAERTFunc +{ +public: + + DAERT (void) + : base_diff_alg_eqn (), DAERTFunc () { } + + DAERT (const ColumnVector& x, const ColumnVector& xdot, double t, + DAERTFunc& f) + : base_diff_alg_eqn (x, xdot, t), DAERTFunc (f) { } + + DAERT (const DAERT& a) + : base_diff_alg_eqn (a), DAERTFunc (a) { } + + DAERT& operator = (const DAERT& a) + { + if (this != &a) + { + base_diff_alg_eqn::operator = (a); + DAERTFunc::operator = (a); + + } + return *this; + } + + ~DAERT (void) { } + + void initialize (const ColumnVector& x, const ColumnVector& xdot, double t) + { + base_diff_alg_eqn::initialize (x, xdot, t); + } +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/DAERTFunc.h @@ -0,0 +1,84 @@ +/* + +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_DAERTFunc_h) +#define octave_DAERTFunc_h 1 + +#include "dMatrix.h" + +class +DAERTFunc : DAEFunc +{ +public: + + typedef ColumnVector (*DAERTConstrFunc) (const ColumnVector& x, double t); + + DAERTFunc (void) + : DAEFunc (), constr (0) { } + + DAERTFunc (DAERHSFunc f) + : DAEFunc (f), constr (0) { } + + DAERTFunc (DAERHSFunc f, DAEJacFunc j) + : DAEFunc (f, j), constr (0) { } + + DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf) + : DAEFunc (f), constr (cf) { } + + DAERTFunc (DAERHSFunc f, DAERTConstrFunc cf, DAEJacFunc j) + : DAEFunc (f, j), constr (cf) { } + + DAERTFunc (const DAERTFunc& a) + : DAEFunc (a), constr (a.constr) { } + + DAERTFunc& operator = (const DAERTFunc& a) + { + if (this != &a) + { + DAEFunc::operator = (a); + constr = a.constr; + } + return *this; + } + + ~DAERTFunc (void) { } + + DAERTConstrFunc constraint_function (void) const { return constr; } + + DAERTFunc& set_constraint_function (DAERTConstrFunc cf) + { + constr = cf; + return *this; + } + +protected: + + DAERTConstrFunc constr; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/DASPK.cc +++ b/liboctave/DASPK.cc @@ -195,7 +195,7 @@ // Fix up the matrix of partial derivatives for daspk. - tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot; + tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot; for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++)
new file mode 100644 --- /dev/null +++ b/liboctave/DASRT.cc @@ -0,0 +1,667 @@ +/* + +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 <fstream.h> + +#include <cstdlib> +#include <cfloat> +#include <cmath> +#include "defun-dld.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" + +// For instantiating the Array<Matrix> object. +#include "Array.h" +#include "Array.cc" + +#include "DASRT.h" +#include "f77-fcn.h" +#include "lo-error.h" + +#ifndef F77_FUNC +#define F77_FUNC(x, X) F77_FCN (x, X) +#endif + +extern "C" +{ + int F77_FUNC (ddasrt, DASRT) (int (*)(const double&, double*, double*, + double*, int&, double*, int*), + const int&, const double&, double*, double*, + const double&, int*, double*, double*, + int&, double*, const int&, int*, + const int&, double*, int*, + int (*)(const double&, double*, + double*, double*, + const double&, double*, int*), + int (*)(const int&, const double&, double*, + const int&, double*, double*, int*), + const int&, int*); +} + +template class Array<Matrix>; + +static DAEFunc::DAERHSFunc user_fsub; +static DAEFunc::DAEJacFunc user_jsub; +static DAERTFunc::DAERTConstrFunc user_csub; +static int nn; + +static int +ddasrt_f (const double& t, double *state, double *deriv, double *delta, + int& ires, double *rpar, int *ipar) +{ + ColumnVector tmp_state (nn); + for (int i = 0; i < nn; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_deriv (nn); + for (int i = 0; i < nn; i++) + tmp_deriv(i) = deriv[i]; + + ColumnVector tmp_fval = user_fsub (tmp_state, tmp_deriv, t, ires); + + if (tmp_fval.length () == 0) + ires = -2; + else + { + for (int i = 0; i < nn; i++) + delta[i] = tmp_fval(i); + } + + return 0; +} + +//typedef int (*efptr) (const double& t, const int& n, double *state, +// double *ework, double *rpar, int *ipar, +// const int& ieform, int& ires); + +//static efptr e_fun; + +static int +ddasrt_j (const double& t, double *state, double *deriv, + double *pdwork, const double& cj, double *rpar, int *ipar) +{ + ColumnVector tmp_state (nn); + for (int i = 0; i < nn; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_deriv (nn); + for (int i = 0; i < nn; i++) + tmp_deriv(i) = deriv[i]; + + // XXX FIXME XXX + + Matrix tmp_dfdxdot (nn, nn); + Matrix tmp_dfdx (nn, nn); + + DAEFunc::DAEJac tmp_jac; + tmp_jac.dfdxdot = &tmp_dfdxdot; + tmp_jac.dfdx = &tmp_dfdx; + + tmp_jac = user_jsub (tmp_state, tmp_deriv, t); + + // Fix up the matrix of partial derivatives for dasrt. + + tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot; + + for (int j = 0; j < nn; j++) + for (int i = 0; i < nn; i++) + pdwork[j*nn+i] = tmp_dfdx.elem (i, j); + + return 0; +} + +static int +ddasrt_g (const int& neq, const double& t, double *state, const int& ng, + double *gout, double *rpar, int *ipar) +{ + int n = neq; + + ColumnVector tmp_state (n); + for (int i = 0; i < n; i++) + tmp_state(i) = state[i]; + + ColumnVector tmp_fval = user_csub (tmp_state, t); + + for (int i = 0; i < ng; i++) + gout[i] = tmp_fval(i); + + return 0; +} + + +DASRT::DASRT (void) + : DAERT () +{ + initialized = false; + restart = false; + + stop_time_set = false; + stop_time = 0.0; + + sanity_checked = false; + + info.resize (30, 0); + + npar = 0; + + liw = 0; + lrw = 0; +} + +DASRT::DASRT (const int& ng, const ColumnVector& state, + const ColumnVector& deriv, double time, DAERTFunc& f) + : DAERT (state, deriv, time, f) +{ + n = size (); + + initialized = false; + restart = false; + + stop_time_set = false; + stop_time = 0.0; + + DAERTFunc::operator = (f); + + sanity_checked = false; + + info.resize (30, 0); + jroot.resize (ng, 1); + + npar = 0; + + rpar.resize (npar+1); + ipar.resize (npar+1); + + info(11) = npar; + + // Also store it here, for communication with user-supplied + // subroutines. + ipar(0) = npar; + + y.resize (n, 1, 0.0); + ydot.resize (n, 1, 0.0); +} + +void +DASRT::init_work_size (int info_zero) +{ + double t; + double *py = y.fortran_vec (); + double *pydot = ydot.fortran_vec (); + double rel_tol = relative_tolerance (); + double abs_tol = absolute_tolerance (); + int *pinfo = info.fortran_vec (); + double *prpar = rpar.fortran_vec (); + int *pipar = ipar.fortran_vec (); + int *pjroot = jroot.fortran_vec (); + int idid; + + // We do not have to lie. + rwork.resize (5000+9*n+n*n, 0.0); + iwork.resize (n+20, 0); + + liw = n+20; + lrw = 5000+9*n+n*n; + + double *prwork = rwork.fortran_vec (); + int *piwork = iwork.fortran_vec (); + + + F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo, + &rel_tol, &abs_tol, idid, prwork, lrw, + piwork, liw, prpar, pipar, ddasrt_j, + ddasrt_g, ng, pjroot); + + int iwadd = iwork(18); + + + if (iwadd > 0) + liw += iwadd; + + info(0) = 0; + + iwork.resize (liw, 0); + + piwork = iwork.fortran_vec (); + + F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo, + &rel_tol, &abs_tol, idid, prwork, lrw, + piwork, liw, prpar, pipar, ddasrt_j, + ddasrt_g, ng, pjroot); + + int rwadd = iwork(19); + + + if (rwadd > 0) + lrw += rwadd; + + rwork.resize (lrw, 0.0); + + info(0) = info_zero; + +} + +void +DASRT::force_restart (void) +{ + restart = true; + integration_error = false; +} + +void +DASRT::set_stop_time (double t) +{ + stop_time_set = true; + stop_time = t; +} + +void +DASRT::set_ng (int the_ng) +{ + ng = the_ng; +} + +int DASRT::get_ng (void) +{ + return ng; +} + +void +DASRT::clear_stop_time (void) +{ + stop_time_set = false; +} + +void +DASRT::integrate (double tout) +{ + DASRT_result retval; + + if (! initialized) + { + info(0) = 0; + + for (int i = 0; i < n; i++) + { + y(i,0) = x(i); + ydot(i,0) = xdot(i); + } + + integration_error = false; + + user_fsub = DAEFunc::function (); + user_jsub = DAEFunc::jacobian_function (); + user_csub = DAERTFunc::constraint_function (); + + if (user_jsub) + info(4) = 1; + else + info(4) = 0; + + if (! sanity_checked) + { + int ires = 0; + + ColumnVector fval = user_fsub (x, xdot, t, ires); + + if (fval.length () != x.length ()) + { + (*current_liboctave_error_handler) + ("dassl: inconsistent sizes for state and residual vectors"); + + integration_error = true; + return; + } + + sanity_checked = true; + } + + + init_work_size (info(0)); + + + + + if (iwork.length () != liw) + iwork.resize (liw); + + if (rwork.length () != lrw) + rwork.resize (lrw); + + abs_tol = absolute_tolerance (); + rel_tol = relative_tolerance (); + + + 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; + + + py = y.fortran_vec (); + pydot = ydot.fortran_vec (); + pinfo = info.fortran_vec (); + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + prpar = rpar.fortran_vec (); + pipar = ipar.fortran_vec (); + pjroot = jroot.fortran_vec (); + + info(5) = 0; + info(8) = 0; + initialized = true; + } + + if (restart) + { + info(0) = 0; + + if (stop_time_set) + { + info(3) = 1; + rwork(0) = stop_time; + } + else + info(3) = 0; + } + + + + + F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo, + &rel_tol, &abs_tol, idid, prwork, lrw, + piwork, liw, prpar, pipar, ddasrt_j, + ddasrt_g, ng, pjroot)); + + if (f77_exception_encountered) + { + integration_error = true; + (*current_liboctave_error_handler) ("unrecoverable error in dassl"); + } + else + { + switch (idid) + { + case 0: // Initial conditions made consistent. + case 1: // A step was successfully taken in intermediate-output + // mode. The code has not yet reached TOUT. + case 2: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping exactly to TOUT. + case 3: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping past TOUT. Y(*) is obtained by + // interpolation. YPRIME(*) is obtained by interpolation. + case 5: // The integration to TSTOP was successfully completed + // (T=TSTOP) by stepping to TSTOP within the + // tolerance. Must restart to continue. + for (int i = 0; i < n; i++) + x(i) = y(i,0); + t = tout; + break; + + case 4: // We've hit the stopping condition. + for (int i = 0; i < n; i++) + x(i) = y(i,0); + break; + + case -1: // A large amount of work has been expended. (~500 steps). + case -2: // The error tolerances are too stringent. + case -3: // The local error test cannot be satisfied because you + // specified a zero component in ATOL and the + // corresponding computed solution component is zero. + // Thus, a pure relative error test is impossible for + // this component. + case -6: // DDASRT had repeated error test failures on the last + // attempted step. + case -7: // The corrector could not converge. + case -8: // The matrix of partial derivatives is singular. + case -9: // The corrector could not converge. There were repeated + // error test failures in this step. + case -10: // The corrector could not converge because IRES was + // equal to minus one. + case -11: // IRES equal to -2 was encountered and control is being + // returned to the calling program. + case -12: // DASSL failed to compute the initial YPRIME. + case -33: // The code has encountered trouble from which it cannot + // recover. A message is printed explaining the trouble + // and control is returned to the calling program. For + // example, this occurs when invalid input is detected. + default: + integration_error = true; + (*current_liboctave_error_handler) + ("ddasrt failed with IDID = %d", idid); + break; + } + } +} + +DASRT_result +DASRT::integrate (const ColumnVector& tout) +{ + DASRT_result retval; + + Matrix x_out; + Matrix xdot_out; + ColumnVector t_out; + + int oldj = 0; + + int n_out = tout.capacity (); + + + if (n_out > 0 && n > 0) + { + x_out.resize (n_out, n); + xdot_out.resize (n_out, n); + t_out.resize (n_out); + + for (int j = 0; j < n_out; j++) + { + integrate (tout(j)); + if (integration_error) + { + retval = DASRT_result (x_out, xdot_out, t_out); + return retval; + } + if (idid == 4) + t_out(j) = t; + else + t_out(j) = tout(j); + + + for (int i = 0; i < n; i++) + { + x_out(j,i) = y(i,0); + xdot_out(j,i) = ydot(i,0); + } + if (idid ==4) + { + oldj = j; + j = n_out; + x_out.resize (oldj+1, n); + xdot_out.resize (oldj+1, n); + t_out.resize (oldj+1); + } + } + } + + retval = DASRT_result (x_out, xdot_out, t_out); + + return retval; +} + +DASRT_result +DASRT::integrate (const ColumnVector& tout, const ColumnVector& tcrit) +{ + DASRT_result retval; + + Matrix x_out; + Matrix xdot_out; + ColumnVector t_outs; + + int n_out = tout.capacity (); + + if (n_out > 0 && n > 0) + { + x_out.resize (n_out, n); + xdot_out.resize (n_out, n); + t_outs.resize (n_out); + + 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 = 1; + 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 = 0; + i_crit++; + do_restart = true; + } + else + { + clear_stop_time (); + t_out = next_out; + save_output = 1; + i_out++; + } + } + else + { + set_stop_time (next_crit); + t_out = next_out; + save_output = 1; + i_out++; + } + + integrate (t_out); + + if (integration_error) + { + retval = DASRT_result (x_out, xdot_out, t_outs); + return retval; + } + if (idid == 4) + t_out = t; + + if (save_output) + { + for (int i = 0; i < n; i++) + { + x_out(i_out-1,i) = y(i,0); + xdot_out(i_out-1,i) = ydot(i,0); + } + t_outs(i_out-1) = t_out; + if (idid ==4) + { + x_out.resize (i_out, n); + xdot_out.resize (i_out, n); + t_outs.resize (i_out); + i_out = n_out; + } + + } + + if (do_restart) + force_restart (); + } + + retval = DASRT_result (x_out, xdot_out, t_outs); + } + else + { + retval = integrate (tout); + + if (integration_error) + return retval; + } + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/DASRT.h @@ -0,0 +1,235 @@ +/* + +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_DASRT_h) +#define octave_DASRT_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include <cfloat> +#include <cmath> + +#include "DAERT.h" + +class +DASRT_options +{ +public: + + DASRT_options (void) { init (); } + + DASRT_options (const DASRT_options& opt) { copy (opt); } + + DASRT_options& operator = (const DASRT_options& opt) + { + if (this != &opt) + copy (opt); + + return *this; + } + + ~DASRT_options (void) { } + + void 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; + x_step_limit = -1; + } + + void copy (const DASRT_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; + x_step_limit = opt.x_step_limit; + } + + void set_default_options (void) { init (); } + + void set_absolute_tolerance (double val) + { x_absolute_tolerance = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); } + + void set_initial_step_size (double val) + { x_initial_step_size = (val >= 0.0) ? val : -1.0; } + + 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 >= 0) ? val : -1; } + + double absolute_tolerance (void) { return x_absolute_tolerance; } + + double initial_step_size (void) { return x_initial_step_size; } + + double maximum_step_size (void) { return x_maximum_step_size; } + + double minimum_step_size (void) { return x_minimum_step_size; } + + double relative_tolerance (void) { return x_relative_tolerance; } + + int step_limit (void) { return x_step_limit; } + +private: + + double x_absolute_tolerance; + double x_initial_step_size; + double x_maximum_step_size; + double x_minimum_step_size; + double x_relative_tolerance; + int x_step_limit; +}; + +class +DASRT_result +{ +public: + + DASRT_result (void) { } + + DASRT_result (const Matrix& xx, const Matrix& xxdot, const ColumnVector& tt) + : x (xx), xdot (xxdot), t (tt) { } + + DASRT_result (const DASRT_result& r) + : x (r.x), xdot (r.xdot), t (r.t) { } + + DASRT_result& operator = (const DASRT_result& r) + { + if (this != &r) + { + x = r.x; + xdot = r.xdot; + t = r.t; + } + return *this; + } + + ~DASRT_result (void) { } + + Matrix state (void) const { return x; } + Matrix deriv (void) const { return xdot; } + ColumnVector times (void) const { return t; } + +private: + + Matrix x; + Matrix xdot; + ColumnVector t; +}; + +class +DASRT : public DAERT, public DASRT_options +{ +public: + + DASRT (void); + + DASRT (const int& ng, const ColumnVector& x, const ColumnVector& xdot, + double time, DAERTFunc& f); + + ~DASRT (void) { } + + void force_restart (void); + + void set_stop_time (double t); + void clear_stop_time (void); + void set_ng (int the_ng); + int get_ng (void); + + DASRT_result integrate (const ColumnVector& tout); + + DASRT_result integrate (const ColumnVector& tout, + const ColumnVector& tcrit); + +private: + + bool initialized; + + bool sanity_checked; + + bool stop_time_set; + double stop_time; + + bool restart; + + bool integration_error; + + int liw; + int lrw; + int idid; + int ieform; + int lun; + + int n; + int npar; + int ng; + + Array<int> info; + Array<int> iwork; + Array<int> ipar; + Array<int> jroot; + + Array<double> rwork; + Array<double> rpar; + + Matrix y; + Matrix ydot; + + double abs_tol; + double rel_tol; + + double *py; + double *pydot; + int *pinfo; + int *piwork; + double *prwork; + double *prpar; + int *pipar; + int *pjroot; + + void init_work_size (int); + + void integrate (double t); +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -181,7 +181,7 @@ // Fix up the matrix of partial derivatives for dassl. - tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot; + tmp_dfdx = tmp_dfdx + cj * tmp_dfdxdot; for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++)
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -44,18 +44,19 @@ vx-rv-cs.h vx-s-ccv.h vx-s-crv.h \ vx-rv-crv.h vx-cv-ccv.h vx-crv-rv.h vx-ccv-cv.h -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 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 \ - lo-utils.h mach-info.h oct-alloc.h oct-cmplx.h oct-env.h \ - oct-fftw.h oct-getopt.h oct-group.h oct-kpse.h oct-passwd.h \ - oct-rl-edit.h oct-rl-hist.h oct-shlib.h oct-syscalls.h oct-time.h \ - pathlen.h pathsearch.h prog-args.h statdefs.h str-vec.h sun-utils.h \ - sysdir.h systime.h syswait.h \ +INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DAERT.h DAERTFunc.h \ + DASPK.h DASRT.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 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 lo-utils.h mach-info.h \ + oct-alloc.h oct-cmplx.h oct-env.h oct-fftw.h oct-getopt.h \ + oct-group.h oct-kpse.h oct-passwd.h oct-rl-edit.h \ + oct-rl-hist.h oct-shlib.h oct-syscalls.h oct-time.h \ + pathlen.h pathsearch.h prog-args.h statdefs.h str-vec.h\ + sun-utils.h sysdir.h systime.h syswait.h \ $(MATRIX_INC) \ $(MX_OP_INC) \ $(VX_OP_INC) @@ -85,10 +86,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 ODES.cc ODESSA.cc Quad.cc \ - Range.cc data-conv.cc dir-ops.cc \ +LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc \ + DASPK.cc DASRT.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 \
new file mode 100644 --- /dev/null +++ b/liboctave/base-dae.h @@ -0,0 +1,86 @@ +/* + +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_base_dae_h) +#define octave_base_dae_h 1 + +#include "base-de.h" + +class +base_diff_alg_eqn : public base_diff_eqn +{ +public: + + base_diff_alg_eqn (void) + : base_diff_eqn (), xdot () { } + + base_diff_alg_eqn (const ColumnVector& xx, double tt) + : base_diff_eqn (xx, tt), xdot (xx.length (), 0.0) { } + + base_diff_alg_eqn (const ColumnVector& xx, const ColumnVector& xxdot, + double tt) + : base_diff_eqn (xx, tt), xdot (xxdot) { } + + base_diff_alg_eqn (const base_diff_alg_eqn& a) + : base_diff_eqn (a), xdot (a.xdot) { } + + virtual ~base_diff_alg_eqn (void) { } + + base_diff_alg_eqn& operator = (const base_diff_alg_eqn& a) + { + if (this != &a) + { + base_diff_eqn::operator = (a); + xdot = a.xdot; + } + return *this; + } + + void initialize (const ColumnVector& x0, double t0) + { + base_diff_eqn::initialize (x0, t0); + xdot.resize (x0.length (), 0.0); + force_restart (); + } + + void initialize (const ColumnVector& x0, const ColumnVector& xdot0, + double t0) + { + base_diff_eqn::initialize (x0, t0); + xdot = xdot0; + force_restart (); + } + + ColumnVector state_derivative (void) { return xdot; } + +protected: + + ColumnVector xdot; +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2002-07-16 John W. Eaton <jwe@bevo.che.wisc.edu> + + * DLD-FUNCTIONS/dasrt.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. + 2002-07-12 John W. Eaton <jwe@bevo.che.wisc.edu> * lex.l (@): Handle new token.
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/dasrt.cc @@ -0,0 +1,788 @@ +/* + +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 <iostream.h> + +#include "DASRT.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 "parse.h" +#include "unwind-prot.h" +#include "utils.h" +#include "variables.h" + +// Global pointers for user defined function required by dassl. +static octave_function *dasrt_f; +static octave_function *dasrt_j; +static octave_function *dasrt_cf; + +static DASRT_options dasrt_opts; + +// Is this a recursive call? +static int call_depth = 0; + +static ColumnVector +dasrt_user_f (const ColumnVector& x, const ColumnVector& xprime, + double t, int& ires) +{ + ColumnVector retval; + + octave_value_list args; + + int n = x.length (); + + if (n > 1) + { + args(0) = x; + args(1) = xprime; + } + else if (n == 1) + { + args(0) = x(0); + args(1) = xprime(0); + } + else + { + args(0) = Matrix (); + args(1) = Matrix (); + } + + args(2) = t; + + if (dasrt_f) + { + octave_value_list tmp = dasrt_f->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("dasrt"); + 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 ("dasrt"); + } + else + gripe_user_supplied_eval ("dasrt"); + } + + return retval; +} + +static ColumnVector +dasrt_user_cf (const ColumnVector& x, double t) +{ + ColumnVector retval; + + octave_value_list args; + + int n = x.length (); + + if (n > 1) + args(0) = x; + else if (n == 1) + args(0) = x(0); + else + args(0) = Matrix (); + + args(1) = t; + + if (dasrt_cf) + { + octave_value_list tmp = dasrt_cf->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("dasrt"); + 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 ("dasrt"); + } + else + gripe_user_supplied_eval ("dasrt"); + } + + return retval; +} + +static ColumnVector +dasrt_dumb_cf (const ColumnVector& x, double t) +{ + ColumnVector retval (1, 1.0); + return retval; +} + +#if 0 +static Matrix +dasrt_user_mf (double t, const ColumnVector& x, const ColumnVector& xprime, + const double& cj, octave_function *mf) +{ + Matrix retval; + + if (mf) + { + octave_value_list args; + + int n = x.length (); + + if (n > 1) + { + args(0) = x; + args(1) = xprime; + args(3) = cj; + } + else if (n == 1) + { + args(0) = x(0); + args(1) = xprime(0); + args(3) = cj; + } + else + { + args(0) = Matrix (); + args(1) = Matrix (); + args(3) = Matrix (); + } + + args(2) = t; + + octave_value_list tmp = mf->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("dasrt"); + 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 ("dasrt"); + } + else + gripe_user_supplied_eval ("dasrt"); + } + + return retval; +} + +static Matrix +dasrt_user_j (const ColumnVector& x, const ColumnVector& xprime, double t) +{ + return dasrt_user_mf (t, x, xprime, dasrt_j); +} +#endif + +#define DASRT_ABORT \ + do \ + { \ + unwind_protect::run_frame ("Fdasrt"); \ + return retval; \ + } \ + while (0) + +#define DASRT_ABORT1(msg) \ + do \ + { \ + ::error ("dasrt: " ## msg); \ + DASRT_ABORT; \ + } \ + while (0) + +#define DASRT_ABORT2(fmt, arg) \ + do \ + { \ + ::error ("dasrt: " ## fmt, arg); \ + DASRT_ABORT; \ + } \ + while (0) + +DEFUN_DLD (dasrt, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{t}] =} dasrt (@var{fj} [, @var{g}], @var{x_0}, @var{xdot_0}, @var{t_out} [, @var{t_crit}])\n\ +Solve a system of differential/algebraic equations with functional\n\ +stopping criteria.\n\ +\n\ +The function to be integrated must be of the form:\n\ +@example\n\ +@var{res} = f (@var{x}, @var{xdot}, @var{t}) = 0\n\ +@end example\n\ +\n\ +The stopping condition must be of the form:\n\ +\n\ +@example\n\ +@var{res} = g (@var{x}, @var{t}) = 0\n\ +@end example\n\ +\n\ +The Jacobian (if included) must be of the form:\n\ +\n\ +@example\n\ +@var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{cj})\n\ + = df/dx + cj*df/dxdot\n\ +@end example\n\ +\n\ +@noindent\n\ +The following inputs are entered:\n\ +\n\ +@table @var\n\ +@item fj\n\ +The string vector containing either @var{f} or both @var{f} and @var{j}.\n\ +\n\ +@item f\n\ +The function to be integrated.\n\ +\n\ +@item g\n\ +The function with the stopping conditions.\n\ +\n\ +@item j\n\ +The optional Jacobian function. If not included, it will be approximated\n\ +by finite differences.\n\ +\n\ +@item x_0\n\ +The initial state.\n\ +\n\ +@item xdot_0\n\ +The time derivative of the initial state.\n\ +\n\ +@item t_out\n\ +The times at which the solution should be returned. This vector should\n\ +include the initial time a the first value.\n\ +\n\ +@item t_crit\n\ +The times at which the function is non-smooth or poorly behaved.\n\ +\n\ +@end table\n\ +\n\ +@noindent\n\ +The following outputs are returned:\n\ +\n\ +@table @var\n\ +@item x\n\ +The states at each time in @var{t}.\n\ +\n\ +@item xdot\n\ +The time derivative of the states at each time in @var{t}.\n\ +\n\ +@item t\n\ +All the times in the requested output time vector up to the stopping\n\ +criteria. The time at which the stopping criteria is achieved is returned\n\ +as the last time in the vector.\n\ +@end table\n\ +\n\ +You can use the function @code{dasrt_options} to set optional\n\ +parameters for @code{dasrt}.\n\ +@end deftypefn") +{ + octave_value_list retval; + + unwind_protect::begin_frame ("Fdasrt"); + + unwind_protect_int (call_depth); + call_depth++; + + if (call_depth > 1) + DASRT_ABORT1 ("invalid recursive call"); + + int argp = 0; + int ng = 0; + + int nargin = args.length (); + + if (nargin < 4 || nargin > 6) + { + print_usage ("dasrt"); + unwind_protect::run_frame ("Fdasrt"); + return retval; + } + + dasrt_f = 0; + dasrt_j = 0; + dasrt_cf = 0; + + // Check all the arguments. Are they the right animals? + + // Here's where I take care of f and j in one shot: + + octave_value f_arg = args(0); + + switch (f_arg.rows ()) + { + case 1: + dasrt_f = extract_function + (args(0), "dasrt", "__dasrt_fcn__", + "function res = __dasrt_fcn__ (x, xdot, t) res = ", + "; endfunction"); + break; + + case 2: + { + string_vector tmp = args(0).all_strings (); + + if (! error_state) + { + dasrt_f = extract_function + (tmp(0), "dasrt", "__dasrt_fcn__", + "function res = __dasrt_fcn__ (x, xdot, t) res = ", + "; endfunction"); + + if (dasrt_f) + { + dasrt_j = extract_function + (tmp(1), "dasrt", "__dasrt_jac__", + "function jac = __dasrt_jac__ (x, xdot, t, cj) jac = ", + "; endfunction"); + + if (! dasrt_j) + dasrt_f = 0; + } + } + } + break; + + default: + DASRT_ABORT1 + ("first arg should be a string or 2-element string array"); + } + + if (error_state || (! dasrt_f)) + DASRT_ABORT; + + DAERTFunc func (dasrt_user_f); + + argp++; + double t0; + ColumnVector x0; + + if (args(1).is_string ()) + { + dasrt_cf = is_valid_function (args(1), "dasrt", true); + + if (! dasrt_cf) + DASRT_ABORT1 ("expecting function name as argument 2"); + else + { + octave_value_list blah; + octave_value_list inputs; + ColumnVector the_ts; + x0 = ColumnVector (args(2).vector_value ()); + + if (error_state) + DASRT_ABORT1 ("bad x0"); + + the_ts = ColumnVector (args(4).vector_value ()); + + if (error_state) + DASRT_ABORT1 ("bad tout"); + + t0 = the_ts(0); + + inputs(0) = x0; + inputs(1) = t0; + + string g_name = args(1).string_value(); + + if (error_state) + DASRT_ABORT1 ("bad function name of stopping condition"); + + blah = feval (g_name, inputs, 2); + ColumnVector gout = ColumnVector (blah(0).vector_value ()); + + ng = gout.length(); + int testg; + for (testg = 0; testg < ng; testg++) + { + if (gout(testg) == 0) + DASRT_ABORT1 ("stopping condition satisfied at initial point"); + } + } + argp++; + + func.set_constraint_function (dasrt_user_cf); + } + else + { + // Now this second argument is not a string. It has to be x0. + // we call the dummy g function and set ng = 1; + + x0 = ColumnVector (args(1).vector_value ()); + + if (error_state) + DASRT_ABORT1 ("bad x0"); + + func.set_constraint_function (dasrt_dumb_cf); + + ng = 1; + } + + ColumnVector state (args(argp++).vector_value ()); + + if (error_state) + DASRT_ABORT2 ("expecting state vector as argument %d", argp); + + ColumnVector stateprime (args(argp++).vector_value()); + + if (error_state) + DASRT_ABORT2 + ("expecting time derivative of state vector as argument %d", argp); + + ColumnVector old_out_times (args(argp++).vector_value ()); + + if (error_state) + DASRT_ABORT2 + ("expecting output time vector as %s argument %d", argp); + + double tzero = old_out_times (0); + + int ol = old_out_times.length (); + int ijk; + + ColumnVector out_times (ol-1, 0.0); + + for (ijk = 1; ijk < ol; ijk++) + out_times(ijk-1) = old_out_times(ijk); + + ColumnVector crit_times; + + bool crit_times_set = false; + + if (argp < nargin) + { + crit_times = ColumnVector (args(argp++).vector_value ()); + + if (error_state) + DASRT_ABORT2 + ("expecting critical time vector as argument %d", argp); + + crit_times_set = true; + } + +#if 0 + if (dasrt_j) + func.set_jacobian_function (dasrt_user_j); +#endif + + DASRT_result output; + + DASRT dae = DASRT (ng, state, stateprime, tzero, func); + + dae.copy (dasrt_opts); + dae.set_ng (ng); + + if (error_state) + DASRT_ABORT1 ("something is wrong"); + + if (crit_times_set) + output = dae.integrate (out_times, crit_times); + else + output = dae.integrate (out_times); + + if (! error_state) + { + ColumnVector old_output_times = output.times (); + + Matrix old_output_deriv = output.deriv (); + Matrix old_output_state = output.state (); + + int lstuff = old_output_times.length (); + int lstate = x0.length (); + + ColumnVector output_times (lstuff+1, 0.0); + + Matrix output_deriv (lstuff+1, lstate, 0.0); + Matrix output_state (lstuff+1, lstate, 0.0); + + output_times(0) = tzero; + + for (ijk = 0; ijk < lstate; ijk++) + { + output_deriv(0,ijk) = stateprime(ijk); + output_state(0,ijk) = state(ijk); + } + + for (ijk = 0; ijk < lstuff; ijk++) + { + output_times(ijk+1) = old_output_times(ijk); + + for (int lmnop=0; lmnop < lstate; lmnop++) + { + output_deriv(ijk+1,lmnop) = old_output_deriv(ijk,lmnop); + output_state(ijk+1,lmnop) = old_output_state(ijk,lmnop); + } + } + + retval(2) = output_times; + retval(1) = output_deriv; + retval(0) = output_state; + } + else + { + DASRT_ABORT1("something wicked has occurred!"); + // print_usage ("dasrt"); + } + + unwind_protect::run_frame ("Fdasrt"); + + return retval; +} + +typedef void (DASRT_options::*d_set_opt_mf) (double); +typedef void (DASRT_options::*i_set_opt_mf) (int); +typedef double (DASRT_options::*d_get_opt_mf) (void); +typedef int (DASRT_options::*i_get_opt_mf) (void); + +#define MAX_TOKENS 3 + +struct DASRT_OPTIONS +{ + const char *keyword; + const 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 DASRT_OPTIONS dasrt_option_table [] = +{ + { "absolute tolerance", + { "absolute", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + &DASRT_options::set_absolute_tolerance, 0, + &DASRT_options::absolute_tolerance, 0, }, + + { "initial step size", + { "initial", "step", "size", 0, }, + { 1, 0, 0, 0, }, 1, + &DASRT_options::set_initial_step_size, 0, + &DASRT_options::initial_step_size, 0, }, + + { "maximum step size", + { "maximum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + &DASRT_options::set_maximum_step_size, 0, + &DASRT_options::maximum_step_size, 0, }, + + { "minimum step size", + { "minimum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + &DASRT_options::set_minimum_step_size, 0, + &DASRT_options::minimum_step_size, 0, }, + + { "relative tolerance", + { "relative", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + &DASRT_options::set_relative_tolerance, 0, + &DASRT_options::relative_tolerance, 0, }, + + { "step limit", + { "step", "limit", 0, 0, }, + { 1, 0, 0, 0, }, 1, + 0, &DASRT_options::set_step_limit, + 0, &DASRT_options::step_limit, }, + + { 0, + { 0, 0, 0, 0, }, + { 0, 0, 0, 0, }, 0, + 0, 0, 0, 0, }, +}; + +static void +print_dasrt_option_list (ostream& os) +{ + print_usage ("dasrt_options", 1); + + os << "\n" + << "Options for dasrt include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + DASRT_OPTIONS *list = dasrt_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + + if (list->d_get_fcn) + { + double val = (dasrt_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } + else + { + int val = (dasrt_opts.*list->i_get_fcn) (); + if (val < 0) + os << "computed automatically"; + else + os << val; + } + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_dasrt_option (const string& keyword, double val) +{ + DASRT_OPTIONS *list = dasrt_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->d_set_fcn) + (dasrt_opts.*list->d_set_fcn) (val); + else + { + if (xisnan (val)) + { + error ("dasrt_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (dasrt_opts.*list->i_set_fcn) (NINT (val)); + } + return; + } + list++; + } + + warning ("dasrt_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_dasrt_option (const string& keyword) +{ + octave_value retval; + + DASRT_OPTIONS *list = dasrt_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->d_get_fcn) + { + double val = (dasrt_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } + else + { + int val = (dasrt_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "computed automatically"; + else + retval = static_cast<double> (val); + } + + return retval; + } + list++; + } + + warning ("dasrt_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (dasrt_options, args, , + "dasrt_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for dasrt. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_dasrt_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_dasrt_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_dasrt_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("dasrt_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -40,7 +40,7 @@ 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 \ + daspk.cc dasrt.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 \