# HG changeset patch # User Olaf Till # Date 1289548385 -3600 # Fri Nov 12 08:53:05 2010 +0100 # Node ID 122ca0f6df7337d6b519fb1db3ee97b86ce00006 # Parent 7797e2f3c69b14550b3aff45e3de87dae86694be 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 7797e2f3c69b -r 122ca0f6df73 liboctave/numeric/LSODE-opts.in --- a/liboctave/numeric/LSODE-opts.in Tue Feb 23 09:29:31 2021 -0800 +++ b/liboctave/numeric/LSODE-opts.in Fri Nov 12 08:53:05 2010 +0100 @@ -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 @qcode{"full"} +The default. All partial derivatives are approximated or used from the +user-supplied Jacobian function. + +@item @qcode{"banded"} +Only the diagonal and the number of lower and upper subdiagonals specified by +the options @qcode{"lower jacobian subdiagonals"} and @qcode{"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 @qcode{"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 @qcode{"jacobian type"} is set to +@qcode{"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 @qcode{"jacobian type"} is set to +@qcode{"banded"}. The default is zero. + + END_DOC_ITEM + TYPE = "octave_idx_type" + INIT_VALUE = "0" + SET_EXPR = "val" +END_OPTION diff -r 7797e2f3c69b -r 122ca0f6df73 liboctave/numeric/LSODE.cc --- a/liboctave/numeric/LSODE.cc Tue Feb 23 09:29:31 2021 -0800 +++ b/liboctave/numeric/LSODE.cc Fri Nov 12 08:53:05 2010 +0100 @@ -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);