Skip to content

Commit 6fe51f2

Browse files
danielwedanielwe
authored
Improve CHOLMOD friendliness (#645)
1 parent bb5ecc0 commit 6fe51f2

File tree

2 files changed

+53
-17
lines changed

2 files changed

+53
-17
lines changed

‎src/solvers/cholmod.jl

Lines changed: 40 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -10,8 +10,8 @@
1010

1111
module CHOLMOD
1212

13-
import Base: (*), convert, copy, eltype, getindex, getproperty, show, size,
14-
IndexStyle, IndexLinear, IndexCartesian, adjoint, axes,
13+
import Base: (*), convert, copy, eltype, getindex, getproperty, propertynames,
14+
show, size, IndexStyle, IndexLinear, IndexCartesian, adjoint, axes,
1515
Matrix, Vector
1616
using Base: require_one_based_indexing
1717

@@ -380,24 +380,25 @@ Base.pointer(x::Dense{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_dense},
380380
Base.pointer(x::Sparse{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_sparse}, x)
381381
Base.pointer(x::Factor{Tv}) where {Tv} = Base.unsafe_convert(Ptr{cholmod_factor}, x)
382382

383+
function _factor_components(is_ll)
384+
is_ll ? (:L, :U, :PtL, :UP) : (:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP)
385+
end
386+
383387
# FactorComponent, for encoding particular factors from a factorization
384388
mutable struct FactorComponent{Tv, S, Ti} <: AbstractMatrix{Tv}
385389
F::Factor{Tv, Ti}
386390

387-
function FactorComponent{Tv, S, Ti}(F::Factor{Tv, Ti}) where {Tv, S, Ti}
388-
s = unsafe_load(pointer(F))
389-
if s.is_ll != 0
390-
if !(S === :L || S === :U || S === :PtL || S === :UP)
391-
throw(CHOLMODException(string(S, " not supported for sparse ",
392-
"LLt matrices; try :L, :U, :PtL, or :UP")))
393-
end
394-
elseif !(S === :L || S === :U || S === :PtL || S === :UP ||
395-
S === :D || S === :LD || S === :DU || S === :PtLD || S === :DUP)
396-
throw(CHOLMODException(string(S, " not supported for sparse LDLt ",
397-
"matrices; try :L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, or :DUP")))
398-
end
399-
new(F)
400-
end
391+
function FactorComponent{Tv,S,Ti}(F::Factor{Tv,Ti}) where {Tv,S,Ti}
392+
s = unsafe_load(pointer(F))
393+
is_ll = (s.is_ll != 0)
394+
components = _factor_components(is_ll)
395+
if !(S in components)
396+
type = is_ll ? "LLt" : "LDLt"
397+
component_list = join(repr.(components), ", ")
398+
throw(CHOLMODException("$(repr(S)) not supported for sparse $type factorizations; try $component_list, or :p"))
399+
end
400+
new(F)
401+
end
401402
end
402403
function FactorComponent{Tv, S}(F::Factor{Tv, Ti}) where {Tv, S, Ti}
403404
return FactorComponent{Tv, S, Ti}(F)
@@ -1285,11 +1286,17 @@ function sparse(FC::FactorComponent{Tv,:L}) where Tv
12851286
F = Factor(FC)
12861287
s = unsafe_load(pointer(F))
12871288
if s.is_ll == 0
1288-
throw(CHOLMODException("sparse: supported only for :LD on LDLt factorizations"))
1289+
_sparse_exception(F)
12891290
end
12901291
sparse(Sparse(F))
12911292
end
12921293
sparse(FC::FactorComponent{Tv,:LD}) where {Tv} = sparse(Sparse(Factor(FC)))
1294+
sparse(FC::FactorComponent{Tv}) where {Tv} = _sparse_exception(Factor(FC))
1295+
function _sparse_exception(F::Factor)
1296+
s = unsafe_load(pointer(F))
1297+
details = (s.is_ll == 0) ? ":LD on LDLt" : ":L on LLt"
1298+
throw(CHOLMODException("sparse: supported only for $details factorizations"))
1299+
end
12931300

12941301
# Calculate the offset into the stype field of the cholmod_sparse_struct and
12951302
# change the value
@@ -1396,6 +1403,11 @@ end
13961403
end
13971404
end
13981405

1406+
function propertynames(F::Factor)
1407+
s = unsafe_load(pointer(F))
1408+
(_factor_components(s.is_ll != 0)..., :p, :ptr)
1409+
end
1410+
13991411
function getLd!(S::SparseMatrixCSC)
14001412
d = Vector{eltype(S)}(undef, size(S, 1))
14011413
fill!(d, 0)
@@ -1555,6 +1567,12 @@ using just `L` without accounting for `P` will give incorrect answers.
15551567
To include the effects of permutation,
15561568
it's typically preferable to extract "combined" factors like `PtL = F.PtL`
15571569
(the equivalent of `P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`).
1570+
The complete list of supported factors is `:L, :PtL, :UP, :U`.
1571+
The permutation vector is available as `F.p`, defined such that `L*L' == A[p, p]`,
1572+
1573+
The `L` component can be materialized as a sparse matrix using `sparse(F.L)`.
1574+
Other components cannot be materialized directly, but can be reconstructed
1575+
from `sparse(F.L)` and `F.p` if needed.
15581576
15591577
When `check = true`, an error is thrown if the decomposition fails.
15601578
When `check = false`, responsibility for checking the decomposition's
@@ -1732,6 +1750,11 @@ To include the effects of permutation, it is typically preferable to extract
17321750
"combined" factors like `PtL = F.PtL` (the equivalent of
17331751
`P'*L`) and `LtP = F.UP` (the equivalent of `L'*P`).
17341752
The complete list of supported factors is `:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP`.
1753+
The permutation vector is available as `F.p`, defined such that `L*D*L' == A[p, p]`,
1754+
1755+
The `LD` component can be materialized as a sparse matrix using `sparse(F.LD)`,
1756+
Other components cannot be materialized directly, but can be reconstructed from
1757+
`sparse(F.LD)` and `F.p` if needed.
17351758
17361759
Unlike the related Cholesky factorization, the ``LDL'`` factorization does not
17371760
require `A` to be positive definite. However, it still requires all leading

‎test/cholmod.jl

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -589,9 +589,13 @@ end
589589

590590
@testset "cholesky, no permutation $Tv" begin
591591
Fs = cholesky(As, perm=[1:3;])
592+
@test sort(collect(propertynames(Fs))) == sort([:L, :U, :PtL, :UP, :p, :ptr])
592593
@test Fs.p == [1:3;]
593594
@test sparse(Fs.L) ≈ Lf
594595
@test sparse(Fs) ≈ As
596+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.U)
597+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.PtL)
598+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :L on LLt factorizations") sparse(Fs.UP)
595599
b = rand(Tv, 3)
596600
bs = sparse(b)
597601
@test Fs\b ≈ Af\b ≈ (Fs\bs)::SparseVector
@@ -645,9 +649,18 @@ end
645649

646650
@testset "ldlt, no permutation" begin
647651
Fs = ldlt(As, perm=[1:3;])
652+
@test sort(collect(propertynames(Fs))) == sort([:L, :U, :PtL, :UP, :D, :LD, :DU, :PtLD, :DUP, :p, :ptr])
648653
@test Fs.p == [1:3;]
649654
@test sparse(Fs.LD) ≈ LDf
650655
@test sparse(Fs) ≈ As
656+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.L)
657+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.U)
658+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.PtL)
659+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.UP)
660+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.D)
661+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.DU)
662+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.PtLD)
663+
@test_throws CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations") sparse(Fs.DUP)
651664
b = rand(Tv, 3)
652665
bs = sparse(b)
653666
@test Fs\b ≈ Af\b ≈ (Fs\bs)::SparseVector

0 commit comments

Comments
 (0)