For banded Jacobian matrices, we need a new Jacobian procedure interface:
!> User supplied subprogram for evaluation of the Jacobian.
subroutine jacobian_sub(n,y,df)
import wp
integer, intent(in) :: n
real(wp), intent(in) :: y(n)
real(wp), intent(inout) :: df(n,n)
end subroutine
There are a few possible choices. One is to create an interface that allows both dense and banded matrices (from LSODA):
C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
The constraint NROWPD >= (ML + MU + 1) is expected to hold. STIFF3 would require an extra option to toggle between dense and banded.
We could also have an exact interface:
subroutine jac(neqn, y, ml, mu, df)
integer, intent(in) :: neqn, ml, mu
real(wp), intent(in) :: y(neqn)
real(wp), intent(out) :: df(-ml:mu,neqn)
end subroutine
Using custom bounds makes the indexing convenient, because the diagonal is placed on row 0. To accomodate the fill-in expected also by the LAPACK banded factorization routines, we could even do:
subroutine jac(neqn, y, ml, mu, df)
integer, intent(in) :: neqn, ml, mu
real(wp), intent(in) :: y(neqn)
real(wp), intent(out) :: df(-ml:*,neqn)
end subroutine
This makes indexing easy, but we lose the ability for strict bounds-checking. In practice we know that atleast ml extra rows are needed, so the upper bound will likely be mu+ml, but it may be higher if we want to use padding for optimization.
For banded Jacobian matrices, we need a new Jacobian procedure interface:
There are a few possible choices. One is to create an interface that allows both dense and banded matrices (from LSODA):
The constraint NROWPD >= (ML + MU + 1) is expected to hold. STIFF3 would require an extra option to toggle between dense and banded.
We could also have an exact interface:
Using custom bounds makes the indexing convenient, because the diagonal is placed on row 0. To accomodate the fill-in expected also by the LAPACK banded factorization routines, we could even do:
This makes indexing easy, but we lose the ability for strict bounds-checking. In practice we know that atleast ml extra rows are needed, so the upper bound will likely be mu+ml, but it may be higher if we want to use padding for optimization.