Skip to content

Commit

Permalink
Make cholmod preserve user-specified permutation
Browse files Browse the repository at this point in the history
  • Loading branch information
Simon Etter committed Apr 22, 2021
1 parent 7fe05db commit 0bef415
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 4 deletions.
2 changes: 2 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,8 @@ Standard library changes

#### SparseArrays

* `cholesky()` now fully preserves the user-specified permutation. ([#40560])


#### Dates

Expand Down
8 changes: 4 additions & 4 deletions stdlib/SuiteSparse/src/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1374,9 +1374,9 @@ end
## Factorization methods

## Compute that symbolic factorization only
function fact_(A::Sparse{<:VTypes};
function symbolic(A::Sparse{<:VTypes};
perm::Union{Nothing,AbstractVector{SuiteSparse_long}}=nothing,
postorder::Bool=true, userperm_only::Bool=true)
postorder::Bool=isnothing(perm)||isempty(perm), userperm_only::Bool=true)

sA = unsafe_load(pointer(A))
sA.stype == 0 && throw(ArgumentError("sparse matrix is not symmetric/Hermitian"))
Expand Down Expand Up @@ -1438,7 +1438,7 @@ function cholesky(A::Sparse; shift::Real=0.0, check::Bool = true,
cm.print = 0

# Compute the symbolic factorization
F = fact_(A; perm = perm)
F = symbolic(A; perm = perm)

# Compute the numerical factorization
cholesky!(F, A; shift = shift, check = check)
Expand Down Expand Up @@ -1607,7 +1607,7 @@ function ldlt(A::Sparse; shift::Real=0.0, check::Bool = true,
cm.supernodal = 0

# Compute the symbolic factorization
F = fact_(A; perm = perm)
F = symbolic(A; perm = perm)

# Compute the numerical factorization
ldlt!(F, A; shift = shift, check = check)
Expand Down
17 changes: 17 additions & 0 deletions stdlib/SuiteSparse/test/cholmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -909,3 +909,20 @@ end
@test C * C' == Sparse(spzeros(3, 3))
@test C' * C == Sparse(spzeros(0, 0))
end

@testset "permutation handling" begin
@testset "default permutation" begin
# Assemble arrow matrix
A = sparse(5I,3,3)
A[:,1] .= 1; A[1,:] .= A[:,1]

# Ensure cholesky eliminates the fill-in
@test cholesky(A).p[1] != 1
end

@testset "user-specified permutation" begin
n = 100
A = sprand(n,n,5/n) |> t -> t't + I
@test cholesky(A, perm=1:n).p == 1:n
end
end

0 comments on commit 0bef415

Please sign in to comment.