Skip to content

Improve CHOLMOD friendliness #645

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.

Sign up for GitHub

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

Merged
merged 4 commits into from
Aug 29, 2025
Merged

Improve CHOLMOD friendliness #645

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f5acf29
Implement propertynames(::Factor)
danielwe Aug 27, 2025
bcbb1cc
Add more informative `sparse` exceptions
danielwe Aug 27, 2025
8a7bdfe
Add CHOLMOD docs about materializing; list all components
danielwe Aug 27, 2025
bd3536e
Add tests
danielwe Aug 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 40 additions & 17 deletions src/solvers/cholmod.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@

module CHOLMOD

import Base: (*), convert, copy, eltype, getindex, getproperty, show, size,
IndexStyle, IndexLinear, IndexCartesian, adjoint, axes,
import Base: (*), convert, copy, eltype, getindex, getproperty, propertynames,
show, size, IndexStyle, IndexLinear, IndexCartesian, adjoint, axes,
Matrix, Vector
using Base: require_one_based_indexing

Expand Down Expand Up @@ -380,24 +380,25 @@ 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)

function _factor_components(is_ll)
is_ll ? (:L, :U, :PtL, :UP) : (:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP)
end

# FactorComponent, for encoding particular factors from a factorization
mutable struct FactorComponent{Tv, S, Ti} <: AbstractMatrix{Tv}
F::Factor{Tv, Ti}

function FactorComponent{Tv, S, Ti}(F::Factor{Tv, Ti}) where {Tv, S, Ti}
s = unsafe_load(pointer(F))
if s.is_ll != 0
if !(S === :L || S === :U || S === :PtL || S === :UP)
throw(CHOLMODException(string(S, " not supported for sparse ",
"LLt matrices; try :L, :U, :PtL, or :UP")))
end
elseif !(S === :L || S === :U || S === :PtL || S === :UP ||
S === :D || S === :LD || S === :DU || S === :PtLD || S === :DUP)
throw(CHOLMODException(string(S, " not supported for sparse LDLt ",
"matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")))
end
new(F)
end
function FactorComponent{Tv,S,Ti}(F::Factor{Tv,Ti}) where {Tv,S,Ti}
s = unsafe_load(pointer(F))
is_ll = (s.is_ll != 0)
components = _factor_components(is_ll)
if !(S in components)
type = is_ll ? "LLt" : "LDLt"
component_list = join(repr.(components), ", ")
throw(CHOLMODException("$(repr(S)) not supported for sparse $type factorizations; try $component_list, or :p"))
end
new(F)
end
end
function FactorComponent{Tv, S}(F::Factor{Tv, Ti}) where {Tv, S, Ti}
return FactorComponent{Tv, S, Ti}(F)
Expand Down Expand Up @@ -1285,11 +1286,17 @@ function sparse(FC::FactorComponent{Tv,:L}) where Tv
F = Factor(FC)
s = unsafe_load(pointer(F))
if s.is_ll == 0
throw(CHOLMODException("sparse: supported only for :LD on LDLt factorizations"))
_sparse_exception(F)
end
sparse(Sparse(F))
end
sparse(FC::FactorComponent{Tv,:LD}) where {Tv} = sparse(Sparse(Factor(FC)))
sparse(FC::FactorComponent{Tv}) where {Tv} = _sparse_exception(Factor(FC))
function _sparse_exception(F::Factor)
s = unsafe_load(pointer(F))
details = (s.is_ll == 0) ? ":LD on LDLt" : ":L on LLt"
throw(CHOLMODException("sparse: supported only for $details factorizations"))
end

# Calculate the offset into the stype field of the cholmod_sparse_struct and
# change the value
Expand Down Expand Up @@ -1396,6 +1403,11 @@ end
end
end

function propertynames(F::Factor)
s = unsafe_load(pointer(F))
(_factor_components(s.is_ll != 0)..., :p, :ptr)
end

function getLd!(S::SparseMatrixCSC)
d = Vector{eltype(S)}(undef, size(S, 1))
fill!(d, 0)
Expand Down Expand Up @@ -1555,6 +1567,12 @@ using just `L` without accounting for `P` will give incorrect answers.
To include the effects of permutation,
it's typically preferable to extract "combined" factors like `PtL = F.PtL`
(the equivalent of `P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`).
The complete list of supported factors is `:L, :PtL, :UP, :U`.
The permutation vector is available as `F.p`, defined such that `L*L' == A[p, p]`,

The `L` component can be materialized as a sparse matrix using `sparse(F.L)`.
Other components cannot be materialized directly, but can be reconstructed
from `sparse(F.L)` and `F.p` if needed.

When `check = true`, an error is thrown if the decomposition fails.
When `check = false`, responsibility for checking the decomposition's
Expand Down Expand Up @@ -1732,6 +1750,11 @@ To include the effects of permutation, it is typically preferable to extract
"combined" factors like `PtL = F.PtL` (the equivalent of
`P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`).
The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
The permutation vector is available as `F.p`, defined such that `L*D*L' == A[p, p]`,

The `LD` component can be materialized as a sparse matrix using `sparse(F.LD)`,
Other components cannot be materialized directly, but can be reconstructed from
`sparse(F.LD)` and `F.p` if needed.

Unlike the related Cholesky factorization, the ``LDL'`` factorization does not
require `A` to be positive definite. However, it still requires all leading
Expand Down
13 changes: 13 additions & 0 deletions test/cholmod.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -589,9 +589,13 @@ end

@testset "cholesky, no permutation $Tv" begin
Fs = cholesky(As, perm=[1:3;])
@test sort(collect(propertynames(Fs))) == sort([:L, :U, :PtL, :UP, :p, :ptr])
@test Fs.p == [1:3;]
@test sparse(Fs.L) ≈ Lf
@test sparse(Fs) ≈ As
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.U)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.PtL)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.UP)
b = rand(Tv, 3)
bs = sparse(b)
@test Fs\b ≈ Af\b ≈ (Fs\bs)::SparseVector
Expand Down Expand Up @@ -645,9 +649,18 @@ end

@testset "ldlt, no permutation" begin
Fs = ldlt(As, perm=[1:3;])
@test sort(collect(propertynames(Fs))) == sort([:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP, :p, :ptr])
@test Fs.p == [1:3;]
@test sparse(Fs.LD) ≈ LDf
@test sparse(Fs) ≈ As
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.L)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.U)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.PtL)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.UP)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.D)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.DU)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.PtLD)
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.DUP)
b = rand(Tv, 3)
bs = sparse(b)
@test Fs\b ≈ Af\b ≈ (Fs\bs)::SparseVector
Expand Down
Loading