# HG changeset patch
# User jwe
# Date 1029444775 0
# Node ID 7b0c139ac8af78a93e8155b63a9dd250e958def3
# Parent  7787c144d5d9ccc2d79c3adcec49c63a79532495
[project @ 2002-08-15 20:52:55 by jwe]

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,8 @@
+2002-08-15  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DASSL.cc (DASSL::do_integrate (double)): Handle more optoins.
+	* DASPK.cc (DASPK::do_integrate (double)): Likewise.
+
 2002-08-15  Paul Kienzle <pkienzle@users.sf.net>
 
 	* DASPK-opts.in, DASPK.h: Move include to .in file.
diff --git a/liboctave/DASPK-opts.in b/liboctave/DASPK-opts.in
--- a/liboctave/DASPK-opts.in
+++ b/liboctave/DASPK-opts.in
@@ -50,6 +50,26 @@
 END_OPTION
 
 OPTION
+  NAME = "algebraic variables"
+  TYPE = "Array<int>"
+  SET_ARG_TYPE = const $TYPE&
+  INIT_BODY
+    $OPTVAR.resize (1);
+    $OPTVAR(0) = 0;
+  END_INIT_BODY
+  SET_CODE
+    void set_$OPT (int val)
+      {
+        $OPTVAR.resize (1);
+        $OPTVAR(0) = val;
+      }
+
+    void set_$OPT (const $TYPE& val)
+      { $OPTVAR = val; }
+  END_SET_CODE
+END_OPTION
+
+OPTION
   NAME = "enforce inequality constraints"
   TYPE = "int"
   INIT_VALUE = "0"
@@ -68,7 +88,7 @@
     void set_$OPT (int val)
       {
         $OPTVAR.resize (1);
-        $OPTVAR(0) = (val > 0.0) ? val : 0;
+        $OPTVAR(0) = val;
       }
 
     void set_$OPT (const $TYPE& val)
@@ -77,55 +97,36 @@
 END_OPTION
 
 OPTION
-  NAME = "exclude algebraic variables in error test"
+  NAME = "exclude algebraic variables from error test"
+  TYPE = "int"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
+  NAME = "use initial condition heuristics"
   TYPE = "int"
   INIT_VALUE = "0"
   SET_EXPR = "val"
 END_OPTION
 
 OPTION
-  NAME = "initial condition maximum step"
-  TYPE = "double"
-  INIT_VALUE = "-1.0"
-  SET_EXPR = "(val >= 0.0) ? val : -1.0"
-END_OPTION
-
-OPTION
-  NAME = "initial condition maximum jacobian evaluations"
-  TYPE = "int"
-  INIT_VALUE = "-1"
-  SET_EXPR = "(val >= 0) ? val : -1"
-END_OPTION
-
-OPTION
-  NAME = "initial condition maximum newton iterations"
-  TYPE = "int"
-  INIT_VALUE = "-1"
-  SET_EXPR = "(val >= 0) ? val : -1"
-END_OPTION
-
-OPTION
-  NAME = "initial condition minimum linesearch step"
-  TYPE = "double"
-  INIT_VALUE = "-1.0"
-  SET_EXPR = "(val >= 0.0) ? val : 0.0"
-END_OPTION
-
-OPTION
-  NAME = "initial condition omit linesearch"
-  TYPE = "int"
-  INIT_VALUE = "0"
+  NAME = "initial condition heuristics"
+  TYPE = "Array<double>"
+  SET_ARG_TYPE = "const $TYPE&"
+  INIT_BODY
+    $OPTVAR.resize (6, 0.0);
+    $OPTVAR(0) = 5.0;
+    $OPTVAR(1) = 6.0;
+    $OPTVAR(2) = 5.0;
+    $OPTVAR(3) = 0.0;
+    $OPTVAR(4) = ::pow (DBL_EPSILON, 2.0/3.0);
+    $OPTVAR(5) = 0.01;
+  END_INIT_BODY
   SET_EXPR = "val"
 END_OPTION
 
 OPTION
-  NAME = "initial condition swing factor"
-  TYPE = "double"
-  INIT_VALUE = "1.0"
-  SET_EXPR = "(val >= 0.0) ? val : -1.0"
-END_OPTION
-
-OPTION
   NAME = "initial step size"
   TYPE = "double"
   INIT_VALUE = "-1.0"
@@ -147,13 +148,6 @@
 END_OPTION
 
 OPTION
-  NAME = "minimum step size"
-  TYPE = "double"
-  INIT_VALUE = "0.0"
-  SET_EXPR = "(val >= 0.0) ? val : 0.0"
-END_OPTION
-
-OPTION
   NAME = "print initial condition info"
   TYPE = "int"
   INIT_VALUE = "0"
diff --git a/liboctave/DASPK.cc b/liboctave/DASPK.cc
--- a/liboctave/DASPK.cc
+++ b/liboctave/DASPK.cc
@@ -68,9 +68,9 @@
 {
   sanity_checked = 0;
 
-  info.resize (15);
+  info.resize (20);
 
-  for (int i = 0; i < 15; i++)
+  for (int i = 0; i < 20; i++)
     info.elem (i) = 0;
 }
 
@@ -176,6 +176,9 @@
 ColumnVector
 DASPK::do_integrate (double tout)
 {
+  // XXX FIXME XXX -- should handle all this option stuff just once
+  // for each new problem.
+
   ColumnVector retval;
 
   if (restart)
@@ -259,27 +262,196 @@
   else
     {
       (*current_liboctave_error_handler)
-	("dassl: inconsistent sizes for tolerance arrays");
+	("daspk: inconsistent sizes for tolerance arrays");
 
       integration_error = true;
       return retval;
     }
 
-  if (initial_step_size () >= 0.0)
+  double hmax = maximum_step_size ();
+  if (hmax >= 0.0)
     {
-      rwork.elem (2) = initial_step_size ();
+      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;
 
-  if (maximum_step_size () >= 0.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 enforce inequality constraints option");
+      integration_error = true;
+      return retval;
+    }
+
+  int ccic = compute_consistent_initial_condition ();
+  if (ccic)
     {
-      rwork.elem (1) = maximum_step_size ();
-      info.elem (6) = 1;
+      if (ccic == 1)
+	{
+	  // XXX FIXME XXX -- this code is duplicated below.
+
+	  Array<int> av = algebraic_variables ();
+
+	  if (av.length () == n)
+	    {
+	      int lid;
+	      if (eiq == 0 || eiq == 2)
+		lid = 40;
+	      else if (eiq == 1 || eiq == 3)
+		lid = 40 + n;
+	      else
+		abort ();
+
+	      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;
+	}
+
+      info(10) = ccic;
     }
-  else
-    info.elem (6) = 0;
+
+  if (exclude_algebraic_variables_from_error_test ())
+    {
+      info(15) = 1;
+
+      // XXX FIXME XXX -- this code is duplicated above.
+
+      Array<int> av = algebraic_variables ();
+
+      if (av.length () == n)
+	{
+	  int lid;
+	  if (eiq == 0 || eiq == 2)
+	    lid = 40;
+	  else if (eiq == 1 || eiq == 3)
+	    lid = 40 + n;
+	  else
+	    abort ();
+
+	  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: invalid initial condition heuristics option");
+	  integration_error = true;
+	  return retval;
+	}
+
+      info(16) = 1;
+    }
+
+  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;
diff --git a/liboctave/DASSL-opts.in b/liboctave/DASSL-opts.in
--- a/liboctave/DASSL-opts.in
+++ b/liboctave/DASSL-opts.in
@@ -43,6 +43,20 @@
 END_OPTION
 
 OPTION
+  NAME = "compute consistent initial condition"
+  TYPE = "int"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
+  NAME = "enforce nonnegativity constraints"
+  TYPE = "int"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
   NAME = "initial step size"
   TYPE = "double"
   INIT_VALUE = "-1.0"
@@ -62,10 +76,3 @@
   INIT_VALUE = "-1.0"
   SET_EXPR = "(val >= 0.0) ? val : -1.0"
 END_OPTION
-
-OPTION
-  NAME = "minimum step size"
-  TYPE = "double"
-  INIT_VALUE = "0.0"
-  SET_EXPR = "(val >= 0.0) ? val : 0.0"
-END_OPTION
diff --git a/liboctave/DASSL.cc b/liboctave/DASSL.cc
--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -161,6 +161,9 @@
 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)
@@ -239,21 +242,46 @@
       return retval;
     }
 
-  if (initial_step_size () >= 0.0)
+  double hmax = maximum_step_size ();
+  if (hmax >= 0.0)
     {
-      rwork.elem (2) = initial_step_size ();
+      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;
 
-  if (maximum_step_size () >= 0.0)
+  int maxord = maximum_order ();
+  if (maxord >= 0)
     {
-      rwork.elem (1) = maximum_step_size ();
-      info.elem (6) = 1;
+      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;
+	}
     }
-  else
-    info.elem (6) = 0;
+
+  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;
diff --git a/mk-opts.pl b/mk-opts.pl
--- a/mk-opts.pl
+++ b/mk-opts.pl
@@ -566,6 +566,8 @@
 {
   local ($i);
 
+  ## XXX FIXME XXX -- determine the width of the table automatically.
+
   print "static void
 print_${class_name} (std::ostream& os)
 {
@@ -573,15 +575,15 @@
 
   os << \"\\n\"
      << \"Options for $CLASS include:\\n\\n\"
-     << \"  keyword                                   value\\n\"
-     << \"  -------                                   -----\\n\";
+     << \"  keyword                                             value\\n\"
+     << \"  -------                                             -----\\n\";
 
   $struct_name *list = $static_table_name;\n\n";
 
   for ($i = 0; $i < $opt_num; $i++)
     {
       print "  {\n    os << \"  \"
-       << std::setiosflags (std::ios::left) << std::setw (40)
+       << std::setiosflags (std::ios::left) << std::setw (50)
        << list[$i].keyword
        << std::resetiosflags (std::ios::left)
        << \"  \";\n\n";