--- a/liboctave/DASPK.cc	Fri Aug 16 02:43:14 2002 +0000
+++ b/liboctave/DASPK.cc	Fri Aug 16 08:54:31 2002 +0000
@@ -64,47 +64,7 @@
 static DAEFunc::DAEJacFunc user_jac;
 static int nn;
 
-DASPK::DASPK (void) : DAE ()
-{
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-DASPK::DASPK (const ColumnVector& state, double time, DAEFunc& f)
-  : DAE (state, time, f)
-{
-  n = size ();
-
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-DASPK::DASPK (const ColumnVector& state, const ColumnVector& deriv,
-	  double time, DAEFunc& f)
-  : DAE (state, deriv, time, f)
-{
-  n = size ();
-
-  DAEFunc::set_function (f.function ());
-  DAEFunc::set_jacobian_function (f.jacobian_function ());
-
-  sanity_checked = 0;
-
-  info.resize (20);
-
-  for (int i = 0; i < 20; i++)
-    info.elem (i) = 0;
-}
-
-int
+static int
 ddaspk_f (const double& time, const double *state, const double *deriv,
 	  const double&, double *delta, int& ires, double *, int *)
 {
@@ -137,7 +97,7 @@
 //NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
 //C                          WP, IWP, B, EPLIN, IER, RPAR, IPAR)
 
-int
+static int
 ddaspk_psol (const int& neq, const double& time, const double *state,
 	     const double *deriv, const double *savr,
 	     const double& cj, const double *wght, double *wp,
@@ -149,7 +109,7 @@
 }
 
 
-int
+static int
 ddaspk_j (const double& time, const double *state, const double *deriv,
 	  double *pd, const double& cj, double *, int *)
 {
@@ -181,178 +141,220 @@
 
   ColumnVector retval;
 
-  if (restart)
+  if (! initialized || restart || DAEFunc::reset|| DASPK_options::reset)
     {
-      restart = false;
-      info.elem (0) = 0;
-    }
+      integration_error = false;
+
+      initialized = true;
+
+      info.resize (20);
 
-  liw = 40 + n;
-  if (info(9) == 1 || info(9) == 3)
-    liw += n;
-  if (info (10) == 1 || info(15) == 1)
-    liw += n;
+      for (int i = 0; i < 20; i++)
+	info(i) = 0;
+
+      pinfo = info.fortran_vec ();
 
-  lrw = 50 + 9*n;
-  if (info(5) == 0)
-    lrw += n*n;
-  if (info(15) == 1)
-    lrw += n;
+      int n = size ();
 
-  if (iwork.length () != liw)
-    iwork.resize (liw);
+      nn = n;
+
+      info(0) = 0;
 
-  if (rwork.length () != lrw)
-    rwork.resize (lrw);
-
-  integration_error = false;
+      if (stop_time_set)
+	{
+	  rwork(0) = stop_time;
+	  info(3) = 1;
+	}
+      else
+	info(3) = 0;
 
-  if (DAEFunc::jacobian_function ())
-    info.elem (4) = 1;
-  else
-    info.elem (4) = 0;
+      px = x.fortran_vec ();
+      pxdot = xdot.fortran_vec ();
+
+      // DAEFunc
+
+      user_fun = DAEFunc::function ();
+      user_jac = DAEFunc::jacobian_function ();
 
-  double *px    = x.fortran_vec ();
-  double *pxdot = xdot.fortran_vec ();
+      if (user_fun)
+	{
+	  int ires = 0;
 
-  nn = n;
-  user_fun = DAEFunc::fun;
-  user_jac = DAEFunc::jac;
+	  ColumnVector res = (*user_fun) (x, xdot, t, ires);
 
-  if (! sanity_checked)
-    {
-      int ires = 0;
+	  if (res.length () != x.length ())
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: inconsistent sizes for state and residual vectors");
 
-      ColumnVector res = (*user_fun) (x, xdot, t, ires);
-
-      if (res.length () != x.length ())
+	      integration_error = true;
+	      return retval;
+	    }
+	}
+      else
 	{
 	  (*current_liboctave_error_handler)
-	    ("daspk: inconsistent sizes for state and residual vectors");
+	    ("daspk: no user supplied RHS subroutine!");
+
+	  integration_error = true;
+	  return retval;
+	}
+  
+      info(4) = user_jac ? 1 : 0;
+
+      DAEFunc::reset = false;
+
+      // DASPK_options
+
+      Array<double> abs_tol = absolute_tolerance ();
+      Array<double> rel_tol = relative_tolerance ();
+
+      int abs_tol_len = abs_tol.length ();
+      int rel_tol_len = rel_tol.length ();
+
+      if (abs_tol_len == 1 && rel_tol_len == 1)
+	{
+	  info(1) = 0;
+	}
+      else if (abs_tol_len == n && rel_tol_len == n)
+	{
+	  info(1) = 1;
+	}
+      else
+	{
+	  (*current_liboctave_error_handler)
+	    ("daspk: inconsistent sizes for tolerance arrays");
 
 	  integration_error = true;
 	  return retval;
 	}
 
-      sanity_checked = 1;
-    }
-  
-  if (stop_time_set)
-    {
-      rwork.elem (0) = stop_time;
-      info.elem (3) = 1;
-    }
-  else
-    info.elem (3) = 0;
-
-  Array<double> abs_tol = absolute_tolerance ();
-  Array<double> rel_tol = relative_tolerance ();
-
-  int abs_tol_len = abs_tol.length ();
-  int rel_tol_len = rel_tol.length ();
-
-  if (abs_tol_len == 1 && rel_tol_len == 1)
-    {
-      info.elem (1) = 0;
-    }
-  else if (abs_tol_len == n && rel_tol_len == n)
-    {
-      info.elem (1) = 1;
-    }
-  else
-    {
-      (*current_liboctave_error_handler)
-	("daspk: inconsistent sizes for tolerance arrays");
+      double hmax = maximum_step_size ();
+      if (hmax >= 0.0)
+	{
+	  rwork(1) = hmax;
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
 
-      integration_error = true;
-      return retval;
-    }
-
-  double hmax = maximum_step_size ();
-  if (hmax >= 0.0)
-    {
-      rwork.elem (1) = hmax;
-      info.elem (6) = 1;
-    }
-  else
-    info.elem (6) = 0;
-
-  double h0 = initial_step_size ();
-  if (h0 >= 0.0)
-    {
-      rwork.elem (2) = h0;
-      info.elem (7) = 1;
-    }
-  else
-    info.elem (7) = 0;
-
-  int maxord = maximum_order ();
-  if (maxord >= 0)
-    {
-      if (maxord > 0 && maxord < 6)
+      double h0 = initial_step_size ();
+      if (h0 >= 0.0)
 	{
-	  info(8) = 1;
-	  iwork(2) = maxord;
+	  rwork(2) = h0;
+	  info(7) = 1;
 	}
       else
+	info(7) = 0;
+
+      int maxord = maximum_order ();
+      if (maxord >= 0)
 	{
+	  if (maxord > 0 && maxord < 6)
+	    {
+	      info(8) = 1;
+	      iwork(2) = maxord;
+	    }
+	  else
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: invalid value for maximum order");
+	      integration_error = true;
+	      return retval;
+	    }
+	}
+
+      int eiq = enforce_inequality_constraints ();
+      switch (eiq)
+	{
+	case 1:
+	case 3:
+	  {
+	    Array<int> ict = inequality_constraint_types ();
+
+	    if (ict.length () == n)
+	      {
+		for (int i = 0; i < n; i++)
+		  {
+		    int val = ict(i);
+		    if (val < -2 || val > 2)
+		      {
+			(*current_liboctave_error_handler)
+			  ("daspk: invalid value for inequality constraint type");
+			integration_error = true;
+			return retval;
+		      }
+		    iwork(40+i) = val;
+		  }
+	      }
+	    else
+	      {
+		(*current_liboctave_error_handler)
+		  ("daspk: inequality constraint types size mismatch");
+		integration_error = true;
+		return retval;
+	      }
+	  }
+	  // Fall through...
+
+	case 2:
+	  info(9) = eiq;
+	  break;
+
+	default:
 	  (*current_liboctave_error_handler)
-	    ("daspk: invalid value for maximum order");
+	    ("daspk: invalid value for enforce inequality constraints option");
 	  integration_error = true;
 	  return retval;
 	}
-    }
+
+      int ccic = compute_consistent_initial_condition ();
+      if (ccic)
+	{
+	  if (ccic == 1)
+	    {
+	      // XXX FIXME XXX -- this code is duplicated below.
+
+	      Array<int> av = algebraic_variables ();
 
-  int eiq = enforce_inequality_constraints ();
-  switch (eiq)
-    {
-    case 1:
-    case 3:
-      {
-	Array<int> ict = inequality_constraint_types ();
+	      if (av.length () == n)
+		{
+		  int lid;
+		  if (eiq == 0 || eiq == 2)
+		    lid = 40;
+		  else if (eiq == 1 || eiq == 3)
+		    lid = 40 + n;
+		  else
+		    abort ();
 
-	if (ict.length () == n)
-	  {
-	    for (int i = 0; i < n; i++)
-	      {
-		int val = ict(i);
-		if (val < -2 || val > 2)
-		  {
-		    (*current_liboctave_error_handler)
-		      ("daspk: invalid value for inequality constraint type");
-		    integration_error = true;
-		    return retval;
-		  }
-		iwork(40+i) = val;
-	      }
-	  }
-	else
-	  {
-	    (*current_liboctave_error_handler)
-	      ("daspk: inequality constraint types size mismatch");
-	    integration_error = true;
-	    return retval;
-	  }
-      }
-      // Fall through...
+		  for (int i = 0; i < n; i++)
+		    iwork(lid+i) = av(i) ? -1 : 1;
+		}
+	      else
+		{
+		  (*current_liboctave_error_handler)
+		    ("daspk: algebraic variables size mismatch");
+		  integration_error = true;
+		  return retval;
+		}
+	    }
+	  else if (ccic != 2)
+	    {
+	      (*current_liboctave_error_handler)
+		("daspk: invalid value for compute consistent initial condition option");
+	      integration_error = true;
+	      return retval;
+	    }
 
-    case 2:
-      info(9) = eiq;
-      break;
+	  info(10) = ccic;
+	}
 
-    default:
-      (*current_liboctave_error_handler)
-	("daspk: invalid value for enforce inequality constraints option");
-      integration_error = true;
-      return retval;
-    }
+      int eavfet = exclude_algebraic_variables_from_error_test ();
+      if (eavfet)
+	{
+	  info(15) = 1;
 
-  int ccic = compute_consistent_initial_condition ();
-  if (ccic)
-    {
-      if (ccic == 1)
-	{
-	  // XXX FIXME XXX -- this code is duplicated below.
+	  // XXX FIXME XXX -- this code is duplicated above.
 
 	  Array<int> av = algebraic_variables ();
 
@@ -369,102 +371,77 @@
 	      for (int i = 0; i < n; i++)
 		iwork(lid+i) = av(i) ? -1 : 1;
 	    }
+	}
+
+      if (use_initial_condition_heuristics ())
+	{
+	  Array<double> ich = initial_condition_heuristics ();
+
+	  if (ich.length () == 6)
+	    {
+	      iwork(31) = NINT (ich(0));
+	      iwork(32) = NINT (ich(1));
+	      iwork(33) = NINT (ich(2));
+	      iwork(34) = NINT (ich(3));
+
+	      rwork(13) = ich(4);
+	      rwork(14) = ich(5);
+	    }
 	  else
 	    {
 	      (*current_liboctave_error_handler)
-		("daspk: algebraic variables size mismatch");
+		("daspk: invalid initial condition heuristics option");
 	      integration_error = true;
 	      return retval;
 	    }
+
+	  info(16) = 1;
 	}
-      else if (ccic != 2)
+
+      int pici = print_initial_condition_info ();
+      switch (pici)
 	{
+	case 0:
+	case 1:
+	case 2:
+	  info(17) = pici;
+	  break;
+
+	default:
 	  (*current_liboctave_error_handler)
-	    ("daspk: invalid value for compute consistent initial condition option");
+	    ("daspk: invalid value for print initial condition info option");
 	  integration_error = true;
 	  return retval;
+	  break;
 	}
 
-      info(10) = ccic;
-    }
+      DASPK_options::reset = false;
 
-  if (exclude_algebraic_variables_from_error_test ())
-    {
-      info(15) = 1;
-
-      // XXX FIXME XXX -- this code is duplicated above.
-
-      Array<int> av = algebraic_variables ();
+      liw = 40 + n;
+      if (eiq == 1 || eiq == 3)
+	liw += n;
+      if (ccic == 1 || eavfet == 1)
+	liw += n;
 
-      if (av.length () == n)
-	{
-	  int lid;
-	  if (eiq == 0 || eiq == 2)
-	    lid = 40;
-	  else if (eiq == 1 || eiq == 3)
-	    lid = 40 + n;
-	  else
-	    abort ();
+      lrw = 50 + 9*n;
+      if (! user_jac)
+	lrw += n*n;
+      if (eavfet == 1)
+	lrw += n;
 
-	  for (int i = 0; i < n; i++)
-	    iwork(lid+i) = av(i) ? -1 : 1;
-	}
+      iwork.resize (liw);
+      rwork.resize (lrw);
+
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
     }
 
-  if (use_initial_condition_heuristics ())
-    {
-      Array<double> ich = initial_condition_heuristics ();
-
-      if (ich.length () == 6)
-	{
-	  iwork(31) = NINT (ich(0));
-	  iwork(32) = NINT (ich(1));
-	  iwork(33) = NINT (ich(2));
-	  iwork(34) = NINT (ich(3));
-
-	  rwork(13) = ich(4);
-	  rwork(14) = ich(5);
-	}
-      else
-	{
-	  (*current_liboctave_error_handler)
-	    ("daspk: invalid initial condition heuristics option");
-	  integration_error = true;
-	  return retval;
-	}
-
-      info(16) = 1;
-    }
+  static double *dummy = 0;
+  static int *idummy = 0;
 
-  int pici = print_initial_condition_info ();
-  switch (pici)
-    {
-    case 0:
-    case 1:
-    case 2:
-      info(17) = pici;
-      break;
-
-    default:
-      (*current_liboctave_error_handler)
-	("daspk: invalid value for print initial condition info option");
-      integration_error = true;
-      return retval;
-      break;
-    }
-
-  double *dummy = 0;
-  int *idummy = 0;
-
-  int *pinfo = info.fortran_vec ();
-  int *piwork = iwork.fortran_vec ();
-  double *prwork = rwork.fortran_vec ();
-  double *pabs_tol = abs_tol.fortran_vec ();
-  double *prel_tol = rel_tol.fortran_vec ();
-
-// again:
-
-  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, n, t, px, pxdot, tout, pinfo,
+  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo,
 			     prel_tol, pabs_tol, istate, prwork, lrw,
 			     piwork, liw, dummy, idummy, ddaspk_j,
 			     ddaspk_psol));
@@ -546,7 +523,9 @@
 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -589,7 +568,9 @@
 		  const ColumnVector& tcrit) 
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {