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

Using FBDF solver to deal with super-large-sparse-matrix problems #906

Closed
ytdHuang opened this issue Oct 11, 2022 · 4 comments
Closed

Using FBDF solver to deal with super-large-sparse-matrix problems #906

ytdHuang opened this issue Oct 11, 2022 · 4 comments

Comments

@ytdHuang
Copy link

ytdHuang commented Oct 11, 2022

Hi,

I have some issues when solving ODE with super-large (but sparse) matrix
The ODE problem I want to solve can be considered as:
$\frac{d\rho(t)}{dt}=L\rho(t)$,
where $\rho$ is a Vector{ComplexF64} and $L$ is a SparseMatrixCSC{ComplexF64, Int64} with size: 9,523,200 × 9,523,200 and 92,553,046 stored entries.

Therefore, I tried to solve this problem with solver FBDF :

using OrdinaryDiffEq
 
function func!(dρ, ρ, L, t)
    @inbounds dρ .= L * ρ
end

ρ = ...  # generate the initial vector
L = ...  # generate the large sparse matrix
tspan = (0, 1000)

sol = solve(
        ODEProblem(func!, ρ, tspan, L), 
        FBDF(autodiff=false);
        abstol = 1e-8, 
        reltol = 1e-6,
        maxiters = 1e5,
        save_everystep = false
    )

And here comes a whole bunch of errors:

OutOfMemoryError()

Stacktrace:
  [1] Array
    @ ./boot.jl:461 [inlined]
  [2] zeromatrix
    @ /usr/local/julia-1.8.0/share/julia/packages/ArrayInterfaceCore/j22dF/src/ArrayInterfaceCore.jl:511 [inlined]
  [3] build_J_W(alg::FBDF{5, 12, false, LinearSolve.LUFactorization{RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, u::Vector{ComplexF64}, uprev::Vector{ComplexF64}, p::SparseMatrixCSC{ComplexF64, Int64}, t::Float64, dt::Float64, f::ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, #unused#::Type{ComplexF64}, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/derivative_utils.jl:889
  [4] build_nlsolver(alg::FBDF{5, 12, false, LinearSolve.LUFactorization{RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, nlalg::NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, u::Vector{ComplexF64}, uprev::Vector{ComplexF64}, p::SparseMatrixCSC{ComplexF64, Int64}, t::Float64, dt::Float64, f::Function, rate_prototype::Vector{ComplexF64}, #unused#::Type{ComplexF64}, #unused#::Type{ComplexF64}, #unused#::Type{Float64}, γ::Float64, c::Float64, α::Int64, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/nlsolve/utils.jl:174
  [5] build_nlsolver
    @ ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/nlsolve/utils.jl:143 [inlined]
  [6] build_nlsolver
    @ ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/nlsolve/utils.jl:134 [inlined]
  [7] alg_cache(alg::FBDF{5, 12, false, LinearSolve.LUFactorization{RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, u::Vector{ComplexF64}, rate_prototype::Vector{ComplexF64}, #unused#::Type{ComplexF64}, #unused#::Type{ComplexF64}, #unused#::Type{Float64}, uprev::Vector{ComplexF64}, uprev2::Vector{ComplexF64}, f::ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::SparseMatrixCSC{ComplexF64, Int64}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/caches/bdf_caches.jl:582
  [8] __init(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::FBDF{5, 12, false, LinearSolve.LUFactorization{RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::DiscreteCallback{DiffEqCallbacks.var"#55#57"{Nothing, Float64, Float64, typeof(DiffEqCallbacks.allDerivPass)}, DiffEqCallbacks.var"#56#58", typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Rational{Int64}, qsteady_max::Rational{Int64}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/solve.jl:316
  [9] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/vfMzV/src/solve.jl:5 [inlined]
 [10] #solve_call#26
    @ /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:472 [inlined]
 [11] solve_up(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{ComplexF64}, p::SparseMatrixCSC{ComplexF64, Int64}, args::FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:maxiters, :callback, :save_everystep, :save_start), Tuple{Float64, DiscreteCallback{DiffEqCallbacks.var"#55#57"{Nothing, Float64, Float64, typeof(DiffEqCallbacks.allDerivPass)}, DiffEqCallbacks.var"#56#58", typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, Bool, Bool}}})
    @ DiffEqBase /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:834
 [12] solve(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:maxiters, :callback, :save_everystep, :save_start), Tuple{Float64, DiscreteCallback{DiffEqCallbacks.var"#55#57"{Nothing, Float64, Float64, typeof(DiffEqCallbacks.allDerivPass)}, DiffEqCallbacks.var"#56#58", typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, Bool, Bool}}})
    @ DiffEqBase /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:801
 [13] __solve(::SteadyStateProblem{Vector{ComplexF64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::SteadyStateDiffEq.DynamicSS{FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, Float64, Float64, Float64}; save_everystep::Bool, save_start::Bool, save_idxs::Nothing, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:maxiters,), Tuple{Float64}}})
    @ SteadyStateDiffEq /usr/local/julia-1.8.0/share/julia/packages/SteadyStateDiffEq/6hOEc/src/solve.jl:81
 [14] solve_call(_prob::SteadyStateProblem{Vector{ComplexF64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::SteadyStateDiffEq.DynamicSS{FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, Float64, Float64, Float64}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:maxiters, :save_everystep), Tuple{Float64, Bool}}})
    @ DiffEqBase /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:472
 [15] solve_up(prob::SteadyStateProblem{Vector{ComplexF64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Vector{ComplexF64}, p::SparseMatrixCSC{ComplexF64, Int64}, args::SteadyStateDiffEq.DynamicSS{FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, Float64, Float64, Float64}; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:maxiters, :save_everystep), Tuple{Float64, Bool}}})
    @ DiffEqBase /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:834
 [16] solve(prob::SteadyStateProblem{Vector{ComplexF64}, true, SparseMatrixCSC{ComplexF64, Int64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Heom.HeomAPI._hierarchy!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::SteadyStateDiffEq.DynamicSS{FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, Float64, Float64, Float64}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:maxiters, :save_everystep), Tuple{Float64, Bool}}})
    @ DiffEqBase /usr/local/julia-1.8.0/share/julia/packages/DiffEqBase/5rKYk/src/solve.jl:801
 [17] SteadyState(M::M_Fermion, ados::ADOs; solver::FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, reltol::Float64, abstol::Float64, maxiters::Float64, save_everystep::Bool, verbose::Bool, SOLVEROptions::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Heom.HeomAPI ~/heom/Heom.jl/src/ADOs/SteadyState.jl:144
 [18] SteadyState(M::M_Fermion, ρ0::Matrix{Int64}; solver::FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, reltol::Float64, abstol::Float64, maxiters::Float64, save_everystep::Bool, verbose::Bool, SOLVEROptions::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Heom.HeomAPI ~/heom/Heom.jl/src/ADOs/SteadyState.jl:90
 [19] top-level scope
    @ ./timing.jl:263 [inlined]
 [20] top-level scope
    @ ./In[21]:0
 [21] eval
    @ ./boot.jl:368 [inlined]
 [22] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1428

I know that I'm working in a very large matrix, but the work station I'm using has totally "3TB" of memory.
The memory should be enough if I'm working under sparse matrix (?)
I've also tried another solver DP5() and it works fine (just need to take a couple days to get the solution).

In Stacktrace [2], I found the keyword zeromatrix.
I think, during the solving process, FBDF tries to construct a zero matrix which has the same size as my input matrix L but in standard matrix type Matrix{ComplexF64}?
If then, is it possible for FBDF to contruct the zero matrix in type: SparseMatrixCSC ?

The reason I'm using other solvers like FBDF is because that I follow the suggestions in the documentation, and try to see if it can speed up the calculations I'm working on.
Or does anyone have other suggested solvers for my problem? I have tried CVODE_BDF in package Sundials.jl at the beginning but it doesn't support for complex numbers.

@ChrisRackauckas
Copy link
Member

@ytdHuang
Copy link
Author

ytdHuang commented Oct 12, 2022

@ChrisRackauckas , thanks for the advice.
I have added the jac_prototype=SparseMatrixCSC{ComplexF64, Int64} and create an ODEFunction.
But I found out that the FBDF() solver is taking much more time than DP5() (maybe not even solving), which was not expected after reading your documentations.
Maybe I should show my code in detail:
Because I also want to handle something in each time step (from a given tlist) during the solving process, I choosed to use the integrator interface:

using SparseArrays
using OrdinaryDiffEq
 
function func!(dρ, ρ, L, t)
    @inbounds dρ .= L * ρ
end

ρ::Vector{ComplexF64}  # a vector with size: 9,523,200.
L::SparseMatrixCSC{ComplexF64, Int64}  # a large sparse matrix with size: 9,523,200 × 9,523,200 and 92,553,046 stored entries.
tlist = 0:1:1000

# setup ode function and solver
dim, = size(L)
ode = ODEFunction(func!; jac_prototype = spzeros(ComplexF64, dim, dim))
solver = FBDF(autodiff=false)

# set up integrator
integrator = init(
    ODEProblem(ode, ρ, (tlist[1], tlist[end]), L),
    solver;
    reltol = 1.0e-6,
    abstol = 1.0e-8,
    maxiters = 1e5,
    save_everystep = false
)

dt_list = diff(tlist)
for dt in dt_list
    step!(integrator, dt, true)

    # Do some other things
    # basicly like printing some message in each iteration 
    # and store the solution to a file at each time step
    # since this is a long calculation (spend couple days)
end

I also tried to change the solver: solver = FBDF(autodiff = false, linsolve = UMFPACKFactorization()),
where UMFPACKFactorization() is given by LinearSolve.jl.
The reason of choosing UMFPACK is that I also need to calculate L \ ρ in another analysis, and UMFPACK work much more faster than KLU in this case.

Maybe this is not the best solver for me to use in this case.
Does anyone have any advice? I have been struggling in this problem for two months and really need some help.

@ytdHuang ytdHuang changed the title OutOfMermory error when using FBDF solver Using FBDF solver to deal with super-large-sparse-matrix jacobian problems Oct 12, 2022
@ChrisRackauckas
Copy link
Member

Is the ODE stiff? If it's non-stiff, then the explicit methods will just be faster. That would be dependent on the eigenvalues of L, if they have values which are orders of magnitude apart.

Note that for linear ODEs, using LinearExponential will be the fastest.

@ytdHuang ytdHuang changed the title Using FBDF solver to deal with super-large-sparse-matrix jacobian problems Using FBDF solver to deal with super-large-sparse-matrix problems Oct 13, 2022
@ytdHuang
Copy link
Author

@ChrisRackauckas thanks for the advice again.

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

2 participants