diff liboctave/LSODE.cc @ 258:1c9a678906fb

[project @ 1993-12-14 06:24:04 by jwe]
author jwe
date Tue, 14 Dec 1993 06:24:21 +0000
parents 780cbbc57b7c
children c23f50e61c58
line wrap: on
line diff
--- a/liboctave/LSODE.cc	Wed Dec 08 23:36:43 1993 +0000
+++ b/liboctave/LSODE.cc	Tue Dec 14 06:24:21 1993 +0000
@@ -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;