fix exception type for sparse LDLᵀ

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
 
 

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

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

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)

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
 

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)

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.PosDefException(1)))
+    check && (issuccess(F) || throw(LinearAlgebra.ZeroPivotException(1)))
     return F
 end
 

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, 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 PosDefException ldlt(C - C[1,1]*I)
-    @test_throws PosDefException 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 PosDefException ldlt(M)
-        @test_throws PosDefException 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