From 5131e03f5d203c49b418c0a2683922b747899449 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 24 Sep 2019 08:33:29 -0400 Subject: [PATCH 1/4] =?UTF-8?q?fix=20exception=20type=20for=20sparse=20LDL?= =?UTF-8?q?=E1=B5=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stdlib/SuiteSparse/src/cholmod.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/SuiteSparse/src/cholmod.jl b/stdlib/SuiteSparse/src/cholmod.jl index 20c24e7c7e2d1..9af38d41ae57b 100644 --- a/stdlib/SuiteSparse/src/cholmod.jl +++ b/stdlib/SuiteSparse/src/cholmod.jl @@ -744,7 +744,7 @@ function solve(sys::Integer, F::Factor{Tv}, B::Dense{Tv}) where Tv<:VTypes if !issuccess(F) s = unsafe_load(pointer(F)) if s.is_ll == 1 - throw(LinearAlgebra.PosDefException(s.minor)) + throw(LinearAlgebra.SingularException(s.minor)) else throw(ArgumentError("factorized matrix has one or more zero pivots. Try using `lu` instead.")) end @@ -1446,7 +1446,7 @@ function ldlt!(F::Factor{Tv}, A::Sparse{Tv}; # Compute the numerical factorization factorize_p!(A, shift, F, cm) - check && (issuccess(F) || throw(LinearAlgebra.PosDefException(1))) + check && (issuccess(F) || throw(LinearAlgebra.SingularException(1))) return F end From 768d014dca270967b55b0ae649d0c4b242283f3d Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 24 Sep 2019 18:48:47 -0400 Subject: [PATCH 2/4] update tests for SingularException --- stdlib/SuiteSparse/test/cholmod.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/stdlib/SuiteSparse/test/cholmod.jl b/stdlib/SuiteSparse/test/cholmod.jl index a19fafd3e9f99..3a6eab67973db 100644 --- a/stdlib/SuiteSparse/test/cholmod.jl +++ b/stdlib/SuiteSparse/test/cholmod.jl @@ -5,7 +5,7 @@ using DelimitedFiles using Test using Random using Serialization -using LinearAlgebra: issuccess, PosDefException +using LinearAlgebra: issuccess, PosDefException, SingularException using SparseArrays using SparseArrays: getcolptr @@ -384,8 +384,8 @@ end b = fill(1., size(A1, 1)) @test_throws PosDefException cholesky(C - 2λmaxC*I) @test_throws PosDefException cholesky(C, shift=-2λmaxC) - @test_throws PosDefException ldlt(C - C[1,1]*I) - @test_throws PosDefException ldlt(C, shift=-real(C[1,1])) + @test_throws SingularException ldlt(C - C[1,1]*I) + @test_throws SingularException ldlt(C, shift=-real(C[1,1])) @test !isposdef(cholesky(C - 2λmaxC*I; check = false)) @test !isposdef(cholesky(C, shift=-2λmaxC; check = false)) @test !issuccess(ldlt(C - C[1,1]*I; check = false)) @@ -849,8 +849,8 @@ end B = sparse(ComplexF64[0 0; 0 0]) for M in (A, B, Symmetric(A), Hermitian(B)) F = ldlt(M; check = false) - @test_throws PosDefException ldlt(M) - @test_throws PosDefException ldlt!(F, M) + @test_throws SingularException ldlt(M) + @test_throws SingularException ldlt!(F, M) @test !issuccess(ldlt(M; check = false)) @test !issuccess(ldlt!(F, M; check = false)) end From 5037f874441419d561d53b9a87339373dc1559c8 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 25 Sep 2019 09:01:22 -0400 Subject: [PATCH 3/4] whoops, is_ll means this is Cholesky --- stdlib/SuiteSparse/src/cholmod.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/SuiteSparse/src/cholmod.jl b/stdlib/SuiteSparse/src/cholmod.jl index 9af38d41ae57b..9a4aed6a8b4fe 100644 --- a/stdlib/SuiteSparse/src/cholmod.jl +++ b/stdlib/SuiteSparse/src/cholmod.jl @@ -744,7 +744,7 @@ function solve(sys::Integer, F::Factor{Tv}, B::Dense{Tv}) where Tv<:VTypes if !issuccess(F) s = unsafe_load(pointer(F)) if s.is_ll == 1 - throw(LinearAlgebra.SingularException(s.minor)) + throw(LinearAlgebra.PosDefException(s.minor)) else throw(ArgumentError("factorized matrix has one or more zero pivots. Try using `lu` instead.")) end From 61c1af74d6bdc449ce9fe5ecb44504c80dc3b5d4 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 27 Sep 2019 11:39:20 -0400 Subject: [PATCH 4/4] ZeroPivotException --- NEWS.md | 2 ++ stdlib/LinearAlgebra/docs/src/index.md | 1 + stdlib/LinearAlgebra/src/exceptions.jl | 19 ++++++++++++++++++- stdlib/LinearAlgebra/src/factorization.jl | 4 +++- stdlib/LinearAlgebra/src/lu.jl | 4 ++-- stdlib/LinearAlgebra/test/lu.jl | 8 ++++---- stdlib/SuiteSparse/src/cholmod.jl | 4 ++-- stdlib/SuiteSparse/test/cholmod.jl | 10 +++++----- 8 files changed, 37 insertions(+), 15 deletions(-) diff --git a/NEWS.md b/NEWS.md index 661e4b85bfcda..fb648764f2292 100644 --- a/NEWS.md +++ b/NEWS.md @@ -46,6 +46,8 @@ Standard library changes * `dot` now admits a 3-argument method `dot(x, A, y)` to compute generalized dot products `dot(x, A*y)`, but without computing and storing the intermediate result `A*y` ([#32739]). +* `ldlt` and non-pivoted `lu` now throw a new `ZeroPivotException` type ([#33372]). + #### Random diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index df83cc81cedd1..d341eb7f22e84 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -318,6 +318,7 @@ Base.:*(::AbstractMatrix, ::AbstractMatrix) Base.:\(::AbstractMatrix, ::AbstractVecOrMat) LinearAlgebra.SingularException LinearAlgebra.PosDefException +LinearAlgebra.ZeroPivotException LinearAlgebra.dot LinearAlgebra.cross LinearAlgebra.factorize diff --git a/stdlib/LinearAlgebra/src/exceptions.jl b/stdlib/LinearAlgebra/src/exceptions.jl index 64aa0f6c7d34f..6704a9ac6ae4d 100644 --- a/stdlib/LinearAlgebra/src/exceptions.jl +++ b/stdlib/LinearAlgebra/src/exceptions.jl @@ -3,7 +3,8 @@ export LAPACKException, SingularException, PosDefException, - RankDeficientException + RankDeficientException, + ZeroPivotException struct LAPACKException <: Exception info::BlasInt @@ -43,3 +44,19 @@ end struct RankDeficientException <: Exception info::BlasInt end + +""" + ZeroPivotException <: Exception + +Exception thrown when a matrix factorization/solve encounters a zero in a pivot (diagonal) +position and cannot proceed. This may *not* mean that the matrix is singular: +it may be fruitful to switch to a diffent factorization such as pivoted LU +that can re-order variables to eliminate spurious zero pivots. +The `info` field indicates the location of (one of) the zero pivot(s). +""" +struct ZeroPivotException <: Exception + info::BlasInt +end +function Base.showerror(io::IO, ex::ZeroPivotException) + print(io, "ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.") +end \ No newline at end of file diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 848dbdc09aa86..e90bffedb547a 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -16,7 +16,9 @@ size(F::Adjoint{<:Any,<:Factorization}) = reverse(size(parent(F))) size(F::Transpose{<:Any,<:Factorization}) = reverse(size(parent(F))) checkpositivedefinite(info) = info == 0 || throw(PosDefException(info)) -checknonsingular(info) = info == 0 || throw(SingularException(info)) +checknonsingular(info, pivoted::Val{true}) = info == 0 || throw(SingularException(info)) +checknonsingular(info, pivoted::Val{false}) = info == 0 || throw(ZeroPivotException(info)) +checknonsingular(info) = checknonsingular(info, Val{true}()) """ issuccess(F::Factorization) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index fcca98ab004f4..c82eae2cbb49b 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -175,7 +175,7 @@ function generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot} = Val(true); end end end - check && checknonsingular(info) + check && checknonsingular(info, Val{Pivot}()) return LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info)) end @@ -552,7 +552,7 @@ function lu!(A::Tridiagonal{T,V}, pivot::Union{Val{false}, Val{true}} = Val(true end end B = Tridiagonal{T,V}(dl, d, du, du2) - check && checknonsingular(info) + check && checknonsingular(info, pivot) return LU{T,Tridiagonal{T,V}}(B, ipiv, convert(BlasInt, info)) end diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 5833c59c8dda8..28f26367e5141 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -207,10 +207,10 @@ end @test_throws SingularException lu!(copy(A); check = true) @test !issuccess(lu(A; check = false)) @test !issuccess(lu!(copy(A); check = false)) - @test_throws SingularException lu(A, Val(false)) - @test_throws SingularException lu!(copy(A), Val(false)) - @test_throws SingularException lu(A, Val(false); check = true) - @test_throws SingularException lu!(copy(A), Val(false); check = true) + @test_throws ZeroPivotException lu(A, Val(false)) + @test_throws ZeroPivotException lu!(copy(A), Val(false)) + @test_throws ZeroPivotException lu(A, Val(false); check = true) + @test_throws ZeroPivotException lu!(copy(A), Val(false); check = true) @test !issuccess(lu(A, Val(false); check = false)) @test !issuccess(lu!(copy(A), Val(false); check = false)) F = lu(A; check = false) diff --git a/stdlib/SuiteSparse/src/cholmod.jl b/stdlib/SuiteSparse/src/cholmod.jl index 9a4aed6a8b4fe..86381eeb3d179 100644 --- a/stdlib/SuiteSparse/src/cholmod.jl +++ b/stdlib/SuiteSparse/src/cholmod.jl @@ -746,7 +746,7 @@ function solve(sys::Integer, F::Factor{Tv}, B::Dense{Tv}) where Tv<:VTypes if s.is_ll == 1 throw(LinearAlgebra.PosDefException(s.minor)) else - throw(ArgumentError("factorized matrix has one or more zero pivots. Try using `lu` instead.")) + throw(LinearAlgebra.ZeroPivotException(s.minor)) end end Dense(ccall((@cholmod_name("solve"),:libcholmod), Ptr{C_Dense{Tv}}, @@ -1446,7 +1446,7 @@ function ldlt!(F::Factor{Tv}, A::Sparse{Tv}; # Compute the numerical factorization factorize_p!(A, shift, F, cm) - check && (issuccess(F) || throw(LinearAlgebra.SingularException(1))) + check && (issuccess(F) || throw(LinearAlgebra.ZeroPivotException(1))) return F end diff --git a/stdlib/SuiteSparse/test/cholmod.jl b/stdlib/SuiteSparse/test/cholmod.jl index 3a6eab67973db..4f6927b14f05a 100644 --- a/stdlib/SuiteSparse/test/cholmod.jl +++ b/stdlib/SuiteSparse/test/cholmod.jl @@ -5,7 +5,7 @@ using DelimitedFiles using Test using Random using Serialization -using LinearAlgebra: issuccess, PosDefException, SingularException +using LinearAlgebra: issuccess, PosDefException, ZeroPivotException using SparseArrays using SparseArrays: getcolptr @@ -384,8 +384,8 @@ end b = fill(1., size(A1, 1)) @test_throws PosDefException cholesky(C - 2λmaxC*I) @test_throws PosDefException cholesky(C, shift=-2λmaxC) - @test_throws SingularException ldlt(C - C[1,1]*I) - @test_throws SingularException ldlt(C, shift=-real(C[1,1])) + @test_throws ZeroPivotException ldlt(C - C[1,1]*I) + @test_throws ZeroPivotException ldlt(C, shift=-real(C[1,1])) @test !isposdef(cholesky(C - 2λmaxC*I; check = false)) @test !isposdef(cholesky(C, shift=-2λmaxC; check = false)) @test !issuccess(ldlt(C - C[1,1]*I; check = false)) @@ -849,8 +849,8 @@ end B = sparse(ComplexF64[0 0; 0 0]) for M in (A, B, Symmetric(A), Hermitian(B)) F = ldlt(M; check = false) - @test_throws SingularException ldlt(M) - @test_throws SingularException ldlt!(F, M) + @test_throws ZeroPivotException ldlt(M) + @test_throws ZeroPivotException ldlt!(F, M) @test !issuccess(ldlt(M; check = false)) @test !issuccess(ldlt!(F, M; check = false)) end