# HG changeset patch
# User Lachlan Andrew <lachlanbis@gmail.com>
# Date 1583865959 14400
#      Tue Mar 10 14:45:59 2020 -0400
# Node ID 2d6e25b9428c591fa349c0eb19ee3e7a3c21114c
# Parent  38a9f6444eb015714ca13ef043e657e403bb0c24
Make Jacobian types offered by dlsode.f accessible by lsode (bug #31626)

* LSODE-opts.in: Add options "jacobian type", "lower jacobian
  subdiagonals", "upper jacobian subdiagonals".
* LSODE.cc (file scope, lsode_j, LSODE::do_integrate (double)):
  Handle new configurable Jacobian types.

diff -r 38a9f6444eb0 -r 2d6e25b9428c liboctave/numeric/LSODE-opts.in
--- a/liboctave/numeric/LSODE-opts.in	Tue Mar 10 07:51:55 2020 -0400
+++ b/liboctave/numeric/LSODE-opts.in	Tue Mar 10 14:45:59 2020 -0400
@@ -164,3 +164,67 @@
   INIT_VALUE = "100000"
   SET_EXPR = "val"
 END_OPTION
+
+OPTION
+  NAME = "jacobian type"
+  DOC_ITEM
+A string specifying the type of Jacobian used with the stiff backward
+differentiation formula (BDF) integration method. Valid values are
+
+@table @asis
+@item \"full\"
+The default. All partial derivatives are approximated or used from the
+user-supplied Jacobian function.
+
+@item \"banded\"
+Only the diagonal and the number of lower and upper subdiagonals specified by
+the options \"lower jacobian subdiagonals\" and \"upper jacobian
+subdiagonals\", respectively, are approximated or used from the user-supplied
+Jacobian function. A user-supplied Jacobian function may set all other
+partial derivatives to arbitrary values.
+
+@item \"diagonal\"
+If a Jacobian function is supplied by the user, this setting has no effect.
+A Jacobian approximated by @code{lsode} is restricted to the diagonal, where
+each partial derivative is computed by applying a finite change to all
+elements of the state together; if the real Jacobian is indeed always diagonal,
+this has the same effect as applying the finite change only to the respective
+element of the state, but is more efficient.
+@end table
+
+  END_DOC_ITEM
+  TYPE = "std::string"
+  SET_ARG_TYPE = "const $TYPE&"
+  INIT_VALUE = ""full""
+  SET_BODY
+    if (val == "full" || val == "banded" || val == "diagonal")
+      $OPTVAR = val;
+    else
+      (*current_liboctave_error_handler)
+        ("lsode_options: jacobian type must be \"full\", \"banded\", or \"diagonal\"");
+  END_SET_BODY
+END_OPTION
+
+OPTION
+  NAME = "lower jacobian subdiagonals"
+  DOC_ITEM
+Number of lower subdiagonals used if option \"jacobian type\" is set to
+\"banded\". The default is zero.
+
+  END_DOC_ITEM
+  TYPE = "octave_idx_type"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
+  NAME = "upper jacobian subdiagonals"
+  DOC_ITEM
+Number of upper subdiagonals used if option \"jacobian type\" is set to
+\"banded\". The default is zero.
+
+  END_DOC_ITEM
+  TYPE = "octave_idx_type"
+  INIT_VALUE = "0"
+  SET_EXPR = "val"
+END_OPTION
diff -r 38a9f6444eb0 -r 2d6e25b9428c liboctave/numeric/LSODE.cc
--- a/liboctave/numeric/LSODE.cc	Tue Mar 10 07:51:55 2020 -0400
+++ b/liboctave/numeric/LSODE.cc	Tue Mar 10 14:45:59 2020 -0400
@@ -55,6 +55,7 @@
 static ODEFunc::ODERHSFunc user_fun;
 static ODEFunc::ODEJacFunc user_jac;
 static ColumnVector *tmp_x;
+static bool user_jac_ignore_ml_mu;
 
 static F77_INT
 lsode_f (const F77_INT& neq, const double& time, double *, double *deriv,
@@ -80,8 +81,9 @@
 }
 
 static F77_INT
-lsode_j (const F77_INT& neq, const double& time, double *, const F77_INT&,
-         const F77_INT&, double *pd, const F77_INT& nrowpd)
+lsode_j (const F77_INT& neq, const double& time, double *, 
+         const F77_INT& ml, const F77_INT& mu,
+         double *pd, const F77_INT& nrowpd)
 {
   Matrix tmp_jac (neq, neq);
 
@@ -91,10 +93,17 @@
 
   tmp_jac = (*user_jac) (*tmp_x, time);
 
-  for (F77_INT j = 0; j < neq; j++)
-    for (F77_INT i = 0; i < neq; i++)
-      pd[nrowpd * j + i] = tmp_jac (i, j);
-
+  if (user_jac_ignore_ml_mu)
+    for (F77_INT j = 0; j < neq; j++)
+      for (F77_INT i = 0; i < neq; i++)
+        pd[nrowpd * j + i] = tmp_jac (i, j);
+  else
+    // upper left ends of subdiagonals in tmp_jac
+    for (F77_INT i = 0, j = mu; i <= ml; j == 0 ? i++ : j--)
+      // walk down the subdiagonal in tmp_jac
+      for (F77_INT k = i, l = j; k < neq && l < neq; k++, l++)
+        pd[nrowpd * l + k + mu - l] = tmp_jac (k, l);    
+      
   return 0;
 }
 
@@ -119,14 +128,43 @@
 
       octave_idx_type max_maxord = 0;
 
+      user_jac_ignore_ml_mu = true;
+
+      iwork(0) = lower_jacobian_subdiagonals (); // 'ML' in dlsode.f
+
+      iwork(1) = upper_jacobian_subdiagonals (); // 'MU' in dlsode.f
+
       if (integration_method () == "stiff")
         {
           max_maxord = 5;
 
           if (jac)
-            method_flag = 21;
+            {
+              if (jacobian_type () == "banded")
+                {
+                  method_flag = 24;
+                  user_jac_ignore_ml_mu = false;
+                }
+              else
+                method_flag = 21;
+            }
           else
-            method_flag = 22;
+            {
+              if (jacobian_type () == "full")
+                method_flag = 22;
+              else if (jacobian_type () == "banded")
+                method_flag = 25;
+              else if (jacobian_type () == "diagonal")
+                method_flag = 23;
+              else
+                {
+                  // should be prevented by lsode_options
+                  (*current_liboctave_error_handler)
+                    ("lsode: internal error, wrong jacobian type");
+                  integration_error = true;
+                  return retval;
+                }
+            }
 
           liw = 20 + n;
           lrw = 22 + n * (9 + n);
