Mercurial > hg > octave-max

changeset 4140:303b28a7a7e4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-11-01 02:53:13 by jwe]
author jwe
date Fri, 01 Nov 2002 02:53:14 +0000
parents 02ca908056e9
children 8c710385c572
files src/ChangeLog src/DLD-FUNCTIONS/daspk.cc src/DLD-FUNCTIONS/dasrt.cc src/DLD-FUNCTIONS/dassl.cc src/DLD-FUNCTIONS/fsolve.cc src/DLD-FUNCTIONS/lsode.cc src/DLD-FUNCTIONS/odessa.cc src/DLD-FUNCTIONS/quad.cc
diffstat 8 files changed, 350 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,6 +1,22 @@
 2002-10-31  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
-	* DLD-FUNCTIONS/fsolve.cc (fsolve_user_function): Print warning if
+	* DLD-FUNCTIONS/fsolve.cc (fsolve_user_function, fsolve_user_jacobian):
+	Print warning if user returns complex value.
+	* DLD-FUNCTIONS/quad.cc (quad_user_function): Likewise.
+	* DLD-FUNCTIONS/lsode.cc (lsode_user_function, lsode_user_jacobian):
+	Likewise.
+	* DLD-FUNCTIONS/dassl.cc (dassl_user_function, dassl_user_jacobian):
+	Likewise.
+	* DLD-FUNCTIONS/dasrt.cc (dasrt_user_f, dasrt_user_cf, dasrt_user_j):
+	Likewise. 
+	* DLD-FUNCTIONS/daspk.cc (daspk_user_function, daspk_user_jacobian):
+	Likewise.
+
+	* DLD-FUNCTIONS/daspk.cc (daspk_user_jacobian): New function.
+	(Fdaspk): Handle extracting Jacobian from function argument.
+
+	* DLD-FUNCTIONS/fsolve.cc (fsolve_user_function): New function.
+	(Ffsolve): Handle extracting Jacobian from function argument.
 
 	* Makefile.in (%.oct): Depend on octave$(EXEEXT) so that octave
 	will be built before any .oct files.
--- a/src/DLD-FUNCTIONS/daspk.cc
+++ b/src/DLD-FUNCTIONS/daspk.cc
@@ -46,6 +46,13 @@
 // Global pointer for user defined function required by daspk.
 static octave_function *daspk_fcn;
 
+// Global pointer for optional user defined jacobian function.
+static octave_function *daspk_jac;
+
+// Have we warned about imaginary values returned from user function?
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -99,6 +106,12 @@
       int tlen = tmp.length ();
       if (tlen > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("daspk: ignoring imaginary part returned from user-supplied function");
+	      warned_fcn_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (tlen > 1)
@@ -114,6 +127,76 @@
   return retval;
 }
 
+Matrix
+daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot,
+		     double t, double cj)
+{
+  Matrix retval;
+
+  int nstates = x.capacity ();
+
+  assert (nstates == xdot.capacity ());
+
+  octave_value_list args;
+
+  args(3) = cj;
+  args(2) = t;
+
+  if (nstates > 1)
+    {
+      Matrix m1 (nstates, 1);
+      Matrix m2 (nstates, 1);
+      for (int i = 0; i < nstates; i++)
+	{
+	  m1 (i, 0) = x (i);
+	  m2 (i, 0) = xdot (i);
+	}
+      octave_value state (m1);
+      octave_value deriv (m2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+  else
+    {
+      double d1 = x (0);
+      double d2 = xdot (0);
+      octave_value state (d1);
+      octave_value deriv (d2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+
+  if (daspk_jac)
+    {
+      octave_value_list tmp = daspk_jac->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("daspk");
+	  return retval;
+	}
+
+      int tlen = tmp.length ();
+      if (tlen > 0 && tmp(0).is_defined ())
+	{
+	  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_jac_imaginary = true;
+	    }
+
+	  retval = tmp(0).matrix_value ();
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("daspk");
+	}
+      else
+	gripe_user_supplied_eval ("daspk");
+    }
+
+  return retval;
+}
+
 #define DASPK_ABORT() \
   do \
     { \
@@ -235,6 +318,9 @@
 {
   octave_value_list retval;
 
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
+
   unwind_protect::begin_frame ("Fdaspk");
 
   unwind_protect_int (call_depth);
@@ -247,12 +333,46 @@
 
   if (nargin > 3 && nargin < 6)
     {
-      daspk_fcn = extract_function
-	(args(0), "daspk", "__daspk_fcn__",
-	 "function res = __daspk_fcn__ (x, xdot, t) res = ",
-	 "; endfunction");
+      daspk_fcn = 0;
+      daspk_jac = 0;
+
+      octave_value f_arg = args(0);
+
+      switch (f_arg.rows ())
+	{
+	case 1:
+	  daspk_fcn = extract_function
+	    (args(0), "daspk", "__daspk_fcn__",
+	     "function res = __daspk_fcn__ (x, xdot, t) res = ",
+	     "; endfunction");
+	  break;
+
+	case 2:
+	  {
+	    string_vector tmp = f_arg.all_strings ();
 
-      if (! daspk_fcn)
+	    if (! error_state)
+	      {
+		daspk_fcn = extract_function
+		  (tmp(0), "daspk", "__daspk_fcn__",
+		   "function res = __daspk_fcn__ (x, xdot, t) res = ",
+		   "; endfunction");
+
+		if (daspk_fcn)
+		  {
+		    daspk_jac = extract_function
+		      (tmp(1), "daspk", "__daspk_jac__",
+		       "function jac = __daspk_jac__ (x, xdot, t, cj) jac = ",
+		       "; endfunction");
+
+		    if (! daspk_jac)
+		      daspk_fcn = 0;
+		  }
+	      }
+	  }
+	}
+
+      if (error_state || ! daspk_fcn)
 	DASPK_ABORT ();
 
       ColumnVector state = ColumnVector (args(1).vector_value ());
@@ -288,6 +408,9 @@
       double tzero = out_times (0);
 
       DAEFunc func (daspk_user_function);
+      if (daspk_jac)
+	func.set_jacobian_function (daspk_user_jacobian);
+
       DASPK dae (state, deriv, tzero, func);
       dae.set_options (daspk_opts);
 
--- a/src/DLD-FUNCTIONS/dasrt.cc
+++ b/src/DLD-FUNCTIONS/dasrt.cc
@@ -48,6 +48,11 @@
 static octave_function *dasrt_j;
 static octave_function *dasrt_cf;
 
+// Have we warned about imaginary values returned from user function?
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
+static bool warned_cf_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -91,6 +96,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("dasrt: ignoring imaginary part returned from user-supplied function");
+	      warned_fcn_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (error_state || retval.length () == 0)
@@ -133,6 +144,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_cf_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("dasrt: ignoring imaginary part returned from user-supplied constraint function");
+	      warned_cf_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (error_state || retval.length () == 0)
@@ -197,6 +214,12 @@
       int tlen = tmp.length ();
       if (tlen > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_jac_imaginary = true;
+	    }
+
 	  retval = tmp(0).matrix_value ();
 
 	  if (error_state || retval.length () == 0)
@@ -366,6 +389,10 @@
 {
   octave_value_list retval;
 
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
+  warned_cf_imaginary = false;
+
   unwind_protect::begin_frame ("Fdasrt");
 
   unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/dassl.cc
+++ b/src/DLD-FUNCTIONS/dassl.cc
@@ -49,6 +49,10 @@
 // Global pointer for optional user defined jacobian function.
 static octave_function *dassl_jac;
 
+// Have we warned about imaginary values returned from user function?
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -102,6 +106,12 @@
       int tlen = tmp.length ();
       if (tlen > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("dassl: ignoring imaginary part returned from user-supplied function");
+	      warned_fcn_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (tlen > 1)
@@ -169,6 +179,12 @@
       int tlen = tmp.length ();
       if (tlen > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_jac_imaginary = true;
+	    }
+
 	  retval = tmp(0).matrix_value ();
 
 	  if (error_state || retval.length () == 0)
@@ -302,6 +318,9 @@
 {
   octave_value_list retval;
 
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
+
   unwind_protect::begin_frame ("Fdassl");
 
   unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/fsolve.cc
+++ b/src/DLD-FUNCTIONS/fsolve.cc
@@ -46,8 +46,12 @@
 // Global pointer for user defined function required by hybrd1.
 static octave_function *fsolve_fcn;
 
+// Global pointer for optional user defined jacobian function.
+static octave_function *fsolve_jac;
+
 // Have we warned about imaginary values returned from user function?
-static bool warned_imaginary = false;
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
 
 // Is this a recursive call?
 static int call_depth = 0;
@@ -91,7 +95,7 @@
 {
   ColumnVector retval;
 
-  int n = x.capacity ();
+  int n = x.length ();
 
   octave_value_list args;
   args.resize (1);
@@ -117,10 +121,10 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
-	  if (! warned_imaginary && tmp(0).is_complex_type ())
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
 	    {
 	      warning ("fsolve: ignoring imaginary part returned from user-supplied function");
-	      warned_imaginary = true;
+	      warned_fcn_imaginary = true;
 	    }
 
 	  retval = ColumnVector (tmp(0).vector_value ());
@@ -135,6 +139,55 @@
   return retval;
 }
 
+Matrix
+fsolve_user_jacobian (const ColumnVector& x)
+{
+  Matrix retval;
+
+  int n = x.length ();
+
+  octave_value_list args;
+  args.resize (1);
+
+  if (n > 1)
+    {
+      Matrix m (n, 1);
+      for (int i = 0; i < n; i++)
+	m (i, 0) = x (i);
+      octave_value vars (m);
+      args(0) = vars;
+    }
+  else
+    {
+      double d = x (0);
+      octave_value vars (d);
+      args(0) = vars;
+    }
+
+  if (fsolve_fcn)
+    {
+      octave_value_list tmp = fsolve_jac->do_multi_index_op (1, args);
+
+      if (tmp.length () > 0 && tmp(0).is_defined ())
+	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_fcn_imaginary = true;
+	    }
+
+	  retval = tmp(0).matrix_value ();
+
+	  if (error_state || retval.length () <= 0)
+	    gripe_user_supplied_eval ("fsolve");
+	}
+      else
+	gripe_user_supplied_eval ("fsolve");
+    }
+
+  return retval;
+}
+
 #define FSOLVE_ABORT() \
   do \
     { \
@@ -172,7 +225,8 @@
 {
   octave_value_list retval;
 
-  warned_imaginary = false;
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
 
   unwind_protect::begin_frame ("Ffsolve");
 
@@ -186,10 +240,46 @@
 
   if (nargin == 2 && nargout < 4)
     {
-      fsolve_fcn = extract_function (args(0), "fsolve", "__fsolve_fcn__",
-				    "function y = __fsolve_fcn__ (x) y = ",
-				    "; endfunction");
-      if (! fsolve_fcn)
+      fsolve_fcn = 0;
+      fsolve_jac = 0;
+
+      octave_value f_arg = args(0);
+
+      switch (f_arg.rows ())
+	{
+	case 1:
+	  fsolve_fcn = extract_function
+	    (f_arg, "fsolve", "__fsolve_fcn__",
+	     "function y = __fsolve_fcn__ (x) y = ",
+	     "; endfunction");
+	  break;
+
+	case 2:
+	  {
+	    string_vector tmp = f_arg.all_strings ();
+
+	    if (! error_state)
+	      {
+		fsolve_fcn = extract_function
+		  (tmp(0), "fsolve", "__fsolve_fcn__",
+		   "function y = __fsolve_fcn__ (x) y = ",
+		   "; endfunction");
+
+		if (fsolve_fcn)
+		  {
+		    fsolve_jac = extract_function
+		      (tmp(1), "fsolve", "__fsolve_jac__",
+		       "function jac = __fsolve_jac__ (x) jac = ",
+		       "; endfunction");
+
+		    if (! fsolve_jac)
+		      fsolve_fcn = 0;
+		  }
+	      }
+	  }
+	}
+
+      if (error_state || ! fsolve_fcn)
 	FSOLVE_ABORT ();
 
       ColumnVector x (args(1).vector_value ());
@@ -204,6 +294,9 @@
 	warning ("fsolve: can't compute path output yet");
 
       NLFunc nleqn_fcn (fsolve_user_function);
+      if (fsolve_jac)
+	nleqn_fcn.set_jacobian_function (fsolve_user_jacobian);
+
       NLEqn nleqn (x, nleqn_fcn);
       nleqn.set_options (fsolve_opts);
 
--- a/src/DLD-FUNCTIONS/lsode.cc
+++ b/src/DLD-FUNCTIONS/lsode.cc
@@ -51,6 +51,10 @@
 // Global pointer for optional user defined jacobian function used by lsode.
 static octave_function *lsode_jac;
 
+// Have we warned about imaginary values returned from user function?
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -82,6 +86,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("lsode: ignoring imaginary part returned from user-supplied function");
+	      warned_fcn_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (error_state || retval.length () == 0)
@@ -122,6 +132,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_jac_imaginary = true;
+	    }
+
 	  retval = tmp(0).matrix_value ();
 
 	  if (error_state || retval.length () == 0)
@@ -268,6 +284,9 @@
 {
   octave_value_list retval;
 
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
+
   unwind_protect::begin_frame ("Flsode");
 
   unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/odessa.cc
+++ b/src/DLD-FUNCTIONS/odessa.cc
@@ -51,6 +51,11 @@
 static octave_function *odessa_j;
 static octave_function *odessa_b;
 
+// Have we warned about imaginary values returned from user function?
+static bool warned_fcn_imaginary = false;
+static bool warned_jac_imaginary = false;
+static bool warned_b_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -92,6 +97,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("odessa: ignoring imaginary part returned from user-supplied function");
+	      warned_fcn_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (error_state || retval.length () == 0)
@@ -143,6 +154,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("odessa: ignoring imaginary part returned from user-supplied jacobian function");
+	      warned_jac_imaginary = true;
+	    }
+
 	  retval = tmp(0).matrix_value ();
 
 	  if (error_state || retval.length () == 0)
@@ -202,6 +219,12 @@
 
       if (tmp.length () > 0 && tmp(0).is_defined ())
 	{
+	  if (! warned_b_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("odessa: ignoring imaginary part returned from user-supplied inhomogeneity function");
+	      warned_b_imaginary = true;
+	    }
+
 	  retval = ColumnVector (tmp(0).vector_value ());
 
 	  if (error_state || retval.length () == 0)
@@ -376,6 +399,10 @@
 {
   octave_value_list retval;
 
+  warned_fcn_imaginary = false;
+  warned_jac_imaginary = false;
+  warned_b_imaginary = false;
+
   unwind_protect::begin_frame ("Fodessa");
 
   unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/quad.cc
+++ b/src/DLD-FUNCTIONS/quad.cc
@@ -51,6 +51,9 @@
 // Global pointer for user defined function required by quadrature functions.
 static octave_function *quad_fcn;
 
+// Have we warned about imaginary values returned from user function?
+static bool warned_imaginary = false;
+
 // Is this a recursive call?
 static int call_depth = 0;
 
@@ -75,6 +78,12 @@
 
       if (tmp.length () && tmp(0).is_defined ())
 	{
+	  if (! warned_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("quad: ignoring imaginary part returned from user-supplied function");
+	      warned_imaginary = true;
+	    }
+
 	  retval = tmp(0).double_value ();
 
 	  if (error_state)
@@ -156,6 +165,8 @@
 {
   octave_value_list retval;
 
+  warned_imaginary = false;
+
   unwind_protect::begin_frame ("Fquad");
 
   unwind_protect_int (call_depth);