fix exception type for sparse LDLᵀ

[stevengj]

Fixes JuliaLang/LinearAlgebra.jl#667.

[andreasnoack]

Unfortunately, the fix here will probably be a bit more complicated. The current fix is still not accurate as sparse([0 1; 1 1]) is not singular but will still fail because of the zero pivot. I suppose the check in

julia/stdlib/SuiteSparse/src/cholmod.jl

Lines 745 to 750 in 4f17558

	s = unsafe_load(pointer(F))
	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."))
	end

was meant to handle this but apparently it doesn't

julia> cholesky(sparse([0 1; 1 1]), check=false)\ones(2)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
 [1] solve(::Int32, ::SuiteSparse.CHOLMOD.Factor{Float64}, ::SuiteSparse.CHOLMOD.Dense{Float64}) at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:747
 [2] \ at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:1678 [inlined]
 [3] \(::SuiteSparse.CHOLMOD.Factor{Float64}, ::Array{Float64,1}) at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:1683
 [4] top-level scope at REPL[256]:1

[stevengj]

cholesky(sparse([0 1; 1 1])

The matrix [0 1; 1 1] is not positive definite, so the PosDefException is correct when you call cholesky. If you use the ldlt function (the subject of this PR) it seems to work fine:

julia> ldlt(sparse([0 1; 1 1]), check=false) \ ones(2)
ERROR: ArgumentError: factorized matrix has one or more zero pivots. Try using `lu` instead.

[andreasnoack]

Let me try again since I clearly wrote the last comment too early in the morning. What I think would still be wrong with this PR is

ldlt(sparse([0 1; 1 1])

which will now throw SingularException but should instead error the same way as

julia> ldlt(sparse([0 1; 1 1]), check=false)\ones(2)
ERROR: ArgumentError: factorized matrix has one or more zero pivots. Try using `lu` instead.
Stacktrace:
 [1] solve(::Int32, ::SuiteSparse.CHOLMOD.Factor{Float64}, ::SuiteSparse.CHOLMOD.Dense{Float64}) at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:749
 [2] \ at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:1678 [inlined]
 [3] \(::SuiteSparse.CHOLMOD.Factor{Float64}, ::Array{Float64,1}) at /Users/andreasnoack/julia/usr/share/julia/stdlib/v1.4/SuiteSparse/src/cholmod.jl:1683
 [4] top-level scope at REPL[283]:1

[stevengj]

Does SuiteSparse allow you to distinguish between the actually singular case and cases where there are just zero pivots for LDLt factorization?

Maybe we need to define a new ZeroPivot exception for ldlt.

[andreasnoack]

Does SuiteSparse allow you to distinguish between the actually singular case and cases where there are just zero pivots for LDLt factorization?

I don't think it has any idea if the matrix is singular or not. It only knows that it has hit a zero pivot.

Maybe we need to define a new ZeroPivot exception for ldlt.

Yes. That might be a good idea. Could also be useful in LU without pivoting.

[stevengj]

By the way, all of the unsafe_load(pointer(A)) code in src/cholmod.jl scares me a bit. Shouldn't this be sprinkled with GC.@preserve calls?

[stevengj]

Added the ZeroPivotException and used it for both ldlt and non-pivoted LU.

[StefanKarpinski]

What's the status here? Does this look good to you now, @andreasnoack?