# Fixed-structure CHOLMOD with in-place \`set\_nonzeros!\`
**URL:** https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922
**Category:** Performance
**Tags:** memory-allocation, cholmod, sparsearrays
**Created:** [July 21, 2025, 8:56pm UTC](https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922 "2025-07-21T20:56:56Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@luke-kiernan](https://discourse.julialang.org/u/luke-kiernan)
#### Post date: [July 21, 2025, 8:56pm UTC](https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922/1 "2025-07-21T20:56:56Z")
I have a large sparse matrix. Its nonzeros change over the course of the problem, but its structure doesn’t. I want to do repeated in-place Cholesky factorizations, re-using the symbolic factorization and modifying the numeric factorization in-place. I can do the first factorization just fine:
```julia
using SparseArrays
using LinearAlgebra
A = sprand(10, 10, 0.3)
B = A + A' + 10 * LinearAlgebra.I(10)
B_cholmod = SparseArrays.CHOLMOD.Sparse(B)
sym = SparseArrays.CHOLMOD.symbolic(B_cholmod)
SparseArrays.CHOLMOD.cholesky!(sym, B_cholmod)
```
However, Julia’s interface does not allow me to write to the array of nonzero values of `B_cholmod`. So I’m implementing my own wrapper. Here’s what I have:
```julia
using SparseArrays
using LinearAlgebra
const VTypes = SparseArrays.CHOLMOD.VRealTypes
const ITypes = SparseArrays.CHOLMOD.ITypes
struct FixedStructureCHOLMOD{T, I} <: AbstractSparseMatrix{T, I}
_mat::SparseArrays.CHOLMOD.Sparse{T, I}
_values::Vector{T}
end
function FixedStructureCHOLMOD(A::SparseMatrixCSC{T, I}) where {T<:VTypes, I<:ITypes}
cholmod_mat = SparseArrays.CHOLMOD.Sparse(A, -1) # -1 for hermitian
chol = unsafe_load(cholmod_mat.ptr)
x_ptr = chol.x
# max capacity may not be the same as the number of nonzeros:
# instead of chol.nzmax, we want chol.p[chol.ncol], in C indexing.
# (what about packed matrices, though?)
nonzero_count = unsafe_load(Ptr{I}(chol.p), chol.ncol + 1)
values = unsafe_wrap(Vector{T}, Ptr{T}(x_ptr), nonzero_count)
return FixedStructureCHOLMOD(cholmod_mat, values)
end
"""
set_values!(mat::FixedStructureCHOLMOD, new_vals::AbstractVector{Float64})
In-place update of the numeric values in the CHOLMOD matrix.
"""
function set_values!(mat::FixedStructureCHOLMOD{T, I}, new_vals::AbstractVector{T}) where {T<:VTypes, I<:ITypes}
@assert length(new_vals) == length(mat._values) "New values must be same length as allocated storage!"
copyto!(mat._values, new_vals)
return mat
end
function symbolic_factor(mat::FixedStructureCHOLMOD{T, I}) where {T<:VTypes, I<:ITypes}
return SparseArrays.CHOLMOD.symbolic(mat._mat)
end
function numeric_factor!(F::SparseArrays.CHOLMOD.Factor{T, I}, mat::FixedStructureCHOLMOD{T, I}) where {T<:VTypes, I<:ITypes}
SparseArrays.CHOLMOD.factorize!(mat._mat, F)
return
end
```
My question: should this struct be mutable, with a proper finalizer, or is it memory-safe as is? edit: I’m asking this in the context described above. So constant sparse structure, one call to `symbolic_factor`, then many calls to `numeric_refactor!`.
Useful links: SuiteSparse’s [cholmod.h](https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/7e737e6aeca7b0ef20f8328fdf004ebab6ff4802/CHOLMOD/Include/cholmod.h#L1168), CHOLMOD’s [User Guide](https://fossies.org/linux/SuiteSparse/CHOLMOD/Doc/CHOLMOD_UserGuide.pdf), and Julia’s [cholmod.jl](https://github.com/JuliaSparse/SparseArrays.jl/blob/30201abcb41e558a4b5cc23b11dc5676c5655c0b/src/solvers/cholmod.jl).
> **Sample test code**
>
> ```julia
> const K = 10
> A = sprand(10, 10, 0.3)
> B = A + A' + K * LinearAlgebra.I(10)
> test_struct = FixedStructureCHOLMOD(B)
> sym = symbolic_factor(test_struct)
> numeric_factor!(sym, test_struct)
> b = rand(10)
> x = sym \ b
> @assert isapprox(B * x, b, atol=1e-6)
>
> # Test the in-place modification of values.
> A2 = deepcopy(A)
> nonzeros(A2) = rand(Float64, size(nonzeros(A2)))
> B2 = A2 + A2' + K * LinearAlgebra.I(10) # same structure
> set_values!(test_struct, SparseArrays.nonzeros(B2))
> numeric_factor!(sym, test_struct)
> x2 = sym \ b
> @assert isapprox(B2 * x2, b, atol=1e-6)
>
> ```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 22, 2025, 1:23am UTC](https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922/2 "2025-07-22T01:23:56Z")
Why not just do `F = cholesky(B)` the first time and then do `cholesky!(F, B)` thereafter?
This is the documented way to re-use the symbolic factorization, and then you can just modify `B` and `F` in-place instead of messing with low-level `CHOLMOD` data structures.
---
### Author:  [@luke-kiernan](https://discourse.julialang.org/u/luke-kiernan)
#### Post date: [July 22, 2025, 2:49am UTC](https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922/3 "2025-07-22T02:49:24Z")
I expect that the conversion to `Sparse` [here](https://github.com/JuliaSparse/SparseArrays.jl/blob/30201abcb41e558a4b5cc23b11dc5676c5655c0b/src/solvers/cholmod.jl#L1520C1-L1526C73) is allocating and has a negative impact on performance. I wrote up a quick-and-dirty script designed to test that assumption:
```julia
using LinearAlgebra
using SparseArrays
using SuiteSparse
function testSolve(F, A)
b = rand(size(A, 1))
x = F \ b
@assert isapprox(A * x, b, atol = 1e-6)
end
function choleskyAllocations(A1, A2)
F = cholesky(A1)
testSolve(F, A1)
result = @allocated cholesky!(F, A2)
testSolve(F, A2)
return result
end
for N in (100, 1000, 10000)
density = 0.001
B = sprand(N, N, density)
A1 = B + B' + 2*N/100 * LinearAlgebra.I(N)
copyto!(nonzeros(B), rand(size(nonzeros(B), 1)))
A2 = B + B' + N/100 * LinearAlgebra.I(N)
println("N = $N:")
println(" cholesky! memory usage: $(choleskyAllocations(A1, A2))")
println(" Sparse memory usage: $(@allocated SuiteSparse.CHOLMOD.Sparse(A2))")
end
```
Output:
```julia
N = 100:
cholesky! memory usage: 5936
Sparse memory usage: 4640
N = 1000:
cholesky! memory usage: 113152
Sparse memory usage: 57024
N = 10000:
cholesky! memory usage: 250266272
Sparse memory usage: 3506880
```
Huh. So for large matrices, converting to `Sparse` has a pretty insignificant memory footprint compared to the library call itself, from half the memory footprint at N = 1,000 to around 1% at N = 10,000.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 22, 2025, 3:03am UTC](https://discourse.julialang.org/t/fixed-structure-cholmod-with-in-place-set-nonzeros/130922/4 "2025-07-22T03:03:22Z")
> [@luke-kiernan](#):
>
> I expect that the conversion to `Sparse` [here](https://github.com/JuliaSparse/SparseArrays.jl/blob/30201abcb41e558a4b5cc23b11dc5676c5655c0b/src/solvers/cholmod.jl#L1520C1-L1526C73) is allocating
Yes, I think you are right — looks like it makes a copy of the data in the constructor [here](https://github.com/JuliaSparse/SparseArrays.jl/blob/30201abcb41e558a4b5cc23b11dc5676c5655c0b/src/solvers/cholmod.jl#L1018-L1042).
> [@luke-kiernan](#):
>
> Huh. So for large matrices, converting to `Sparse` has a pretty insignificant memory footprint compared to the library call itself, from half the memory footprint at N = 1,000 to around 1% at N = 10,000.
It’s certainly true that the sparse Cholesky factor typically requires a lot more memory than the matrix, due to fill-in. However, be aware that this depends strongly on the sparsity pattern, so using `sprand` may be misleading if your actual matrix is more structured.
(I’m also not sure if Julia is tracking all of the allocations that happen in the external library here?)