RFC: Delay throwing when CHOLMOD factorizations fail

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
@@ -464,11 +469,12 @@ end
 isposdef(C::Cholesky) = C.info == 0
 
 function det(C::Cholesky)
+    isposdef(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
     end
-    @assertposdef dd C.info
+    return dd
 end
 
 function logdet(C::Cholesky)

base/linalg/factorization.jl

@@ -16,6 +16,25 @@ 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
+```
+"""
+issuccess(F::Factorization)
+
 function logdet(F::Factorization)
     d, s = logabsdet(F)
     return d + log(s)

base/linalg/linalg.jl

@@ -91,6 +91,7 @@ export
     ishermitian,
     isposdef,
     isposdef!,
+    issuccess,
     issymmetric,
     istril,
     istriu,

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
 
@@ -789,6 +789,14 @@ function solve(sys::Integer, F::Factor{Tv}, B::Dense{Tv}) where Tv<:VTypes
         throw(DimensionMismatch("LHS and RHS should have the same number of rows. " *
             "LHS has $(size(F,1)) rows, but RHS has $(size(B,1)) rows."))
     end
+    if !issuccess(F)
+        s = unsafe_load(F.p)
+        if s.is_ll == 1
+            throw(LinAlg.PosDefException(s.minor))
+        else
+            throw(ArgumentError("factorized matrix has one or more zero pivots. Try using lufact instead."))
+        end
+    end
     d = Dense(ccall((@cholmod_name("solve", SuiteSparse_long),:libcholmod), Ptr{C_Dense{Tv}},
             (Cint, Ptr{C_Factor{Tv}}, Ptr{C_Dense{Tv}}, Ptr{UInt8}),
                 sys, get(F.p), get(B.p), common()))
@@ -1195,6 +1203,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 +1365,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 +1404,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 +1448,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 +1503,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 +1704,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{eltype(A), SuiteSparse_long}, A)), B)
+            end
         end
     end
 end
@@ -1750,17 +1760,27 @@ 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
+            # try conversion to LLt
+            change_factor!(eltype(F), true, s.is_super, true, s.is_monotonic, F)
+            b = issuccess(F)
+            # convert back
+            change_factor!(eltype(F), false, s.is_super, true, s.is_monotonic, F)
+            return b
+        end
+    else
         return false
     end
-    true
 end
 
 function ishermitian(A::Sparse{Float64})

base/sparse/linalg.jl

@@ -890,33 +890,30 @@ function factorize(A::SparseMatrixCSC)
             return UpperTriangular(A)
         end
         if ishermitian(A)
-            try
-                return cholfact(Hermitian(A))
-            catch e
-                isa(e, PosDefException) || rethrow(e)
-                return ldltfact(Hermitian(A))
-            end
+            return factorize(Hermitian(A))
         end
         return lufact(A)
     else
         return qrfact(A)
     end
 end
 
-function factorize(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) where Ti
-    try
-        return cholfact(A)
-    catch e
-        isa(e, PosDefException) || rethrow(e)
-        return ldltfact(A)
-    end
-end
-function factorize(A::Hermitian{Complex{Float64}, SparseMatrixCSC{Complex{Float64},Ti}}) where Ti
-    try
-        return cholfact(A)
-    catch e
-        isa(e, PosDefException) || rethrow(e)
-        return ldltfact(A)
+# function factorize(A::Symmetric{Float64,SparseMatrixCSC{Float64,Ti}}) where Ti
+#     F = cholfact(A)
+#     if LinAlg.issuccess(F)
+#         return F
+#     else
+#         ldltfact!(F, A)
+#         return F
+#     end
+# end
+function factorize(A::LinAlg.RealHermSymComplexHerm{Float64,<:SparseMatrixCSC})
+    F = cholfact(A)
+    if LinAlg.issuccess(F)
+        return F
+    else
+        ldltfact!(F, A)
+        return F
     end
 end
 

doc/src/stdlib/linalg.md

@@ -95,6 +95,7 @@ Base.LinAlg.logm
 Base.LinAlg.sqrtm
 Base.LinAlg.lyap
 Base.LinAlg.sylvester
+Base.LinAlg.issuccess
 Base.LinAlg.issymmetric
 Base.LinAlg.isposdef
 Base.LinAlg.isposdef!