Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Analytical banded Jacobian with CVODE_BDF #302

Open
Nicholaswogan opened this issue Mar 28, 2021 · 3 comments
Open

Analytical banded Jacobian with CVODE_BDF #302

Nicholaswogan opened this issue Mar 28, 2021 · 3 comments

Comments

@Nicholaswogan
Copy link

I am trying to use an analytical banded jacobian with CVODE_BDF. My goal is to provide CVODE_BDF with a jacobian function which returns the jacobian in band storage (not a dense matrix). It doesn't seem like this feature is available based on my post on Julialang discourse

@JianghuiDu
Copy link

JianghuiDu commented Apr 15, 2021

Would be great to have! My current workaround is to used a sparse jacobian type with :KLU despite the jacobian being banded, which does improve performance. And I hope using a banded jacobian type with :Band it will be even better!

@JianghuiDu
Copy link

JianghuiDu commented Oct 12, 2021

I think a function converting BandedMatrix to SUNMatrix is missing. There's already one for Matrix

function Base.convert(::Type{Matrix}, J::SUNMatrix)
    _sunmat = unsafe_load(J)
    _mat = convert(SUNMatrixContent_Dense, _sunmat.content)
    mat = unsafe_load(_mat)
    # own is false as memory is allocated by sundials
    unsafe_wrap(Array, mat.data, (mat.M, mat.N), own = false)
end

This conversion is needed here:

function cvodejac(
    t::realtype,
    u::N_Vector,
    du::N_Vector,
    J::SUNMatrix,
    funjac::AbstractFunJac{Nothing},
    tmp1::N_Vector,
    tmp2::N_Vector,
    tmp3::N_Vector,
)

    funjac.u = unsafe_wrap(Vector{Float64}, N_VGetArrayPointer_Serial(u), length(funjac.u))
    _u = funjac.u
    funjac.jac(convert(Matrix, J), _u, funjac.p, t)
    return CV_SUCCESS
end

I don't know enough to write Base.convert(::Type{BandedMatrix}, J::SUNMatrix)

@ChrisRackauckas
Copy link
Member

You'd need to dig into the Sundials implementation of the banded matrix and figure out how to map that over to the BandedMatrices.jl version. I can't say I've looked too deeply at it either. It could be as easy as defining the Arrays directly from the pointers, or it might take an inversion or something more.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants