BUG: integrate: Fix banded jacobian handling in the "vode" and "zvode" solvers.

[WarrenWeckesser]

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.

**That's a summary of the problem. The following describes the changes in this PR, but those changes won't work, so you can stop reading. :) **

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 #3206

[coveralls]

Changes Unknown when pulling 56b4f7d on WarrenWeckesser:vode-bjac into * on scipy:master*.

[pv]

I have some trouble understanding why this would work.

jac points to a Fortran-order array, with strides (sizeof(double), nrowpd*sizeof(double)).
Specifying jac(ipar, n) makes f2py interpret the memory pointer as an array with strides (sizeof(double), ipar*sizeof(double)).
Do the values returned by the callback then end up in the correct locations?

This would probably be OK, if the array was in C-order, but in Fortran order it seems dangerous...

[WarrenWeckesser]

@pv: Thanks. You're correct, it doesn't work.

Back to the drawing board...