Duplicating integrator with KLUFactorization(; reuse_symbolic = true)

General Usage
, ,
1

I have a problem that uses multithreading and I’ve set it up to copy some integrators defining a differential equation problem. I’m running into issues with KLUFactorization(), namely my method fails when KLUFactorization(; reuse_symbolic = true) and is correct when reuse_symbolic = true. Below is an example that shows the issue (I don’t use atomics in the actual code, it’s just to make the point):

using OrdinaryDiffEq, LinearSolve, SparseArrays, Base.Threads
function ode_fnc(du::AbstractVector{T}, u, p, t) where {T}
    J, λ = p 
    fill!(du, zero(T))
    for i in eachindex(du)
        for j in J[i, :].nzind 
            du[i] += u[j] 
        end
        du[i] *=  λ * u[i]
    end
    return nothing 
end
function sum_integrator_serial(integrator, λ_rng)
    reinit!(integrator)
    s = 0.0
    for λ in λ_rng 
        integrator.p[2] = λ
        solve!(integrator)
        for j in eachindex(integrator.sol.u)
            s += sum(integrator.sol.u[j])
        end
        reinit!(integrator)
    end
    return s 
end
function sum_integrator_parallel(integrator, λ_rng)
    reinit!(integrator)
    nt = Threads.nthreads()
    integrators = [deepcopy(integrator) for _ in 1:nt]
    s = Threads.Atomic{Float64}(0.0)
    Threads.@threads for λ in λ_rng 
        id = Threads.threadid()
        integrators[id].p[2] = λ 
        solve!(integrators[id])
        for j in eachindex(integrators[id].sol.u)
            Threads.atomic_add!(s, sum(integrators[id].sol.u[j]))
        end
        reinit!(integrators[id])
    end
    return s[]
end
n, p = 200, 0.02 
J = sprandn(n,n,p)
u0, tspan = ones(n), (0.0, 0.1)
jac_prototype = deepcopy(J)
ode_fnc_2 = ODEFunction(ode_fnc; jac_prototype)
params = [J, 0.1]
prob = ODEProblem(ode_fnc_2, u0, tspan, params)
alg = TRBDF2(linsolve = KLUFactorization(; reuse_symbolic=false))
integrator = init(prob, alg; saveat = 0.01)
λ_rng = 0.1:0.001:0.5

s_serial = sum_integrator(integrator, λ_rng)
s_parallel = sum_integrator_parallel(integrator, λ_rng)

This gives s_serial and s_parallel to be essentially the same value. But instead if I have reuse_symbolic = true:

ERROR: TaskFailedException
Stacktrace:
 [1] [cut some stuff here for the character limit]

    nested task error: LinearAlgebra.SingularException(0)
    Stacktrace:
      [1] kluerror
        @ C:\Users\licer\.julia\packages\KLU\uZoWv\src\KLU.jl:105 [inlined]
      [2] kluerror(common::KLU.LibKLU.klu_l_common_struct)
        @ KLU C:\Users\licer\.julia\packages\KLU\uZoWv\src\KLU.jl:117
      [3] klu!(K::KLU.KLUFactorization{Float64, Int64}, nzval::Vector{Float64})
        @ KLU C:\Users\licer\.julia\packages\KLU\uZoWv\src\KLU.jl:616
      [4] klu!(K::KLU.KLUFactorization{Float64, Int64}, S::SparseMatrixCSC{Float64, Int64})
        @ KLU C:\Users\licer\.julia\packages\KLU\uZoWv\src\KLU.jl:640
      [5] solve(cache::LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}, alg::KLUFactorization; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:reltol,), Tuple{Float64}}})
        @ LinearSolve C:\Users\licer\.julia\packages\LinearSolve\dY1ir\src\factorization.jl:334
      [6] #solve#5
        @ C:\Users\licer\.julia\packages\LinearSolve\dY1ir\src\common.jl:155 [inlined]
      [7] #dolinsolve#3
        @ C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\misc_utils.jl:111 [inlined]
      [8] compute_step!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, integrator::OrdinaryDiffEq.ODEIntegrator{TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
        @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\nlsolve\newton.jl:302
      [9] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, 
Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, integrator::OrdinaryDiffEq.ODEIntegrator{TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, 
Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, 
Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, 
Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, repeat_step::Bool)
        @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\nlsolve\nlsolve.jl:48
     [10] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, 
Float64}, true, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, repeat_step::Bool)
        @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\perform_step\sdirk_perform_step.jl:481
     [11] perform_step!
        @ C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\perform_step\sdirk_perform_step.jl:458 [inlined]
     [12] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Vector{Float64}, Nothing, Float64, Vector{Any}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Any}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, TRBDF2{12, true, KLUFactorization, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Any}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
        @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\SRiFK\src\solve.jl:514
     [13] macro expansion
        @ c:\Users\licer\Documents\intg_testing.jl:34 [inlined]
     [14] [cut some stuff here]
        @ Main .\threadingconstructs.jl:84
     [15] #76#threadsfor_fun
        @ .\threadingconstructs.jl:51 [inlined]
     [16] [cut some stuff here for character limit]

I get the error above (I cut some of the errors at the start and at [16] and [14] due to character limit). I thought deepcopy should duplicate any of the caches used so that there is no communication. How can I fix it so that I can use reuse_symbolic = true?

2

There might be an issue with a ccall there. Open an issue. @rayegun make note of this.

1 Like
3

Will do. Would this be a LinearSolve.jl issue, then?

4

Yes, or maybe KLU.jl. It boils down to the fact that the symbolic factorization is memory that is part of the SuiteSparse factorization and that memory is I think created in C. So deepcopying must be just copying the pointer? We’d need to dig in but it’s some wrapper issue that hopefully we can isolate, but it at least isn’t a differential equation thing as much as it’s something in the KLU wrappers.

1 Like
5

I see. Thanks, I put it here Deepcopying KLUFactorization is different for reuse_symbolic=true vs reuse_symbolic=false · Issue #12 · JuliaSparse/KLU.jl · GitHub.