Skip to content

Commit

Permalink
UMFPACK and KLU both throw errors with singular matrices that are not…
Browse files Browse the repository at this point in the history
… especially easy to silence
  • Loading branch information
oscardssmith committed Jan 7, 2024
1 parent 49183f7 commit 7df519c
Show file tree
Hide file tree
Showing 3 changed files with 17 additions and 10 deletions.
3 changes: 2 additions & 1 deletion src/default.jl
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,8 @@ end
@static if INCLUDE_SPARSE
function defaultalg(A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b,
assump::OperatorAssumptions{Bool}) where {Ti}
if assump.issq
#TODO: find a way to supress the errors that KLU and UMFPACK.jl throw for singluar matrices

Check warning on line 94 in src/default.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"supress" should be "suppress".

Check warning on line 94 in src/default.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"singluar" should be "singular".
if false && assump.issq
if length(b) <= 10_000 && length(nonzeros(A)) / length(A) < 2e-4
DefaultLinearSolver(DefaultAlgorithmChoice.KLUFactorization)
else
Expand Down
8 changes: 5 additions & 3 deletions src/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -748,19 +748,20 @@ function SciMLBase.solve!(cache::LinearCache, alg::UMFPACKFactorization; kwargs.
SparseArrays.decrement(SparseArrays.getrowval(A)) ==
cacheval.rowval)
fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
nonzeros(A)))
nonzeros(A)), check=false)
else
fact = lu!(cacheval,
SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
nonzeros(A)))
nonzeros(A)), check=false)
end
else
fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)))
fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), check=false)

Check warning on line 758 in src/factorization.jl

View check run for this annotation

Codecov / codecov/patch

src/factorization.jl#L758

Added line #L758 was not covered by tests
end
cache.cacheval = fact
cache.isfresh = false
end

# TODO gaurd this against SingularExceptions

Check warning on line 764 in src/factorization.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"gaurd" should be "guard" or "gourd".
y = ldiv!(cache.u, @get_cacheval(cache, :UMFPACKFactorization), cache.b)
SciMLBase.build_linear_solution(alg, y, nothing, cache)
end
Expand Down Expand Up @@ -808,6 +809,7 @@ function init_cacheval(alg::KLUFactorization, A::AbstractSparseArray, b, u, Pl,
nonzeros(A)))
end

# TODO: guard this against errors
function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...)
A = cache.A
A = convert(AbstractMatrix, A)
Expand Down
16 changes: 10 additions & 6 deletions test/default_algs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,19 @@ solve(prob)
LinearSolve.OperatorAssumptions(false)).alg ===
LinearSolve.DefaultAlgorithmChoice.QRFactorization

@test LinearSolve.defaultalg(sprand(10^4, 10^4, 1e-5) + I, zeros(1000)).alg ===
LinearSolve.DefaultAlgorithmChoice.KLUFactorization
prob = LinearProblem(sprand(1000, 1000, 0.5), zeros(1000))

A = spzeros(100, 100)
A[1,1]=1
prob = LinearProblem(A, ones(100))
# test that solving a singluar problem doesn't error

Check warning on line 42 in test/default_algs.jl

View workflow job for this annotation

GitHub Actions / Spell Check with Typos

"singluar" should be "singular".
solve(prob)

@test LinearSolve.defaultalg(sprand(11000, 11000, 0.001), zeros(11000)).alg ===
# test_broken because these would be faster, but are currently not used by default
# because they solvers throw errors for singular problems
@test_broken LinearSolve.defaultalg(sprand(10^4, 10^4, 1e-5) + I, zeros(1000)).alg ===
LinearSolve.DefaultAlgorithmChoice.KLUFactorization
@test_broken LinearSolve.defaultalg(sprand(11000, 11000, 0.001), zeros(11000)).alg ===
LinearSolve.DefaultAlgorithmChoice.UMFPACKFactorization
prob = LinearProblem(sprand(11000, 11000, 0.5), zeros(11000))
solve(prob)

# Test inference
A = rand(4, 4)
Expand Down

0 comments on commit 7df519c

Please sign in to comment.