From octave-sources-request at bevo dot che dot wisc dot edu  Wed Nov 24 13:56:23 1999
Subject: Modified lsode that passes extra params to xdot function
From: David Thompson <dcthomp at mail dot utexas dot edu>
To: octave-sources at bevo dot che dot wisc dot edu, dcthomp@mail.utexas.edu
Date: Wed, 24 Nov 1999 14:07:52 -0600

This is a multi-part message in MIME format.
--------------7A0668A37A343E9BB6E578A3
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

I've enclosed a modified version of lsode.cc that accepts
extra parameters and passes them to the rate of state user
function. I did this to avoid global variables because I wanted
to pass constant parameters to a dld module version of a state
equation; looking up global variable values at each time step
and/or recompiling every time I changed constants seemed like a
mess.

The new function is named lsode_extra so that you can use
the existing lsode function at the same time. To compile
this code:
	mkoctfile lsode_extra.cc -DHAVE_CONFIG_H
or, for the stable version of octave (2.0.14), try:
	mkoctfile lsode_extra.cc -DHAVE_CONFIG_H -DOCTAVE_STABLE_REV
I haven't tested with Octave 2.0.14. Let me know if it works.

		David Thompson
--------------7A0668A37A343E9BB6E578A3
Content-Type: text/plain; charset=us-ascii;
 name="lsode_extra.cc"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="lsode_extra.cc"

/*

Copyright (C) 1996, 1997 John W. Eaton
Modifications (C) 1999 David Thompson <dc at thompson dot cx>, assigned to J. W. Eaton.

This 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.

This 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 <octave/config.h>
#endif

#include <string>

#include <iostream.h>

#include <octave/LSODE.h>
#include <octave/lo-mappers.h>

#include <octave/defun-dld.h>
#include <octave/error.h>
#include <octave/gripes.h>
#include <octave/oct-obj.h>
#include <octave/pager.h>
#include <octave/utils.h>
#include <octave/variables.h>

#ifdef OCTAVE_STABLE_REV
#include <octave/pt-fvc.h>
#else
#include <octave/ov-fcn.h>
#include <octave/unwind-prot.h>
#endif // ! OCTAVE_STABLE_REV

#ifdef OCTAVE_STABLE_REV
// Global pointer for user defined function required by lsode.
static tree_fvc* lsodex_fcn ;
// Global pointer for optional user defined jacobian function used by lsode.
static tree_fvc* lsodex_jac ;
#else
// Global pointer for user defined function required by lsode.
static octave_function *lsodex_fcn;
// Global pointer for optional user defined jacobian function used by lsode.
static octave_function *lsodex_jac;
#endif

static LSODE_options lsodex_opts;

#ifndef OCTAVE_STABLE_REV
static int call_depth = 0;	// Is this a recursive call?
#endif // ! OCTAVE_STABLE_REV

static octave_value_list extras;
static int numextras;

ColumnVector
lsodex_user_function (const ColumnVector& x, double t)
{
  ColumnVector retval;

  int nstates = x.capacity ();

  octave_value_list args;
  args.resize (2+numextras);
  args(1) = t;

  Matrix m (nstates, 1);
  for (int i = 0; i < nstates; i++)
    m (i, 0) = x (i);
  octave_value state (m);
  args(0) = state;

  if (numextras) {
    for(int i=0;i<numextras;i++) {
      args(i+2)=extras(i);
    }
  }

  if (lsodex_fcn)
    {
#ifdef OCTAVE_STABLE_REV
      octave_value_list tmp = lsodex_fcn->eval( 1, args ) ;
#else
      octave_value_list tmp = lsodex_fcn->do_index_op (1, args);
#endif

      if (error_state)
	{
	  gripe_user_supplied_eval ("lsode_extra");
	  return retval;
	}

      if (tmp.length () > 0 && tmp(0).is_defined ())
	{
	  retval = tmp(0).vector_value ();

	  if (error_state || retval.length () == 0)
	    gripe_user_supplied_eval ("lsode_extra");
	}
      else
	gripe_user_supplied_eval ("lsode_extra");
    }

  return retval;
}

Matrix
lsodex_user_jacobian (const ColumnVector& x, double t)
{
  Matrix retval;

  int nstates = x.capacity ();

  octave_value_list args;
  args.resize (2+numextras);
  args(1) = t;

  Matrix m (nstates, 1);
  for (int i = 0; i < nstates; i++)
    m (i, 0) = x (i);
  octave_value state (m);
  args(0) = state;

  if (numextras) {
    for(int i=0;i<numextras;i++) {
      args(i+2)=extras(i);
    }
  }

  if (lsodex_jac)
    {
#ifdef OCTAVE_STABLE_REV
      octave_value_list tmp = lsodex_jac->eval( 1, args ) ;
#else
      octave_value_list tmp = lsodex_jac->do_index_op (1, args);
#endif

      if (error_state)
	{
	  gripe_user_supplied_eval ("lsode_extra");
	  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 ("lsode_extra");
	}
      else
	gripe_user_supplied_eval ("lsode_extra");
    }

  return retval;
}

#ifdef OCTAVE_STABLE_REV
#define LSODE_ABORT() \
  do \
    { \
      return retval; \
    } \
  while (0)
#else // OCTAVE_STABLE_REV !
#define LSODE_ABORT() \
  do \
    { \
      unwind_protect::run_frame ("Flsodex"); \
      return retval; \
    } \
  while (0)
#endif // ! OCTAVE_STABLE_REV
 
#define LSODE_ABORT1(msg) \
  do \
    { \
      ::error ("lsode_extra: " ## msg); \
      LSODE_ABORT (); \
    } \
  while (0)

#define LSODE_ABORT2(fmt, arg) \
  do \
    { \
      ::error ("lsode_extra: " ## fmt, arg); \
      LSODE_ABORT (); \
    } \
  while (0)

DEFUN_DLD (lsode_extra, args, nargout,
  "XDOT = lsode_extra (F, X0, T_OUT, T_CRIT [, EXTRA...] )\n\
\n\
The first argument is the name of the function to call to\n\
compute the vector of right hand sides.  It must have the form\n\
\n\
  XDOT = F (X, t [, EXTRA...] )\n\
\n\
where XDOT and X are vectors and t is a scalar.\n\
The optional EXTRA arguments will be passed to the function F\n\
and are intended for passing time-constant parameters to the\n\
differential equation without using global variables, which are evil.\n\
")
{
  octave_value_list retval;

#ifndef OCTAVE_STABLE_REV
  unwind_protect::begin_frame ("Flsodex");

  unwind_protect_int (call_depth);
  call_depth++;

  if (call_depth > 1)
    LSODE_ABORT1 ("invalid recursive call");
#endif // ! OCTAVE_STABLE_REV

  int nargin = args.length ();

  if (nargin > 2 && nargout < 2)
    {
      octave_value f_arg = args(0);

      switch (f_arg.rows ())
	{
	case 1:
#ifdef OCTAVE_STABLE_REV
	  lsode_fcn = is_valid_function(args(0), "lsode_extra", 1) ;
#else
	  lsodex_fcn = extract_function
	    (args(0), "lsode_extra", "__lsode_extra_fcn__",
	     "function xdot = __lsode_extra_fcn__ (x, t, p) xdot = ",
	     "; endfunction");
#endif // ! OCTAVE_STABLE_REV
	  break;

	case 2:
	  {
	    string_vector tmp = args(0).all_strings ();

	    if (! error_state)
	      {
#ifdef OCTAVE_STABLE_REV
		lsode_fcn = is_valid_function(tmp(0), "lsode_extra", 1) ;
#else
		lsodex_fcn = extract_function
		  (tmp(0), "lsode_extra", "__lsode_extra_fcn__",
		   "function xdot = __lsode_extra_fcn__ (x, t, p) xdot = ",
		   "; endfunction");
#endif // ! OCTAVE_STABLE_REV

		if (lsodex_fcn)
		  {
#ifdef OCTAVE_STABLE_REV
		    lsode_jac = is_valid_function(tmp(1), "lsode_extra", 1) ;
#else
		    lsodex_jac = extract_function
		      (tmp(1), "lsode_extra", "__lsode_extra_jac__",
		       "function jac = __lsode_extra_jac__ (x, t, p) jac = ",
		       "; endfunction");
#endif // ! OCTAVE_STABLE_REV

		    if (! lsodex_jac)
		      lsodex_fcn = 0;
		  }
	      }
	  }
	  break;

	default:
	  LSODE_ABORT1
	    ("first arg should be a string or 2-element string array");
	}

      if (error_state || ! lsodex_fcn)
	LSODE_ABORT ();

      ColumnVector state = args(1).vector_value ();

      if (error_state)
	LSODE_ABORT1 ("expecting state vector as second argument");

      ColumnVector out_times = args(2).vector_value ();

      if (error_state)
	LSODE_ABORT1 ("expecting output time vector as third argument");

      ColumnVector crit_times;

      int crit_times_set = 0;
      if (nargin > 3)
	{
	  crit_times = args(3).vector_value ();

	  if (error_state)
	    LSODE_ABORT1 ("expecting critical time vector as fourth argument");

	  if ( crit_times.capacity() > 0 )
	    crit_times_set = 1;
	}

      if (nargin > 4)
	{ // extras will pass the extra arguments through
	  extras.resize(nargin-4);
	  numextras=nargin-4;

	  for(int i=4;i<nargin;i++)
	    extras(i-4)=args(i);

	} 
      else 
	{
	  extras.resize(0);
	  numextras=0;
	}

      double tzero = out_times (0);
      int nsteps = out_times.capacity ();

      ODEFunc func (lsodex_user_function);
      if (lsodex_jac)
	func.set_jacobian_function (lsodex_user_jacobian);

      LSODE ode (state, tzero, func);

      ode.copy (lsodex_opts);

      int nstates = state.capacity ();
      Matrix output (nsteps, nstates + 1);

      if (crit_times_set)
	output = ode.integrate (out_times, crit_times);
      else
	output = ode.integrate (out_times);

      if (! error_state)
	{
	  retval.resize (1);
	  retval(0) = output;
	}
    }
  else
    print_usage ("lsode_extra");

#ifndef OCTAVE_STABLE_REV
  unwind_protect::run_frame ("Flsodex");
#endif // ! OCTAVE_STABLE_REV

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/

--------------7A0668A37A343E9BB6E578A3--