changeset 2850:b7f43611d1e8

[project @ 1997-03-28 21:36:35 by jwe]
author jwe
date Fri, 28 Mar 1997 21:39:35 +0000
parents 5338beb20eb9
children b960bd6cbfdf
files liboctave/ChangeLog liboctave/LSODE.cc liboctave/LSODE.h src/ChangeLog src/lsode.cc
diffstat 5 files changed, 192 insertions(+), 48 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,13 @@
+Fri Mar 28 15:37:09 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* LSODE.h (x_step_limit): New field.
+	(LSODE_options::init): Initialize it.
+	(LSODE_options::copy): Copy it.
+	(LSODE_options::set_step_limit, LSODE_options::step_limit):
+	New functions.
+	(LSODE::working_too_hard): Delete.
+	* LSODE.cc (LSODE::do_integrate): Handle step limit.
+
 Wed Mar 26 15:31:57 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* MArray-b.cc: Delete.
--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -239,11 +239,12 @@
   rwork.elem (5) = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0;
   rwork.elem (6) = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0;
 
+  if (step_limit () > 0)
+    iwork.elem (5) = step_limit ();
+
   int *piwork = iwork.fortran_vec ();
   double *prwork = rwork.fortran_vec ();
 
-  working_too_hard = 0;
-
  again:
 
   F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol,
@@ -268,12 +269,11 @@
 	  break;
 
 	case -1:  // excess work done on this call (perhaps wrong mf).
-	  working_too_hard++;
-	  if (working_too_hard > 20)
+	  if (step_limit () > 0)
 	    {
 	      (*current_liboctave_error_handler)
 		("giving up after more than %d steps attempted in lsode",
-		 iwork.elem (5) * 20);
+		 step_limit ());
 	      integration_error = 1;
 	    }
 	  else
--- a/liboctave/LSODE.h
+++ b/liboctave/LSODE.h
@@ -64,6 +64,7 @@
       x_maximum_step_size = -1.0;
       x_minimum_step_size = 0.0;
       x_relative_tolerance = sqrt_eps;
+      x_step_limit = 500;
     }
 
   void copy (const LSODE_options& opt)
@@ -73,6 +74,7 @@
       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 (); }
@@ -92,6 +94,9 @@
   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; }
+
   double absolute_tolerance (void)
     { return x_absolute_tolerance; }
 
@@ -107,6 +112,9 @@
   double relative_tolerance (void)
     {  return x_relative_tolerance; }
 
+  int step_limit (void)
+    { return x_step_limit; }
+
 private:
 
   double x_absolute_tolerance;
@@ -114,6 +122,8 @@
   double x_maximum_step_size;
   double x_minimum_step_size;
   double x_relative_tolerance;
+
+  int x_step_limit;
 };
 
 class
@@ -164,7 +174,6 @@
   int iopt;
   int liw;
   int lrw;
-  int working_too_hard;
   int sanity_checked;
 
   friend int lsode_f (int *neq, double *t, double *y, double *ydot);
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,14 @@
+Fri Mar 28 15:33:11 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* lsode.cc (struct LSODE_OPTIONS): Handle integer options.
+	(print_lsode_option_list, set_lsode_option, show_lsode_option): Ditto.
+	(lsode_option_table): Add element for step limit.
+	(lsode_user_jacobian): New function.
+	(Flsode): Allow function name arg to be a 2-element string array
+	specifying the function and jacobian function.
+
+	* variables.cc (get_global_value, set_global_value): New functions.
+
 Wed Mar 26 17:08:27 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	Implement static variable declaration:
--- a/src/lsode.cc
+++ b/src/lsode.cc
@@ -43,6 +43,9 @@
 // Global pointer for user defined function required by lsode.
 static tree_fvc *lsode_fcn;
 
+// Global pointer for optional user defined jacobian function used by lsode.
+static tree_fvc *lsode_jac;
+
 static LSODE_options lsode_opts;
 
 ColumnVector
@@ -55,20 +58,11 @@
   octave_value_list args;
   args(1) = t;
 
-  if (nstates > 1)
-    {
-      Matrix m (nstates, 1);
-     for (int i = 0; i < nstates; i++)
-	m (i, 0) = x (i);
-      octave_value state (m);
-      args(0) = state;
-    }
-  else
-    {
-      double d = x (0);
-      octave_value state (d);
-      args(0) = state;
-    }
+  Matrix m (nstates, 1);
+  for (int i = 0; i < nstates; i++)
+    m (i, 0) = x (i);
+  octave_value state (m);
+  args(0) = state;
 
   if (lsode_fcn)
     {
@@ -94,6 +88,46 @@
   return retval;
 }
 
+Matrix
+lsode_user_jacobian (const ColumnVector& x, double t)
+{
+  Matrix retval;
+
+  int nstates = x.capacity ();
+
+  octave_value_list args;
+  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 (lsode_jac)
+    {
+      octave_value_list tmp = lsode_jac->eval (0, 1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("lsode");
+	  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");
+	}
+      else
+	gripe_user_supplied_eval ("lsode");
+    }
+
+  return retval;
+}
+
 DEFUN_DLD (lsode, args, nargout,
   "lsode (F, X0, T_OUT, T_CRIT)\n\
 \n\
@@ -114,12 +148,48 @@
       return retval;
     }
 
-  lsode_fcn = extract_function
-    (args(0), "lsode", "__lsode_fcn__",
-     "function xdot = __lsode_fcn__ (x, t) xdot = ",
-     "; endfunction");
+  octave_value f_arg = args(0);
+
+  switch (f_arg.rows ())
+    {
+    case 1:
+      lsode_fcn = extract_function
+	(args(0), "lsode", "__lsode_fcn__",
+	 "function xdot = __lsode_fcn__ (x, t) xdot = ",
+	 "; endfunction");
+      break;
+
+    case 2:
+      {
+	string_vector tmp = args(0).all_strings ();
 
-  if (! lsode_fcn)
+	if (! error_state)
+	  {
+	    lsode_fcn = extract_function
+	      (tmp(0), "lsode", "__lsode_fcn__",
+	       "function xdot = __lsode_fcn__ (x, t) xdot = ",
+	       "; endfunction");
+
+	    if (lsode_fcn)
+	      {
+		lsode_jac = extract_function
+		  (tmp(1), "lsode", "__lsode_jac__",
+		   "function jac = __lsode_jac__ (x, t) jac = ",
+		   "; endfunction");
+
+		if (! lsode_jac)
+		  lsode_fcn = 0;
+	      }
+	  }
+      }
+      break;
+
+    default:
+      error ("lsode: second arg should be a string or 2-element string array");
+      break;
+    }
+
+  if (error_state || ! lsode_fcn)
     return retval;
 
   ColumnVector state = args(1).vector_value ();
@@ -158,7 +228,11 @@
   int nsteps = out_times.capacity ();
 
   ODEFunc func (lsode_user_function);
+  if (lsode_jac)
+    func.set_jacobian_function (lsode_user_jacobian);
+
   LSODE ode (state, tzero, func);
+
   ode.copy (lsode_opts);
 
   int nstates = state.capacity ();
@@ -179,7 +253,9 @@
 }
 
 typedef void (LSODE_options::*d_set_opt_mf) (double);
+typedef void (LSODE_options::*i_set_opt_mf) (int);
 typedef double (LSODE_options::*d_get_opt_mf) (void);
+typedef int (LSODE_options::*i_get_opt_mf) (void);
 
 #define MAX_TOKENS 3
 
@@ -190,7 +266,9 @@
   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 LSODE_OPTIONS lsode_option_table [] =
@@ -198,37 +276,43 @@
   { "absolute tolerance",
     { "absolute", "tolerance", 0, 0, },
     { 1, 0, 0, 0, }, 1,
-    LSODE_options::set_absolute_tolerance,
-    LSODE_options::absolute_tolerance, },
+    LSODE_options::set_absolute_tolerance, 0,
+    LSODE_options::absolute_tolerance, 0, },
 
   { "initial step size",
     { "initial", "step", "size", 0, },
     { 1, 0, 0, 0, }, 1,
-    LSODE_options::set_initial_step_size,
-    LSODE_options::initial_step_size, },
+    LSODE_options::set_initial_step_size, 0,
+    LSODE_options::initial_step_size, 0, },
 
   { "maximum step size",
     { "maximum", "step", "size", 0, },
     { 2, 0, 0, 0, }, 1,
-    LSODE_options::set_maximum_step_size,
-    LSODE_options::maximum_step_size, },
+    LSODE_options::set_maximum_step_size, 0,
+    LSODE_options::maximum_step_size, 0, },
 
   { "minimum step size",
     { "minimum", "step", "size", 0, },
     { 2, 0, 0, 0, }, 1,
-    LSODE_options::set_minimum_step_size,
-    LSODE_options::minimum_step_size, },
+    LSODE_options::set_minimum_step_size, 0,
+    LSODE_options::minimum_step_size, 0, },
 
   { "relative tolerance",
     { "relative", "tolerance", 0, 0, },
     { 1, 0, 0, 0, }, 1,
-    LSODE_options::set_relative_tolerance,
-    LSODE_options::relative_tolerance, },
+    LSODE_options::set_relative_tolerance, 0,
+    LSODE_options::relative_tolerance, 0, },
+
+  { "step limit",
+    { "step", "limit", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    0, LSODE_options::set_step_limit,
+    0, LSODE_options::step_limit, },
 
   { 0,
     { 0, 0, 0, 0, },
     { 0, 0, 0, 0, }, 0,
-    0, 0, },
+    0, 0, 0, 0, },
 };
 
 static void
@@ -247,13 +331,22 @@
   while ((keyword = list->keyword) != 0)
     {
       os.form ("  %-40s ", keyword);
-
-      double val = (lsode_opts.*list->d_get_fcn) ();
-      if (val < 0.0)
-	os << "computed automatically";
+      if (list->d_get_fcn)
+	{
+	  double val = (lsode_opts.*list->d_get_fcn) ();
+	  if (val < 0.0)
+	    os << "computed automatically";
+	  else
+	    os << val;
+	}
       else
-	os << val;
-
+	{
+	  int val = (lsode_opts.*list->i_get_fcn) ();
+	  if (val < 0)
+	    os << "infinite";
+	  else
+	    os << val;
+	}
       os << "\n";
       list++;
     }
@@ -271,8 +364,18 @@
       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
 				list->min_toks_to_match, MAX_TOKENS))
 	{
-	  (lsode_opts.*list->d_set_fcn) (val);
-
+	  if (list->d_set_fcn)
+	    (lsode_opts.*list->d_set_fcn) (val);
+	  else
+	    {
+	      if (xisnan (val))
+		{
+		  error ("lsode_options: %s: expecting integer, found NaN",
+			 keyword.c_str ());
+		}
+	      else
+		(lsode_opts.*list->i_set_fcn) (NINT (val));
+	    }
 	  return;
 	}
       list++;
@@ -293,11 +396,22 @@
       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
 				list->min_toks_to_match, MAX_TOKENS))
 	{
-	  double val = (lsode_opts.*list->d_get_fcn) ();
-	  if (val < 0.0)
-	    retval = "computed automatically";
+	  if (list->d_get_fcn)
+	    {
+	      double val = (lsode_opts.*list->d_get_fcn) ();
+	      if (val < 0.0)
+		retval = "computed automatically";
+	      else
+		retval = val;
+	    }
 	  else
-	    retval = val;
+	    {
+	      int val = (lsode_opts.*list->i_get_fcn) ();
+	      if (val < 0)
+		retval = "infinite";
+	      else
+		retval = static_cast<double> (val);
+	    }
 
 	  return retval;
 	}