changeset 258:1c9a678906fb

[project @ 1993-12-14 06:24:04 by jwe]
author jwe
date Tue, 14 Dec 1993 06:24:21 +0000
parents 126791334c68
children 23866011a5f2
files liboctave/LSODE.cc liboctave/ODE.h src/lsode.cc
diffstat 3 files changed, 48 insertions(+), 11 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/LSODE.cc
+++ b/liboctave/LSODE.cc
@@ -31,9 +31,12 @@
 
 extern "C"
 {
-  int F77_FCN (lsode) (int (*)(), int *, double *, double *, double *,
+  int F77_FCN (lsode) (int (*)(int*, double*, double*, double*, int*),
+		       int *, double *, double *, double *,
 		       int *, double *, double *, int *, int *, int *,
-		       double *, int *, int *, int *, int (*)(), int *);
+		       double *, int *, int *, int *,
+		       int (*)(int*, double*, double*, int*, int*,
+			       double*, int*), int *);
 }
 
 static ColumnVector (*user_fun) (ColumnVector&, double);
@@ -51,6 +54,7 @@
   stop_time_set = 0;
   stop_time = 0.0;
 
+  integration_error = 0;
   restart = 1;
 
   istate = 1;
@@ -84,6 +88,7 @@
   stop_time_set = 0;
   stop_time = 0.0;
 
+  integration_error = 0;
   restart = 1;
 
   istate = 1;
@@ -118,6 +123,7 @@
   stop_time_set = 0;
   stop_time = 0.0;
 
+  integration_error = 0;
   restart = 1;
 
   istate = 1;
@@ -147,7 +153,7 @@
 }
 
 int
-lsode_f (int *neq, double *time, double *state, double *deriv)
+lsode_f (int *neq, double *time, double *state, double *deriv, int *ierr)
 {
   int nn = *neq;
   ColumnVector tmp_deriv (nn);
@@ -159,8 +165,13 @@
    */
   tmp_deriv = (*user_fun) (*tmp_x, *time);
 
-  for (int i = 0; i < nn; i++)
-    deriv [i] = tmp_deriv.elem (i);
+  if (tmp_deriv.length () == 0)
+    *ierr = -1;
+  else
+    {
+      for (int i = 0; i < nn; i++)
+	deriv [i] = tmp_deriv.elem (i);
+    }
 
   return 0;
 }
@@ -194,6 +205,8 @@
   else
     method_flag = 21;
 
+  integration_error = 0;
+
   double *xp = x.fortran_vec ();
 
 // NOTE: this won't work if LSODE passes copies of the state vector.
@@ -236,17 +249,16 @@
 
   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.)
-      break;
     case -5: // repeated convergence failures (perhaps bad jacobian
 	     // supplied or wrong choice of mf or tolerances).
-      break;
     case -4: // repeated error test failures (check all inputs).
-      break;
     case -3: // illegal input detected (see printed message).
-      break;
     case -2: // excess accuracy requested (tolerances too small).
+      integration_error = 1;
+      return ColumnVector ();
       break;
     case -1: // excess work done on this call (perhaps wrong mf).
       working_too_hard++;
@@ -255,6 +267,7 @@
 	  (*current_liboctave_error_handler)
 	    ("giving up after more than %d steps attempted in lsode",
 	     iwork[5] * 20);
+	  integration_error = 1;
 	  return ColumnVector ();
 	}
       else
@@ -317,6 +330,10 @@
       for (int j = 1; j < n_out; j++)
 	{
 	  ColumnVector x_next = integrate (tout.elem (j));
+
+	  if (integration_error)
+	    return retval;
+
 	  for (i = 0; i < n; i++)
 	    retval.elem (j, i) = x_next.elem (i);
 	}
@@ -394,6 +411,9 @@
 
 	      ColumnVector x_next = integrate (t_out);
 
+	      if (integration_error)
+		return retval;
+
 	      if (save_output)
 		{
 		  for (i = 0; i < n; i++)
@@ -405,7 +425,12 @@
 	    }
 	}
       else
-	retval = integrate (tout);
+	{
+	  retval = integrate (tout);
+
+	  if (integration_error)
+	    return retval;
+	}
     }
 
   return retval;
--- a/liboctave/ODE.h
+++ b/liboctave/ODE.h
@@ -76,6 +76,7 @@
 
 private:
 
+  int integration_error;
   int restart;
   int method_flag;
   int *iwork;
--- a/src/lsode.cc
+++ b/src/lsode.cc
@@ -75,17 +75,28 @@
   if (lsode_fcn != NULL_TREE)
     {
       tree_constant *tmp = lsode_fcn->eval (args, 3, 1, 0);
+
       delete [] args;
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("lsode");
+	  return retval;
+	}
+
       if (tmp != NULL_TREE_CONST && tmp[0].is_defined ())
 	{
 	  retval = tmp[0].to_vector ();
+
 	  delete [] tmp;
+
+	  if (retval.length () == 0)
+	    gripe_user_supplied_eval ("lsode");
 	}
       else
 	{
 	  delete [] tmp;
 	  gripe_user_supplied_eval ("lsode");
-	  jump_to_top_level ();
 	}
     }