Problem
Solving against an existing sparse Cholesky factorization in a loop allocates on every solve, even with a preallocated solution and RHS and a reused factorization. For workloads with many solves against a fixed factorization, e.g. time-stepping, PDE-constrained optimization, eigensolvers, these per-solve allocations add up.
Example:
using SparseArrays, LinearAlgebra
N = 5000
A = sprand(N, N, 1e-3)
A = A*A' + N*I # sparse SPD
F = cholesky(A)
b = rand(N)
x = similar(b)
ldiv!(x, F, b) # warm-up
@allocated ldiv!(x, F, b) # Allocates
The allocations come from CHOLMOD allocating and freeing its internal Y/E workspace on each solve. cholmod_solve2 lets the caller pass persistent buffers instead.
This can be solved by calling cholmod_solve2 or cholmod_l_solve2 directly, see https://github.com/lorenzogambichler/CHOLMODSolve.jl as a quick proof of concept. But since the SuiteSparse version is tied to the Julia version, this solution is a bit fragile. This can be seen from CholmodSolve2.jl, a previous solution attempt, which has not been updated and is therefore unusable.
Proposed direction
Possibly a reusable workspace object that ldiv! or solve! can optionally accept. Thereby the internal Y and E buffers can be pre-allocated once before solving, similar to the repo above.
This issue has already come up a few times in the past, e.g. in the julialang discourse.
Problem
Solving against an existing sparse Cholesky factorization in a loop allocates on every solve, even with a preallocated solution and RHS and a reused factorization. For workloads with many solves against a fixed factorization, e.g. time-stepping, PDE-constrained optimization, eigensolvers, these per-solve allocations add up.
Example:
The allocations come from CHOLMOD allocating and freeing its internal
Y/Eworkspace on each solve.cholmod_solve2lets the caller pass persistent buffers instead.This can be solved by calling
cholmod_solve2orcholmod_l_solve2directly, see https://github.com/lorenzogambichler/CHOLMODSolve.jl as a quick proof of concept. But since the SuiteSparse version is tied to the Julia version, this solution is a bit fragile. This can be seen fromCholmodSolve2.jl, a previous solution attempt, which has not been updated and is therefore unusable.Proposed direction
Possibly a reusable workspace object that
ldiv!orsolve!can optionally accept. Thereby the internalYandEbuffers can be pre-allocated once before solving, similar to the repo above.This issue has already come up a few times in the past, e.g. in the julialang discourse.