From 7df519c62622f04709c2c9d6de0785f708fe6c36 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sun, 7 Jan 2024 15:47:16 -0500 Subject: [PATCH] UMFPACK and KLU both throw errors with singular matrices that are not especially easy to silence --- src/default.jl | 3 ++- src/factorization.jl | 8 +++++--- test/default_algs.jl | 16 ++++++++++------ 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/default.jl b/src/default.jl index 8972b2a38..6317daf48 100644 --- a/src/default.jl +++ b/src/default.jl @@ -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 + if false && assump.issq if length(b) <= 10_000 && length(nonzeros(A)) / length(A) < 2e-4 DefaultLinearSolver(DefaultAlgorithmChoice.KLUFactorization) else diff --git a/src/factorization.jl b/src/factorization.jl index 9d9b64804..ffd12eb7f 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -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) end cache.cacheval = fact cache.isfresh = false end + # TODO gaurd this against SingularExceptions y = ldiv!(cache.u, @get_cacheval(cache, :UMFPACKFactorization), cache.b) SciMLBase.build_linear_solution(alg, y, nothing, cache) end @@ -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) diff --git a/test/default_algs.jl b/test/default_algs.jl index 5aad1a30c..4d7a766ca 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -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 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)