scipy/scipy — pull request #3218 by WarrenWeckesser (2014-01-15 23:00:45+00:00)
Title: BUG: integrate: Fix banded jacobian handling in the "vode" and "zvode" solvers.

The "vode" solver in the ode class is a wrapper of the Fortran routine
DVODE. DVODE accepts a JAC argument, which is the user's implementation
of the Jacobian function. JAC has the signature

SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD, RPAR, IPAR)
DOUBLE PRECISION T, Y, PD, RPAR
DIMENSION Y(NEQ), PD(NROWPD,NEQ)

In the banded case, JAC is called with NROWPD = 2*ML + MU + 1.
JAC must fill in PD as described in the DVODE documentation:

In the band matrix case (MITER = 4), the elements
within the band are to be loaded into PD in columnwise
manner, with diagonal lines of df/dy loaded into the rows
of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j).
ML and MU are the half-bandwidth parameters.

Even though PD has dimension (NROWPD, NEQ), and NROWPD is
2*ML + MU + 1, the number of rows in PD to be filled in by JAC
is ML + MU + 1. So there are ML extra rows of PD that can be
ignored in JAC. This is fine in Fortran--it is call-by-reference,
so JAC just fills in the first ML + MU + 1 rows and returns.

The .pyf file used by F2PY to create a wrapper for the user-defined
python Jacobian function specifies the interface for JAC as

subroutine jac(n,t,y,ml,mu,jac,nrowpd,rpar,ipar)
integer intent(hide) :: n
double precision :: t
double precision dimension(n),intent(c,in) :: y
integer intent(hide) :: ml,mu
integer intent(hide):: nrowpd
double precision intent(out) :: jac(nrowpd, n)
double precision intent(hide) :: rpar
integer intent(hide) :: ipar
end subroutine jac

Note that the argument `jac` is declared intent(out), and it has
dimensions (nrowpd, n) (which matches the Fortran signature).
This means that even though there are only ml + mu + 1 nonzero
diagonals (so one would expect to create an array with shape
(ml + mu + 1, n)), the user's Python function must return an array
with shape (2*ml + mu + 1, n). Otherwise, the F2Py wrapper generates
an error.

This is a pointless requirement, since the Python wrapper takes the
returned array and copies it into the Fortran function's JAC argument
So the Python code has to pad its array with rows that are not used.

This commit fixes the problem by making a few small changes to the .pyf
interface file. The declaration of the argument `jac` is changed to

double precision intent(out) :: jac(ipar, n)

and in the declaration of the interface to the dvode function,
the declaration of `ipar` is changed from

integer intent(hide) :: ipar = 0

to
integer intent(hide) :: ipar = (((mf == 14) || (mf == 24)) ? (iwork[0] + iwork[1] + 1) : neq)

The banded case is indicated by the digit 4 in `mf`. `iwork[0]` and
`iwork[1]` hold `ml` and `mu`, respectively. So we are using the
previously unused integer parameter `ipar` to hold the *actual*
required number of rows in the Jacobian matrix.

The "zvode" solver has the same problem. Both "vode" and "zvode"
are fixed in this commit.

The problem was pointed out by @bjodah in https://github.com/scipy/scipy/pull/3206