From 004ed926c1fc692ee6d51f8694c4afa815473fb2 Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Sat, 20 Jun 2026 13:17:44 +0200 Subject: [PATCH 1/6] Add CholmodWorkspace for allocation-free CHOLMOD solves. Pre-allocates cholmod_dense_struct wrappers and Y/E buffers so cholmod_solve2 reuses them across calls instead of reallocating --- src/solvers/cholmod.jl | 114 ++++++++++++++++++++++++++++++++++++++++- test/cholmod.jl | 42 +++++++++++++++ 2 files changed, 155 insertions(+), 1 deletion(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 47082bcc..5d95d76a 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -25,6 +25,7 @@ import LinearAlgebra: (\), AdjointFactorization, using SparseArrays using SparseArrays: getcolptr, AbstractSparseVecOrMat export + CholmodWorkspace, Dense, Factor, Sparse @@ -1915,7 +1916,55 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO throw(ArgumentError("self-adjoint sparse system solve not implemented for sparse rhs B," * " consider to convert B to a dense array")) -# in-place ldiv! +""" + CholmodWorkspace() + +Reusable workspace for allocation-free calls to `ldiv!(x, L, b, ws)` with a +CHOLMOD [`Factor`](@ref). + +`cholmod_solve2` allocates two temporary dense matrices (Y and E) on every solve. +A `CholmodWorkspace` holds those buffers across calls so that CHOLMOD reuses them +instead of reallocating. It also pre-allocates the `cholmod_dense_struct` wrappers +for `x` and `b`, eliminating all Julia-side heap allocations per solve. + +The workspace is finalizer-protected: the CHOLMOD-managed Y and E buffers are freed +automatically when the workspace is garbage-collected. A single workspace can be +reused across calls with different arrays or a different number of RHS columns; +CHOLMOD will resize Y and E as needed. + +# Example +```julia +F = cholesky(A) +ws = CholmodWorkspace() +for b in rhs_list + ldiv!(x, F, b, ws) +end +``` +""" +mutable struct CholmodWorkspace + dense_x::cholmod_dense_struct + dense_b::cholmod_dense_struct + X::Ref{Ptr{cholmod_dense_struct}} + Y::Ref{Ptr{cholmod_dense_struct}} + E::Ref{Ptr{cholmod_dense_struct}} + + function CholmodWorkspace() + ws = new( + cholmod_dense_struct(), # all fields written in ldiv! before use + cholmod_dense_struct(), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)) + ) + finalizer(ws) do w + w.Y[] != C_NULL && free!(w.Y[]) + w.E[] != C_NULL && free!(w.E[]) + end + ws + end +end + +# in-place ldiv! and optional allocation-free ldiv! for TI in IndexTypes @eval function ldiv!(x::StridedVecOrMat{T}, L::Factor{T, $TI}, @@ -1973,6 +2022,69 @@ for TI in IndexTypes return x end + + @eval function ldiv!(x::StridedVecOrMat{T}, + L::Factor{T, $TI}, + b::StridedVecOrMat{T}, + ws::CholmodWorkspace) where {T<:VTypes} + if x === b + throw(ArgumentError("output array must not be aliased with input array")) + end + if size(L, 1) != size(b, 1) + throw(DimensionMismatch("Factorization and RHS should have the same number of rows. " * + "Factorization has $(size(L, 2)) rows, but RHS has $(size(b, 1)) rows.")) + end + if size(L, 2) != size(x, 1) + throw(DimensionMismatch("Factorization and solution should match sizes. " * + "Factorization has $(size(L, 1)) columns, but solution has $(size(x, 1)) rows.")) + end + if size(x, 2) != size(b, 2) + throw(DimensionMismatch("Solution and RHS should have the same number of columns. " * + "Solution has $(size(x, 2)) columns, but RHS has $(size(b, 2)) columns.")) + end + if !issuccess(L) + s = unsafe_load(pointer(L)) + if s.is_ll == 1 + throw(LinearAlgebra.PosDefException(s.minor)) + else + throw(LinearAlgebra.ZeroPivotException(s.minor)) + end + end + + # Update all fields of the reused cholmod_dense_structs on every call, + # since T, dimensions, and data pointers may all change between calls. + ws.dense_x.nrow = size(x, 1) + ws.dense_x.ncol = size(x, 2) + ws.dense_x.nzmax = length(x) + ws.dense_x.d = stride(x, 2) + ws.dense_x.x = pointer(x) + ws.dense_x.z = C_NULL + ws.dense_x.xtype = xtyp(T) + ws.dense_x.dtype = dtyp(T) + + ws.dense_b.nrow = size(b, 1) + ws.dense_b.ncol = size(b, 2) + ws.dense_b.nzmax = length(b) + ws.dense_b.d = stride(b, 2) + ws.dense_b.x = pointer(b) + ws.dense_b.z = C_NULL + ws.dense_b.xtype = xtyp(T) + ws.dense_b.dtype = dtyp(T) + + ws.X[] = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_x)) + Bptr = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_b)) + status = GC.@preserve x b ws begin + $(cholname(:solve2, TI))( + CHOLMOD_A, L, + Bptr, C_NULL, + ws.X, C_NULL, + ws.Y, ws.E, + getcommon($TI)) + end + @assert !iszero(status) + + return x + end end ## Other convenience methods diff --git a/test/cholmod.jl b/test/cholmod.jl index 95772144..98f4a9df 100644 --- a/test/cholmod.jl +++ b/test/cholmod.jl @@ -315,6 +315,48 @@ end @test_throws DimensionMismatch ldiv!(x2, factor, B) end +@testset "ldiv! with CholmodWorkspace $Tv $Ti" begin + local A, x, x2, b, X, X2, B + A = sprand(10, 10, 0.1) + A = I + A * A' + A = convert(SparseMatrixCSC{Tv,Ti}, A) + factor = cholesky(A) + ws = CHOLMOD.CholmodWorkspace() + + x = fill(Tv(1), 10) + b = A * x + x2 = zero(x) + ldiv!(x2, factor, b, ws) + @test x2 ≈ x + + # reuse workspace + X = fill(Tv(1), 10, 5) + B = A * X + X2 = zero(X) + ldiv!(X2, factor, B, ws) + @test X2 ≈ X + + # workspace is reusable across calls + fill!(x2, 0) + ldiv!(x2, factor, b, ws) + @test x2 ≈ x + + # allocation reduction + allocs = @allocated ldiv!(x2, factor, b, ws) + @test allocs <= 16 + + c = fill(Tv(1), size(x, 1) + 1) + C = fill(Tv(1), size(X, 1) + 1, size(X, 2)) + y = fill(Tv(1), size(x, 1) + 1) + Y = fill(Tv(1), size(X, 1) + 1, size(X, 2)) + @test_throws DimensionMismatch ldiv!(y, factor, b, ws) + @test_throws DimensionMismatch ldiv!(Y, factor, B, ws) + @test_throws DimensionMismatch ldiv!(x2, factor, c, ws) + @test_throws DimensionMismatch ldiv!(X2, factor, C, ws) + @test_throws DimensionMismatch ldiv!(X2, factor, b, ws) + @test_throws DimensionMismatch ldiv!(x2, factor, B, ws) +end + end #end for Ti ∈ itypes for Tv ∈ (Float32, Float64) From d20ab334cbdac0603e826f5c74025cd598242e55 Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Sat, 20 Jun 2026 13:28:25 +0200 Subject: [PATCH 2/6] minor spelling changes --- src/solvers/cholmod.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 5d95d76a..47758940 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -1919,28 +1919,29 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO """ CholmodWorkspace() -Reusable workspace for allocation-free calls to `ldiv!(x, L, b, ws)` with a -CHOLMOD [`Factor`](@ref). +Reusable workspace for allocation-free calls to `ldiv!(x, L, b, ws)`. `cholmod_solve2` allocates two temporary dense matrices (Y and E) on every solve. A `CholmodWorkspace` holds those buffers across calls so that CHOLMOD reuses them instead of reallocating. It also pre-allocates the `cholmod_dense_struct` wrappers -for `x` and `b`, eliminating all Julia-side heap allocations per solve. +for `x` and `b`. +This eliminates all Julia-side heap allocations per solve. -The workspace is finalizer-protected: the CHOLMOD-managed Y and E buffers are freed -automatically when the workspace is garbage-collected. A single workspace can be -reused across calls with different arrays or a different number of RHS columns; -CHOLMOD will resize Y and E as needed. +A finalizer frees Y adn E automatically when the worspace is garbage-collected. + +A single workspace can be reused across calls with different arrays or a different number of RHS columns since +CHOLMOD will resize Y and E. # Example ```julia -F = cholesky(A) +F = cholesky(A) ws = CholmodWorkspace() for b in rhs_list ldiv!(x, F, b, ws) end ``` """ + mutable struct CholmodWorkspace dense_x::cholmod_dense_struct dense_b::cholmod_dense_struct @@ -1950,7 +1951,7 @@ mutable struct CholmodWorkspace function CholmodWorkspace() ws = new( - cholmod_dense_struct(), # all fields written in ldiv! before use + cholmod_dense_struct(), # all fields will be written in ldiv! before use cholmod_dense_struct(), Ref(Ptr{cholmod_dense_struct}(C_NULL)), Ref(Ptr{cholmod_dense_struct}(C_NULL)), @@ -2051,8 +2052,7 @@ for TI in IndexTypes end end - # Update all fields of the reused cholmod_dense_structs on every call, - # since T, dimensions, and data pointers may all change between calls. + # Update all cholmod_dense_structs fields on every call in case dimensions or T change ws.dense_x.nrow = size(x, 1) ws.dense_x.ncol = size(x, 2) ws.dense_x.nzmax = length(x) From d1b785c886f2be75ef496f70f9de398dcfb75afd Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Tue, 7 Jul 2026 20:35:16 +0200 Subject: [PATCH 3/6] Add workspace to Factor to keep 3-arg API. Redesign CholmodWorkspace similar to UmfpackWS inside UmfpackLU --- src/solvers/cholmod.jl | 271 +++++++++++++++-------------------------- test/cholmod.jl | 47 +++---- 2 files changed, 113 insertions(+), 205 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 47758940..83d2689a 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -25,7 +25,6 @@ import LinearAlgebra: (\), AdjointFactorization, using SparseArrays using SparseArrays: getcolptr, AbstractSparseVecOrMat export - CholmodWorkspace, Dense, Factor, Sparse @@ -259,8 +258,50 @@ function Sparse(p::Ptr{cholmod_sparse}) Sparse{jlxtype(s.xtype, s.dtype)}(p) end +""" + CholmodWorkspace() + +Reusable workspace for allocation-free CHOLMOD solves. + +`cholmod_solve2` allocates two temporary dense matrices (Y and E) on every solve. +A `CholmodWorkspace` holds those buffers across calls so that CHOLMOD reuses them +instead of reallocating. It also pre-allocates the `cholmod_dense_struct` wrappers +for `x` and `b`, eliminating all Julia-side heap allocations per solve. + +Each [`Factor`](@ref) embeds its own workspace and uses it automatically in +`ldiv!(x, F, b)`. For multi-thread solves with a single factor, use `copy(F)` +for every thread so each copy gets its own independent workspace and lock. + +The workspace is finalizer-protected. The Y and E buffers are freed automatically +when the workspace is garbage-collected. +""" +mutable struct CholmodWorkspace + dense_x::cholmod_dense_struct + dense_b::cholmod_dense_struct + X::Ref{Ptr{cholmod_dense_struct}} + Y::Ref{Ptr{cholmod_dense_struct}} + E::Ref{Ptr{cholmod_dense_struct}} + + function CholmodWorkspace() + ws = new( + cholmod_dense_struct(), # all fields will be written in ldiv! before use + cholmod_dense_struct(), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)) + ) + finalizer(ws) do w + w.Y[] != C_NULL && free!(w.Y[]) + w.E[] != C_NULL && free!(w.E[]) + end + ws + end +end + mutable struct Factor{Tv<:VTypes, Ti<:ITypes} <: Factorization{Tv} ptr::Ptr{cholmod_factor} + workspace::CholmodWorkspace + lock::ReentrantLock function Factor{Tv, Ti}(ptr::Ptr{cholmod_factor}, register_finalizer = true) where {Tv, Ti} if ptr == C_NULL throw(ArgumentError("factorization construction failed for " * @@ -277,7 +318,7 @@ mutable struct Factor{Tv<:VTypes, Ti<:ITypes} <: Factorization{Tv} free!(ptr, Ti) throw(CHOLMODException("dtype=$(dtyp(Tv)) not supported")) end - F = new(ptr) + F = new(ptr, CholmodWorkspace(), ReentrantLock()) if register_finalizer finalizer(free!, F) end @@ -312,6 +353,8 @@ Base.pointer(x::Dense{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_dense}, Base.pointer(x::Sparse{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_sparse}, x) Base.pointer(x::Factor{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_factor}, x) +getworkspace(F::Factor) = F.workspace + function _factor_components(is_ll) is_ll ? (:L, :U, :PtL, :UP) : (:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP) end @@ -858,32 +901,6 @@ function Base.convert(::Type{Dense{Tnew}}, A::Dense{T}) where {Tnew, T} end Base.convert(::Type{Dense{T}}, A::Dense{T}) where T = A -# Just calling Dense(x) or Dense(b) will allocate new -# `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse -# the existing memory. We can do this by creating new -# `cholmod_dense_struct`s and filling them manually. -function wrap_dense_and_ptr(x::StridedVecOrMat{T}) where {T <: VTypes} - dense_x = cholmod_dense_struct() - dense_x.nrow = size(x, 1) - dense_x.ncol = size(x, 2) - dense_x.nzmax = length(x) - dense_x.d = stride(x, 2) - dense_x.x = pointer(x) - dense_x.z = C_NULL - dense_x.xtype = xtyp(eltype(x)) - dense_x.dtype = dtyp(eltype(x)) - return dense_x, pointer_from_objref(dense_x) -end -# We need to use a special handling for the case of `Dense` -# input arrays since the `pointer` refers to the pointer to the -# `cholmod_dense`, not to the array values themselves as for -# standard arrays. -function wrap_dense_and_ptr(x::Dense{T}) where {T <: VTypes} - dense_x_ptr = x.ptr - dense_x = unsafe_load(dense_x_ptr) - return dense_x, pointer_from_objref(dense_x) -end - # This constructor assumes zero based colptr and rowval function Sparse(m::Integer, n::Integer, colptr0::Vector{Ti}, rowval0::Vector{Ti}, @@ -1326,8 +1343,8 @@ end @inline function getproperty(F::Factor, sym::Symbol) if sym === :p return get_perm(F) - elseif sym === :ptr - return getfield(F, :ptr) + elseif sym === :ptr || sym === :workspace || sym === :lock # add workspace and lock + return getfield(F, sym) else return FactorComponent(F, sym) end @@ -1434,11 +1451,11 @@ end function cholesky!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0, check::Bool = true) where Tv - # Compute the numerical factorization - @cholmod_param final_ll = true begin - factorize_p!(A, shift, F) + @lock F.lock begin + @cholmod_param final_ll = true begin + factorize_p!(A, shift, F) + end end - check && (issuccess(F) || throw(LinearAlgebra.PosDefException(1))) return F end @@ -1609,12 +1626,10 @@ LinearAlgebra._cholesky(A::Union{SparseMatrixCSC{T}, SparseMatrixCSC{Complex{T}} function ldlt!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0, check::Bool = true) where Tv - # Makes it an LDLt - change_factor!(F, false, false, true, false) - - # Compute the numerical factorization - factorize_p!(A, shift, F) - + @lock F.lock begin + change_factor!(F, false, false, true, false) + factorize_p!(A, shift, F) + end check && (issuccess(F) || throw(LinearAlgebra.ZeroPivotException(1))) return F end @@ -1916,118 +1931,53 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO throw(ArgumentError("self-adjoint sparse system solve not implemented for sparse rhs B," * " consider to convert B to a dense array")) -""" - CholmodWorkspace() - -Reusable workspace for allocation-free calls to `ldiv!(x, L, b, ws)`. - -`cholmod_solve2` allocates two temporary dense matrices (Y and E) on every solve. -A `CholmodWorkspace` holds those buffers across calls so that CHOLMOD reuses them -instead of reallocating. It also pre-allocates the `cholmod_dense_struct` wrappers -for `x` and `b`. -This eliminates all Julia-side heap allocations per solve. - -A finalizer frees Y adn E automatically when the worspace is garbage-collected. - -A single workspace can be reused across calls with different arrays or a different number of RHS columns since -CHOLMOD will resize Y and E. - -# Example -```julia -F = cholesky(A) -ws = CholmodWorkspace() -for b in rhs_list - ldiv!(x, F, b, ws) -end -``` -""" - -mutable struct CholmodWorkspace - dense_x::cholmod_dense_struct - dense_b::cholmod_dense_struct - X::Ref{Ptr{cholmod_dense_struct}} - Y::Ref{Ptr{cholmod_dense_struct}} - E::Ref{Ptr{cholmod_dense_struct}} +# Julia array -> fill ws.dense_b fields and return pointer to it +# CHOLMOD Dense object -> struct is already set up by CHOLMOD, return b.ptr directly +@inline function _setup_bptr(b::StridedVecOrMat{T}, ws::CholmodWorkspace) where T<:VTypes + ws.dense_b.nrow = size(b, 1) + ws.dense_b.ncol = size(b, 2) + ws.dense_b.nzmax = length(b) + ws.dense_b.d = stride(b, 2) + ws.dense_b.x = pointer(b) + ws.dense_b.z = C_NULL + ws.dense_b.xtype = xtyp(T) + ws.dense_b.dtype = dtyp(T) + return Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_b)) +end +@inline _setup_bptr(b::Dense{<:VTypes}, ::CholmodWorkspace) = b.ptr - function CholmodWorkspace() - ws = new( - cholmod_dense_struct(), # all fields will be written in ldiv! before use - cholmod_dense_struct(), - Ref(Ptr{cholmod_dense_struct}(C_NULL)), - Ref(Ptr{cholmod_dense_struct}(C_NULL)), - Ref(Ptr{cholmod_dense_struct}(C_NULL)) - ) - finalizer(ws) do w - w.Y[] != C_NULL && free!(w.Y[]) - w.E[] != C_NULL && free!(w.E[]) - end - ws - end -end - -# in-place ldiv! and optional allocation-free ldiv! for TI in IndexTypes - @eval function ldiv!(x::StridedVecOrMat{T}, - L::Factor{T, $TI}, - b::StridedVecOrMat{T}) where {T<:VTypes} - if x === b - throw(ArgumentError("output array must not be aliased with input array")) - end - if size(L, 1) != size(b, 1) - throw(DimensionMismatch("Factorization and RHS should have the same number of rows. " * - "Factorization has $(size(L, 2)) rows, but RHS has $(size(b, 1)) rows.")) - end - if size(L, 2) != size(x, 1) - throw(DimensionMismatch("Factorization and solution should match sizes. " * - "Factorization has $(size(L, 1)) columns, but solution has $(size(x, 1)) rows.")) - end - if size(x, 2) != size(b, 2) - throw(DimensionMismatch("Solution and RHS should have the same number of columns. " * - "Solution has $(size(x, 2)) columns, but RHS has $(size(b, 2)) columns.")) - end - if !issuccess(L) - s = unsafe_load(pointer(L)) - if s.is_ll == 1 - throw(LinearAlgebra.PosDefException(s.minor)) - else - throw(LinearAlgebra.ZeroPivotException(s.minor)) + @eval function solve!(x::StridedVecOrMat{T}, + L::Factor{T, $TI}, + b::StridedVecOrMat{T}; + ws = getworkspace(L)) where {T<:VTypes} + @lock L.lock begin + ws.dense_x.nrow = size(x, 1) + ws.dense_x.ncol = size(x, 2) + ws.dense_x.nzmax = length(x) + ws.dense_x.d = stride(x, 2) + ws.dense_x.x = pointer(x) + ws.dense_x.z = C_NULL + ws.dense_x.xtype = xtyp(T) + ws.dense_x.dtype = dtyp(T) + + ws.X[] = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_x)) + Bptr = _setup_bptr(b, ws) + status = GC.@preserve x b ws begin + $(cholname(:solve2, TI))( + CHOLMOD_A, L, + Bptr, C_NULL, + ws.X, C_NULL, + ws.Y, ws.E, + getcommon($TI)) end + @assert !iszero(status) end - - # Just calling Dense(x) or Dense(b) will allocate new - # `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse - # the existing memory. We can do this by creating new - # `cholmod_dense_struct`s and filling them manually. - dense_x, dense_x_ptr = wrap_dense_and_ptr(x) - dense_b, dense_b_ptr = wrap_dense_and_ptr(b) - - X_Handle = Ptr{cholmod_dense_struct}(dense_x_ptr) - Y_Handle = Ptr{cholmod_dense_struct}(C_NULL) - E_Handle = Ptr{cholmod_dense_struct}(C_NULL) - status = GC.@preserve x dense_x b dense_b begin - $(cholname(:solve2, TI))( - CHOLMOD_A, L, - Ref(dense_b), C_NULL, - Ref(X_Handle), C_NULL, - Ref(Y_Handle), - Ref(E_Handle), - getcommon($TI)) - end - if Y_Handle != C_NULL - free!(Y_Handle) - end - if E_Handle != C_NULL - free!(E_Handle) - end - @assert !iszero(status) - - return x end @eval function ldiv!(x::StridedVecOrMat{T}, L::Factor{T, $TI}, - b::StridedVecOrMat{T}, - ws::CholmodWorkspace) where {T<:VTypes} + b::StridedVecOrMat{T}) where {T<:VTypes} if x === b throw(ArgumentError("output array must not be aliased with input array")) end @@ -2051,38 +2001,7 @@ for TI in IndexTypes throw(LinearAlgebra.ZeroPivotException(s.minor)) end end - - # Update all cholmod_dense_structs fields on every call in case dimensions or T change - ws.dense_x.nrow = size(x, 1) - ws.dense_x.ncol = size(x, 2) - ws.dense_x.nzmax = length(x) - ws.dense_x.d = stride(x, 2) - ws.dense_x.x = pointer(x) - ws.dense_x.z = C_NULL - ws.dense_x.xtype = xtyp(T) - ws.dense_x.dtype = dtyp(T) - - ws.dense_b.nrow = size(b, 1) - ws.dense_b.ncol = size(b, 2) - ws.dense_b.nzmax = length(b) - ws.dense_b.d = stride(b, 2) - ws.dense_b.x = pointer(b) - ws.dense_b.z = C_NULL - ws.dense_b.xtype = xtyp(T) - ws.dense_b.dtype = dtyp(T) - - ws.X[] = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_x)) - Bptr = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_b)) - status = GC.@preserve x b ws begin - $(cholname(:solve2, TI))( - CHOLMOD_A, L, - Bptr, C_NULL, - ws.X, C_NULL, - ws.Y, ws.E, - getcommon($TI)) - end - @assert !iszero(status) - + solve!(x, L, b) return x end end diff --git a/test/cholmod.jl b/test/cholmod.jl index 98f4a9df..7fdbeeb0 100644 --- a/test/cholmod.jl +++ b/test/cholmod.jl @@ -303,6 +303,15 @@ end @inferred ldiv!(X2, factor, B) @test X2 ≈ X + # reuse across multiple calls (Y/E buffers kept in workspace) + fill!(x2, 0) + ldiv!(x2, factor, b) + @test x2 ≈ x + + # Y/E buffers are reused across calls, remaining 16 bytes come from CHOLMOD internals + allocs = @allocated ldiv!(x2, factor, b) + @test allocs <= 16 + c = fill(Tv(1), size(x, 1) + 1) C = fill(Tv(1), size(X, 1) + 1, size(X, 2)) y = fill(Tv(1), size(x, 1) + 1) @@ -315,46 +324,26 @@ end @test_throws DimensionMismatch ldiv!(x2, factor, B) end -@testset "ldiv! with CholmodWorkspace $Tv $Ti" begin - local A, x, x2, b, X, X2, B +@testset "copy(Factor) workspace isolation $Tv $Ti" begin + local A, x, b, x2, x3 A = sprand(10, 10, 0.1) A = I + A * A' A = convert(SparseMatrixCSC{Tv,Ti}, A) factor = cholesky(A) - ws = CHOLMOD.CholmodWorkspace() + factor2 = copy(factor) x = fill(Tv(1), 10) b = A * x x2 = zero(x) - ldiv!(x2, factor, b, ws) - @test x2 ≈ x + x3 = zero(x) - # reuse workspace - X = fill(Tv(1), 10, 5) - B = A * X - X2 = zero(X) - ldiv!(X2, factor, B, ws) - @test X2 ≈ X - - # workspace is reusable across calls - fill!(x2, 0) - ldiv!(x2, factor, b, ws) + ldiv!(x2, factor, b) + ldiv!(x3, factor2, b) @test x2 ≈ x + @test x3 ≈ x - # allocation reduction - allocs = @allocated ldiv!(x2, factor, b, ws) - @test allocs <= 16 - - c = fill(Tv(1), size(x, 1) + 1) - C = fill(Tv(1), size(X, 1) + 1, size(X, 2)) - y = fill(Tv(1), size(x, 1) + 1) - Y = fill(Tv(1), size(X, 1) + 1, size(X, 2)) - @test_throws DimensionMismatch ldiv!(y, factor, b, ws) - @test_throws DimensionMismatch ldiv!(Y, factor, B, ws) - @test_throws DimensionMismatch ldiv!(x2, factor, c, ws) - @test_throws DimensionMismatch ldiv!(X2, factor, C, ws) - @test_throws DimensionMismatch ldiv!(X2, factor, b, ws) - @test_throws DimensionMismatch ldiv!(x2, factor, B, ws) + # Verify each copy has its own independent workspace + @test factor.workspace !== factor2.workspace end end #end for Ti ∈ itypes From 6b028732a63e722a1c4d3c5e53b58f623967225d Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Wed, 8 Jul 2026 13:13:51 +0200 Subject: [PATCH 4/6] Remove NULL check before free! --- src/solvers/cholmod.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 83d2689a..1e609c8a 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -291,8 +291,8 @@ mutable struct CholmodWorkspace Ref(Ptr{cholmod_dense_struct}(C_NULL)) ) finalizer(ws) do w - w.Y[] != C_NULL && free!(w.Y[]) - w.E[] != C_NULL && free!(w.E[]) + free!(w.Y[]) + free!(w.E[]) end ws end From 143e0698d1d2c4d7113cf00f6c468a862c4aef92 Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Wed, 8 Jul 2026 19:04:23 +0200 Subject: [PATCH 5/6] Add CholmodWorkspace fields directly to Factor --- src/solvers/cholmod.jl | 126 +++++++++++++++++------------------------ test/cholmod.jl | 6 +- 2 files changed, 54 insertions(+), 78 deletions(-) diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 1e609c8a..9d5ce1dc 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -258,49 +258,15 @@ function Sparse(p::Ptr{cholmod_sparse}) Sparse{jlxtype(s.xtype, s.dtype)}(p) end -""" - CholmodWorkspace() - -Reusable workspace for allocation-free CHOLMOD solves. - -`cholmod_solve2` allocates two temporary dense matrices (Y and E) on every solve. -A `CholmodWorkspace` holds those buffers across calls so that CHOLMOD reuses them -instead of reallocating. It also pre-allocates the `cholmod_dense_struct` wrappers -for `x` and `b`, eliminating all Julia-side heap allocations per solve. - -Each [`Factor`](@ref) embeds its own workspace and uses it automatically in -`ldiv!(x, F, b)`. For multi-thread solves with a single factor, use `copy(F)` -for every thread so each copy gets its own independent workspace and lock. - -The workspace is finalizer-protected. The Y and E buffers are freed automatically -when the workspace is garbage-collected. -""" -mutable struct CholmodWorkspace +# Factor stores its own temporary CHOLMOD Y/E buffers for use in ldiv! +# and pre-allocates cholmod_dense_struct wrappers +mutable struct Factor{Tv<:VTypes, Ti<:ITypes} <: Factorization{Tv} + ptr::Ptr{cholmod_factor} dense_x::cholmod_dense_struct dense_b::cholmod_dense_struct X::Ref{Ptr{cholmod_dense_struct}} Y::Ref{Ptr{cholmod_dense_struct}} E::Ref{Ptr{cholmod_dense_struct}} - - function CholmodWorkspace() - ws = new( - cholmod_dense_struct(), # all fields will be written in ldiv! before use - cholmod_dense_struct(), - Ref(Ptr{cholmod_dense_struct}(C_NULL)), - Ref(Ptr{cholmod_dense_struct}(C_NULL)), - Ref(Ptr{cholmod_dense_struct}(C_NULL)) - ) - finalizer(ws) do w - free!(w.Y[]) - free!(w.E[]) - end - ws - end -end - -mutable struct Factor{Tv<:VTypes, Ti<:ITypes} <: Factorization{Tv} - ptr::Ptr{cholmod_factor} - workspace::CholmodWorkspace lock::ReentrantLock function Factor{Tv, Ti}(ptr::Ptr{cholmod_factor}, register_finalizer = true) where {Tv, Ti} if ptr == C_NULL @@ -318,9 +284,15 @@ mutable struct Factor{Tv<:VTypes, Ti<:ITypes} <: Factorization{Tv} free!(ptr, Ti) throw(CHOLMODException("dtype=$(dtyp(Tv)) not supported")) end - F = new(ptr, CholmodWorkspace(), ReentrantLock()) + F = new(ptr, + cholmod_dense_struct(), + cholmod_dense_struct(), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + Ref(Ptr{cholmod_dense_struct}(C_NULL)), + ReentrantLock()) if register_finalizer - finalizer(free!, F) + finalizer(free!, F) # includes Y/E buffers end return F end @@ -353,8 +325,6 @@ Base.pointer(x::Dense{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_dense}, Base.pointer(x::Sparse{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_sparse}, x) Base.pointer(x::Factor{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_factor}, x) -getworkspace(F::Factor) = F.workspace - function _factor_components(is_ll) is_ll ? (:L, :U, :PtL, :UP) : (:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP) end @@ -1258,7 +1228,11 @@ end free!(A::Dense) = free!(pointer(A)) free!(A::Sparse{<:Any, Ti}) where Ti = free!(pointer(A), Ti) -free!(F::Factor{<:Any, Ti}) where Ti = free!(pointer(F), Ti) +function free!(F::Factor{<:Any, Ti}) where Ti + free!(getfield(F, :Y)[]) + free!(getfield(F, :E)[]) + free!(pointer(F), Ti) +end nnz(F::Factor) = nnz(Sparse(F)) @@ -1343,7 +1317,7 @@ end @inline function getproperty(F::Factor, sym::Symbol) if sym === :p return get_perm(F) - elseif sym === :ptr || sym === :workspace || sym === :lock # add workspace and lock + elseif sym === :ptr || sym === :lock return getfield(F, sym) else return FactorComponent(F, sym) @@ -1931,44 +1905,46 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO throw(ArgumentError("self-adjoint sparse system solve not implemented for sparse rhs B," * " consider to convert B to a dense array")) -# Julia array -> fill ws.dense_b fields and return pointer to it +# Julia array -> fill dense_b fields and return pointer to it # CHOLMOD Dense object -> struct is already set up by CHOLMOD, return b.ptr directly -@inline function _setup_bptr(b::StridedVecOrMat{T}, ws::CholmodWorkspace) where T<:VTypes - ws.dense_b.nrow = size(b, 1) - ws.dense_b.ncol = size(b, 2) - ws.dense_b.nzmax = length(b) - ws.dense_b.d = stride(b, 2) - ws.dense_b.x = pointer(b) - ws.dense_b.z = C_NULL - ws.dense_b.xtype = xtyp(T) - ws.dense_b.dtype = dtyp(T) - return Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_b)) -end -@inline _setup_bptr(b::Dense{<:VTypes}, ::CholmodWorkspace) = b.ptr +@inline function _setup_bptr(b::StridedVecOrMat{T}, dense_b::cholmod_dense_struct) where T<:VTypes + dense_b.nrow = size(b, 1) + dense_b.ncol = size(b, 2) + dense_b.nzmax = length(b) + dense_b.d = stride(b, 2) + dense_b.x = pointer(b) + dense_b.z = C_NULL + dense_b.xtype = xtyp(T) + dense_b.dtype = dtyp(T) + return Ptr{cholmod_dense_struct}(pointer_from_objref(dense_b)) +end +@inline _setup_bptr(b::Dense{<:VTypes}, ::cholmod_dense_struct) = b.ptr for TI in IndexTypes - @eval function solve!(x::StridedVecOrMat{T}, - L::Factor{T, $TI}, - b::StridedVecOrMat{T}; - ws = getworkspace(L)) where {T<:VTypes} + @eval function solve!(x::StridedVecOrMat{T}, L::Factor{T, $TI}, b::StridedVecOrMat{T}) where {T<:VTypes} @lock L.lock begin - ws.dense_x.nrow = size(x, 1) - ws.dense_x.ncol = size(x, 2) - ws.dense_x.nzmax = length(x) - ws.dense_x.d = stride(x, 2) - ws.dense_x.x = pointer(x) - ws.dense_x.z = C_NULL - ws.dense_x.xtype = xtyp(T) - ws.dense_x.dtype = dtyp(T) - - ws.X[] = Ptr{cholmod_dense_struct}(pointer_from_objref(ws.dense_x)) - Bptr = _setup_bptr(b, ws) - status = GC.@preserve x b ws begin + dense_x = getfield(L, :dense_x) + X = getfield(L, :X) + Y = getfield(L, :Y) + E = getfield(L, :E) + + dense_x.nrow = size(x, 1) + dense_x.ncol = size(x, 2) + dense_x.nzmax = length(x) + dense_x.d = stride(x, 2) + dense_x.x = pointer(x) + dense_x.z = C_NULL + dense_x.xtype = xtyp(T) + dense_x.dtype = dtyp(T) + + X[] = Ptr{cholmod_dense_struct}(pointer_from_objref(dense_x)) + Bptr = _setup_bptr(b, getfield(L, :dense_b)) + status = GC.@preserve x b L begin $(cholname(:solve2, TI))( CHOLMOD_A, L, Bptr, C_NULL, - ws.X, C_NULL, - ws.Y, ws.E, + X, C_NULL, + Y, E, getcommon($TI)) end @assert !iszero(status) diff --git a/test/cholmod.jl b/test/cholmod.jl index 7fdbeeb0..8f7969ec 100644 --- a/test/cholmod.jl +++ b/test/cholmod.jl @@ -324,7 +324,7 @@ end @test_throws DimensionMismatch ldiv!(x2, factor, B) end -@testset "copy(Factor) workspace isolation $Tv $Ti" begin +@testset "copy(Factor) buffer isolation $Tv $Ti" begin local A, x, b, x2, x3 A = sprand(10, 10, 0.1) A = I + A * A' @@ -342,8 +342,8 @@ end @test x2 ≈ x @test x3 ≈ x - # Verify each copy has its own independent workspace - @test factor.workspace !== factor2.workspace + # Verify each copy has its own independent buffers + @test getfield(factor, :Y) !== getfield(factor2, :Y) end end #end for Ti ∈ itypes From 974a065c42775bfddd367e2c6a4d00e3d21946f5 Mon Sep 17 00:00:00 2001 From: Lorenzo Gambichler Date: Fri, 24 Jul 2026 16:56:55 +0200 Subject: [PATCH 6/6] add memory-leak test --- test/cholmod.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/test/cholmod.jl b/test/cholmod.jl index 8f7969ec..d6b91b1e 100644 --- a/test/cholmod.jl +++ b/test/cholmod.jl @@ -324,6 +324,25 @@ end @test_throws DimensionMismatch ldiv!(x2, factor, B) end +@testset "ldiv! no memory leak $Tv $Ti" begin + local A, b, x, F + A = sprand(10, 10, 0.1) + A = I + A * A' + A = convert(SparseMatrixCSC{Tv,Ti}, A) + F = cholesky(A) + b = A * fill(Tv(1), 10) + x = zero(b) + + ldiv!(x, F, b) # allocate buffers + GC.gc() + before = getcommon(Ti)[].memory_inuse + for _ in 1:1000 + ldiv!(x, F, b) + end + after = getcommon(Ti)[].memory_inuse + @test before == after +end + @testset "copy(Factor) buffer isolation $Tv $Ti" begin local A, x, b, x2, x3 A = sprand(10, 10, 0.1)