Delay throwing when CHOLMOD factorizations fail

Introduce issuccess function to test if a Factorization succeeded.

Fixes #22335

base/linalg/bunchkaufman.jl

@@ -88,8 +88,15 @@ size(B::BunchKaufman, d::Integer) = size(B.LD, d)
 issymmetric(B::BunchKaufman) = B.symmetric
 ishermitian(B::BunchKaufman) = !B.symmetric
 
+issuccess(B::BunchKaufman) = B.info == 0
+
+function Base.show(io::IO, F::BunchKaufman)
+    println(io, summary(F))
+    print(io, "successful: $(issuccess(F))")
+end
+
 function inv(B::BunchKaufman{<:BlasReal})
-    if B.info > 0
+    if !issuccess(B)
         throw(SingularException(B.info))
     end
 
@@ -101,7 +108,7 @@ function inv(B::BunchKaufman{<:BlasReal})
 end
 
 function inv(B::BunchKaufman{<:BlasComplex})
-    if B.info > 0
+    if !issuccess(B)
         throw(SingularException(B.info))
     end
 
@@ -121,7 +128,7 @@ function inv(B::BunchKaufman{<:BlasComplex})
 end
 
 function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
-    if B.info > 0
+    if !issuccess(B)
         throw(SingularException(B.info))
     end
 
@@ -132,7 +139,7 @@ function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
     end
 end
 function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasComplex
-    if B.info > 0
+    if !issuccess(B)
         throw(SingularException(B.info))
     end
 
@@ -161,7 +168,7 @@ function logabsdet(F::BunchKaufman)
     p = F.ipiv
     n = size(F.LD, 1)
 
-    if F.info > 0
+    if !issuccess(F)
         return eltype(F)(-Inf), zero(eltype(F))
     end
     s = one(real(eltype(F)))

base/linalg/cholesky.jl

@@ -410,8 +410,13 @@ function getindex(C::CholeskyPivoted{T}, d::Symbol) where T<:BlasFloat
     throw(KeyError(d))
 end
 
-show(io::IO, C::Cholesky{<:Any,<:AbstractMatrix}) =
-    (println(io, "$(typeof(C)) with factor:");show(io,C[:UL]))
+issuccess(C::Cholesky) = C.info == 0
+
+function show(io::IO, C::Cholesky{<:Any,<:AbstractMatrix})
+    println(io, "$(typeof(C)) with factor:")
+    show(io, C[:UL])
+    print(io, "\nsuccessful: $(issuccess(C))")
+end
 
 A_ldiv_B!(C::Cholesky{T,<:AbstractMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
     @assertposdef LAPACK.potrs!(C.uplo, C.factors, B) C.info
@@ -462,7 +467,7 @@ function A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix)
 end
 
 function det(C::Cholesky)
-    C.info == 0 || throw(PosDefException(C.info))
+    issuccess(C) || throw(PosDefException(C.info))
     dd = one(real(eltype(C)))
     @inbounds for i in 1:size(C.factors,1)
         dd *= real(C.factors[i,i])^2
@@ -471,7 +476,7 @@ function det(C::Cholesky)
 end
 
 function logdet(C::Cholesky)
-    C.info == 0 || throw(PosDefException(C.info))
+    issuccess(C) || throw(PosDefException(C.info))
     dd = zero(real(eltype(C)))
     @inbounds for i in 1:size(C.factors,1)
         dd += log(real(C.factors[i,i]))

base/linalg/factorization.jl

@@ -16,6 +16,24 @@ macro assertnonsingular(A, info)
    :($(esc(info)) == 0 ? $(esc(A)) : throw(SingularException($(esc(info)))))
 end
 
+"""
+    issuccess(F::Factorization)
+
+Test that a factorization of a matrix succeeded.
+
+```jldoctest
+julia> cholfact([1 0; 0 1])
+Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
+[1.0 0.0; 0.0 1.0]
+successful: true
+
+julia> cholfact([1 0; 0 -1])
+Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
+[1.0 0.0; 0.0 -1.0]
+successful: false
+```
+"""
+
 function logdet(F::Factorization)
     d, s = logabsdet(F)
     return d + log(s)

base/linalg/lu.jl

@@ -147,7 +147,7 @@ function lufact(A::AbstractMatrix{T}) where T
     AA = similar(A, S, size(A))
     copy!(AA, A)
     F = lufact!(AA, Val{false})
-    if F.info == 0
+    if issuccess(F)
         return F
     else
         AA = similar(A, S, size(A))
@@ -227,11 +227,14 @@ function getindex(F::LU{T,<:StridedMatrix}, d::Symbol) where T
     end
 end
 
-function show(io::IO, C::LU)
-    println(io, "$(typeof(C)) with factors L and U:")
-    show(io, C[:L])
+issuccess(F::LU) = F.info == 0
+
+function show(io::IO, F::LU)
+    println(io, "$(typeof(F)) with factors L and U:")
+    show(io, F[:L])
     println(io)
-    show(io, C[:U])
+    show(io, F[:U])
+    print(io, "\nsuccessful: $(issuccess(F))")
 end
 
 A_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
@@ -271,31 +274,31 @@ Ac_ldiv_Bc(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasComple
     @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info
 Ac_ldiv_Bc(A::LU, B::StridedVecOrMat) = Ac_ldiv_B(A, ctranspose(B))
 
-function det(A::LU{T}) where T
-    n = checksquare(A)
-    A.info > 0 && return zero(T)
+function det(F::LU{T}) where T
+    n = checksquare(F)
+    issuccess(F) || return zero(T)
     P = one(T)
     c = 0
     @inbounds for i = 1:n
-        P *= A.factors[i,i]
-        if A.ipiv[i] != i
-                c += 1
+        P *= F.factors[i,i]
+        if F.ipiv[i] != i
+            c += 1
         end
     end
     s = (isodd(c) ? -one(T) : one(T))
     return P * s
 end
 
-function logabsdet(A::LU{T}) where T  # return log(abs(det)) and sign(det)
-    n = checksquare(A)
-    A.info > 0 && return log(zero(real(T))), log(one(T))
+function logabsdet(F::LU{T}) where T  # return log(abs(det)) and sign(det)
+    n = checksquare(F)
+    issuccess(F) || return log(zero(real(T))), log(one(T))
     c = 0
     P = one(T)
     abs_det = zero(real(T))
     @inbounds for i = 1:n
-        dg_ii = A.factors[i,i]
+        dg_ii = F.factors[i,i]
         P *= sign(dg_ii)
-        if A.ipiv[i] != i
+        if F.ipiv[i] != i
             c += 1
         end
         abs_det += log(abs(dg_ii))

base/sparse/cholmod.jl

@@ -7,7 +7,7 @@ import Base: (*), convert, copy, eltype, get, getindex, show, size,
 
 import Base.LinAlg: (\), A_mul_Bc, A_mul_Bt, Ac_ldiv_B, Ac_mul_B, At_ldiv_B, At_mul_B,
                  cholfact, cholfact!, det, diag, ishermitian, isposdef,
-                 issymmetric, ldltfact, ldltfact!, logdet
+                 issuccess, issymmetric, ldltfact, ldltfact!, logdet
 
 importall ..SparseArrays
 
@@ -1195,6 +1195,7 @@ function showfactor(io::IO, F::Factor)
     @printf(io, "method: %10s\n", s.is_super!=0 ? "supernodal" : "simplicial")
     @printf(io, "maxnnz: %10d\n", Int(s.nzmax))
     @printf(io, "nnz: %13d\n", nnz(F))
+    @printf(io, "success: %9s\n", "$(s.minor == size(F, 1))")
 end
 
 # getindex not defined for these, so don't use the normal array printer
@@ -1356,8 +1357,6 @@ function cholfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv
     # Compute the numerical factorization
     factorize_p!(A, shift, F, cm)
 
-    s = unsafe_load(get(F.p))
-    s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor))
     return F
 end
 
@@ -1397,8 +1396,6 @@ function cholfact(A::Sparse; shift::Real=0.0,
     # Compute the numerical factorization
     cholfact!(F, A; shift = shift)
 
-    s = unsafe_load(get(F.p))
-    s.minor < size(A, 1) && throw(Base.LinAlg.PosDefException(s.minor))
     return F
 end
 
@@ -1443,13 +1440,17 @@ cholfact(A::Union{SparseMatrixCSC{T}, SparseMatrixCSC{Complex{T}},
 
 
 function ldltfact!(F::Factor{Tv}, A::Sparse{Tv}; shift::Real=0.0) where Tv
-    cm = common()
+    cm = defaults(common())
+    set_print_level(cm, 0)
+
+    # Makes it an LDLt
+    unsafe_store!(common_final_ll, 0)
+    # Really make sure it's an LDLt by avoiding supernodal factorization
+    unsafe_store!(common_supernodal, 0)
 
     # Compute the numerical factorization
     factorize_p!(A, shift, F, cm)
 
-    s = unsafe_load(get(F.p))
-    s.minor < size(A, 1) && throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots"))
     return F
 end
 
@@ -1494,10 +1495,6 @@ function ldltfact(A::Sparse; shift::Real=0.0,
     # Compute the numerical factorization
     ldltfact!(F, A; shift = shift)
 
-    s = unsafe_load(get(F.p))
-    if s.minor < size(A, 1)
-        throw(Base.LinAlg.ArgumentError("matrix has one or more zero pivots"))
-    end
     return F
 end
 
@@ -1699,11 +1696,16 @@ for f in (:\, :Ac_ldiv_B)
     @eval function ($f)(A::Union{Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
                           Hermitian{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}},
                           Hermitian{Complex{Float64},SparseMatrixCSC{Complex{Float64},SuiteSparse_long}}}, B::StridedVecOrMat)
-        try
-            return ($f)(cholfact(A), B)
-        catch e
-            isa(e, LinAlg.PosDefException) || rethrow(e)
-            return ($f)(ldltfact(A) , B)
+        F = cholfact(A)
+        if issuccess(F)
+            return ($f)(F, B)
+        else
+            ldltfact!(F, A)
+            if issuccess(F)
+                return ($f)(F, B)
+            else
+                return ($f)(lufact(convert(SparseMatrixCSC{Float64, SuiteSparse_long}, A)), B)
+            end
         end
     end
 end
@@ -1750,17 +1752,22 @@ end
 
 det(L::Factor) = exp(logdet(L))
 
-function isposdef(A::SparseMatrixCSC{<:VTypes,SuiteSparse_long})
-    if !ishermitian(A)
-        return false
-    end
-    try
-        f = cholfact(A)
-    catch e
-        isa(e, LinAlg.PosDefException) || rethrow(e)
+function issuccess(F::Factor)
+    s = unsafe_load(F.p)
+    return s.minor == size(F, 1)
+end
+
+function isposdef(F::Factor)
+    if issuccess(F)
+        s = unsafe_load(F.p)
+        if s.is_ll == 1
+            return true
+        else
+            return :maybe
+        end
+    else
         return false
     end
-    true
 end
 
 function ishermitian(A::Sparse{Float64})

test/linalg/bunchkaufman.jl

@@ -50,6 +50,7 @@ bimg  = randn(n,2)/2
 
                 @testset "(Automatic) Bunch-Kaufman factor of indefinite matrix" begin
                     bc1 = factorize(asym)
+                    @test LinAlg.issuccess(bc1)
                     @test logabsdet(bc1)[1] ≈ log(abs(det(bc1)))
                     if eltya <: Real
                         @test logabsdet(bc1)[2] == sign(det(bc1))
@@ -72,6 +73,7 @@ bimg  = randn(n,2)/2
                 @testset "Bunch-Kaufman factors of a pos-def matrix" begin
                     @testset for rook in (false, true)
                         bc2 = bkfact(apd, :U, issymmetric(apd), rook)
+                        @test LinAlg.issuccess(bc2)
                         @test logdet(bc2) ≈ log(det(bc2))
                         @test logabsdet(bc2)[1] ≈ log(abs(det(bc2)))
                         @test logabsdet(bc2)[2] == sign(det(bc2))
@@ -103,6 +105,7 @@ end
 
                 @testset for rook in (false, true)
                     F = bkfact(As, :U, issymmetric(As), rook)
+                    @test !LinAlg.issuccess(F)
                     @test det(F) == 0
                     @test_throws LinAlg.SingularException inv(F)
                     @test_throws LinAlg.SingularException F \ ones(size(As, 1))

test/linalg/cholesky.jl

@@ -51,6 +51,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
             @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j]))
         end
         @test apd*inv(capd) ≈ eye(n)
+        @test LinAlg.issuccess(capd)
         @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit
         @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow
 
@@ -73,7 +74,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
                 @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n
             end
             ulstring = sprint(show,capds[:UL])
-            @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring"
+            @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring\nsuccessful: true"
         else
             capdh = cholfact(apdh)
             @test inv(capdh)*apdh ≈ eye(n)
@@ -85,7 +86,7 @@ using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted, PosDefException
             @test inv(capdh)*apdh ≈ eye(n)
             @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
             ulstring = sprint(show,capdh[:UL])
-            @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring"
+            @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring\nsuccessful: true"
         end
 
         # test chol of 2x2 Strang matrix
@@ -270,6 +271,7 @@ end
     for T in (Float32, Float64, Complex64, Complex128)
         A = T[1 2; 2 1]; B = T[1, 1]
         C = cholfact(A)
+        @test !LinAlg.issuccess(C)
         @test_throws PosDefException C\B
         @test_throws PosDefException det(C)
         @test_throws PosDefException logdet(C)

test/linalg/lu.jl

@@ -67,7 +67,7 @@ debug && println("(Automatic) Square LU decomposition. eltya: $eltya, eltyb: $el
 
         lstring = sprint(show,l)
         ustring = sprint(show,u)
-        @test sprint(show,lua) == "$(typeof(lua)) with factors L and U:\n$lstring\n$ustring"
+        @test sprint(show,lua) == "$(typeof(lua)) with factors L and U:\n$lstring\n$ustring\nsuccessful: true"
         let Bs = b, Cs = c
             for atype in ("Array", "SubArray")
                 if atype == "Array"
@@ -97,6 +97,7 @@ debug && println("(Automatic) Square LU decomposition. eltya: $eltya, eltyb: $el
 debug && println("Tridiagonal LU")
         κd    = cond(Array(d),1)
         lud   = lufact(d)
+        @test LinAlg.issuccess(lud)
         @test lufact(lud) == lud
         @test_throws KeyError lud[:Z]
         @test lud[:L]*lud[:U] ≈ lud[:P]*Array(d)
@@ -139,7 +140,7 @@ debug && println("Tridiagonal LU")
             du[1] = zero(eltya)
             dl[1] = zero(eltya)
             zT = Tridiagonal(dl,dd,du)
-            @test lufact(zT).info == 1
+            @test !LinAlg.issuccess(lufact(zT))
         end
 
 debug && println("Thin LU")
@@ -211,3 +212,6 @@ end
 
 # Issue 21453.
 @test_throws ArgumentError LinAlg._cond1Inf(lufact(randn(5,5)), 2, 2.0)
+
+# Singular LU
+@test !LinAlg.issuccess(lufact(zeros(3,3)))