diff liboctave/DASSL.cc @ 4049:a35a3c5d4740

[project @ 2002-08-16 08:54:31 by jwe]
author jwe
date Fri, 16 Aug 2002 08:54:31 +0000
parents 7b0c139ac8af
children b79da8779a0e
line wrap: on
line diff
--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -56,53 +56,7 @@
 
 static int nn;
 
-DASSL::DASSL (void) : DAE ()
-{
-  liw = 0;
-  lrw = 0;
-
-  sanity_checked = false;
-
-  info.resize (15);
-
-  for (int i = 0; i < 15; i++)
-    info.elem (i) = 0;
-}
-
-DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
-  : DAE (state, time, f)
-{
-  n = size ();
-
-  liw = 20 + n;
-  lrw = 40 + 9*n + n*n;
-
-  sanity_checked = false;
-
-  info.resize (15, 0);
-}
-
-DASSL::DASSL (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 ());
-
-  liw = 20 + n;
-  lrw = 40 + 9*n + n*n;
-
-  sanity_checked = false;
-
-  info.resize (15);
-
-  for (int i = 0; i < 15; i++)
-    info.elem (i) = 0;
-}
-
-int
+static int
 ddassl_f (const double& time, const double *state, const double *deriv,
 	  double *delta, int& ires, double *, int *)
 {
@@ -134,7 +88,7 @@
   return 0;
 }
 
-int
+static int
 ddassl_j (const double& time, const double *state, const double *deriv,
 	  double *pd, const double& cj, double *, int *)
 {
@@ -161,140 +115,158 @@
 ColumnVector
 DASSL::do_integrate (double tout)
 {
-  // XXX FIXME XXX -- should handle all this option stuff just once
-  // for each new problem.
-
   ColumnVector retval;
 
-  if (restart)
+  if (! initialized || restart || DAEFunc::reset|| DASSL_options::reset)
     {
-      restart = false;
-      info.elem (0) = 0;
-    }
+      integration_error = false;
+
+      initialized = true;
+
+      info.resize (15);
+
+      for (int i = 0; i < 15; i++)
+	info(i) = 0;
 
-  if (iwork.length () != liw)
-    iwork.resize (liw);
+      pinfo = info.fortran_vec ();
+
+      int n = size ();
+
+      liw = 20 + n;
+      lrw = 40 + 9*n + n*n;
 
-  if (rwork.length () != lrw)
-    rwork.resize (lrw);
+      nn = n;
 
-  integration_error = false;
+      iwork.resize (liw);
+      rwork.resize (lrw);
 
-  user_fun = DAEFunc::fun;
-  user_jac = DAEFunc::jac;
+      info(0) = 0;
 
-  if (user_jac)
-    info.elem (4) = 1;
-  else
-    info.elem (4) = 0;
+      if (stop_time_set)
+	{
+	  rwork(0) = stop_time;
+	  info(3) = 1;
+	}
+      else
+	info(3) = 0;
 
-  double *px = x.fortran_vec ();
-  double *pxdot = xdot.fortran_vec ();
+      px = x.fortran_vec ();
+      pxdot = xdot.fortran_vec ();
+
+      piwork = iwork.fortran_vec ();
+      prwork = rwork.fortran_vec ();
+
+      restart = false;
+
+      // DAEFunc
 
-  nn = n;
+      user_fun = DAEFunc::function ();
+      user_jac = DAEFunc::jacobian_function ();
+
+      if (user_fun)
+	{
+	  int ires = 0;
+
+	  ColumnVector res = (*user_fun) (x, xdot, t, ires);
 
-  if (! sanity_checked)
-    {
-      int ires = 0;
+	  if (res.length () != x.length ())
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: 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)
-	    ("dassl: inconsistent sizes for state and residual vectors");
+	    ("dassl: no user supplied RHS subroutine!");
 
 	  integration_error = true;
 	  return retval;
 	}
 
-      sanity_checked = true;
-    }
+      info(4) = user_jac ? 1 : 0;
   
-  if (stop_time_set)
-    {
-      rwork.elem (0) = stop_time;
-      info.elem (3) = 1;
-    }
-  else
-    info.elem (3) = 0;
+      DAEFunc::reset = false;
 
-  Array<double> abs_tol = absolute_tolerance ();
-  Array<double> rel_tol = relative_tolerance ();
+      // DASSL_options
 
-  int abs_tol_len = abs_tol.length ();
-  int rel_tol_len = rel_tol.length ();
+      double hmax = maximum_step_size ();
+      if (hmax >= 0.0)
+	{
+	  rwork(1) = hmax;
+	  info(6) = 1;
+	}
+      else
+	info(6) = 0;
 
-  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)
-	("dassl: inconsistent sizes for tolerance arrays");
+      double h0 = initial_step_size ();
+      if (h0 >= 0.0)
+	{
+	  rwork(2) = h0;
+	  info(7) = 1;
+	}
+      else
+	info(7) = 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;
+      int maxord = maximum_order ();
+      if (maxord >= 0)
+	{
+	  if (maxord > 0 && maxord < 6)
+	    {
+	      info(8) = 1;
+	      iwork(2) = maxord;
+	    }
+	  else
+	    {
+	      (*current_liboctave_error_handler)
+		("dassl: invalid value for maximum order");
+	      integration_error = true;
+	      return retval;
+	    }
+	}
 
-  double h0 = initial_step_size ();
-  if (h0 >= 0.0)
-    {
-      rwork.elem (2) = h0;
-      info.elem (7) = 1;
-    }
-  else
-    info.elem (7) = 0;
+      int enc = enforce_nonnegativity_constraints ();
+      info(9) = enc ? 1 : 0;
+
+      int ccic = compute_consistent_initial_condition ();
+      info(10) = ccic ? 1 : 0;
+
+      abs_tol = absolute_tolerance ();
+      rel_tol = relative_tolerance ();
 
-  int maxord = maximum_order ();
-  if (maxord >= 0)
-    {
-      if (maxord > 0 && maxord < 6)
+      int abs_tol_len = abs_tol.length ();
+      int rel_tol_len = rel_tol.length ();
+
+      if (abs_tol_len == 1 && rel_tol_len == 1)
 	{
-	  info(8) = 1;
-	  iwork(2) = maxord;
+	  info(1) = 0;
+	}
+      else if (abs_tol_len == n && rel_tol_len == n)
+	{
+	  info(1) = 1;
 	}
       else
 	{
 	  (*current_liboctave_error_handler)
-	    ("dassl: invalid value for maximum order");
+	    ("dassl: inconsistent sizes for tolerance arrays");
+
 	  integration_error = true;
 	  return retval;
 	}
+
+      pabs_tol = abs_tol.fortran_vec ();
+      prel_tol = rel_tol.fortran_vec ();
+
+      DASSL_options::reset = false;
     }
 
-  int enc = enforce_nonnegativity_constraints ();
-  info (9) = enc ? 1 : 0;
-
-  int ccic = compute_consistent_initial_condition ();
-  info(10) = ccic ? 1 : 0;
-
-  double *dummy = 0;
-  int *idummy = 0;
+  static double *dummy = 0;
+  static 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 (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
+  F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo,
 			     prel_tol, pabs_tol, istate, prwork, lrw,
 			     piwork, liw, dummy, idummy, ddassl_j));
 
@@ -366,7 +338,9 @@
 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {
@@ -409,7 +383,9 @@
 		  const ColumnVector& tcrit) 
 {
   Matrix retval;
+
   int n_out = tout.capacity ();
+  int n = size ();
 
   if (n_out > 0 && n > 0)
     {