# HG changeset patch # User Lachlan Andrew # Date 1448701174 -39600 # Sat Nov 28 19:59:34 2015 +1100 # Node ID e9938ad656aa0d3f3c28baf77105524976c580f7 # Parent ed708a7ce0a248d10cc4cbaee95506944b8201f9 Make Jacobian types offered by dlsode.f accessible by lsode (patch by Olaf Til) * 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 ed708a7ce0a2 -r e9938ad656aa liboctave/numeric/LSODE-opts.in --- a/liboctave/numeric/LSODE-opts.in Fri Nov 27 16:00:30 2015 -0500 +++ b/liboctave/numeric/LSODE-opts.in Sat Nov 28 19:59:34 2015 +1100 @@ -157,3 +157,67 @@ Maximum number of steps allowed (default 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 ed708a7ce0a2 -r e9938ad656aa liboctave/numeric/LSODE.cc --- a/liboctave/numeric/LSODE.cc Fri Nov 27 16:00:30 2015 -0500 +++ b/liboctave/numeric/LSODE.cc Sat Nov 28 19:59:34 2015 +1100 @@ -59,6 +59,7 @@ extern "C" static ODEFunc::ODERHSFunc user_fun; static ODEFunc::ODEJacFunc user_jac; static ColumnVector *tmp_x; +static bool user_jac_ignore_ml_mu; static octave_idx_type lsode_f (const octave_idx_type& neq, const double& time, double *, @@ -89,8 +90,8 @@ lsode_f (const octave_idx_type& neq, con static octave_idx_type lsode_j (const octave_idx_type& neq, const double& time, double *, - const octave_idx_type&, const octave_idx_type&, double *pd, - const octave_idx_type& nrowpd) + const octave_idx_type& ml, const octave_idx_type& mu, + double *pd, const octave_idx_type& nrowpd) { BEGIN_INTERRUPT_WITH_EXCEPTIONS; @@ -102,9 +103,16 @@ lsode_j (const octave_idx_type& neq, con tmp_jac = (*user_jac) (*tmp_x, time); - for (octave_idx_type j = 0; j < neq; j++) - for (octave_idx_type i = 0; i < neq; i++) - pd[nrowpd * j + i] = tmp_jac (i, j); + if (user_jac_ignore_ml_mu) + for (octave_idx_type j = 0; j < neq; j++) + for (octave_idx_type i = 0; i < neq; i++) + pd[nrowpd * j + i] = tmp_jac (i, j); + else + // upper left ends of subdiagonals in tmp_jac + for (octave_idx_type i = 0, j = mu; i <= ml; j == 0 ? i++ : j--) + // walk down the subdiagonal in tmp_jac + for (octave_idx_type k = i, l = j; k < neq && l < neq; k++, l++) + pd[nrowpd * l + k + mu - l] = tmp_jac (k, l); END_INTERRUPT_WITH_EXCEPTIONS; @@ -132,14 +140,43 @@ LSODE::do_integrate (double tout) 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);