Add det(Symmetric/Hermitian) based on bkfact.

[andreasnoack]

Fixes JuliaLang/LinearAlgebra.jl#252.

[fredrikekre]

Should this not try cholfact first @andreasnoack ? Since this is the generic factorization method. Currently:

julia> A = rand(Complex128, 4, 4); A = A'A; H = Hermitian(A);

julia> typeof(factorize(A))
Base.LinAlg.Cholesky{Complex{Float64},Array{Complex{Float64},2}}

julia> typeof(factorize(H))
Base.LinAlg.BunchKaufman{Complex{Float64},Array{Complex{Float64},2}}

and then change this and this to bkfact(Hermitian(A)) and bkfact(Symmetric(A)) ?

[tkelman]

cholfact can fail, bunch kaufman won't, so why would starting with cholfact be better?

[fredrikekre]

This type of check is already happening in factorize? Since factorize is a general method which essentially should return the "best" factorization I think the same thing should apply if you have a Hermitian input matrix?

[fredrikekre]

(And since #21976 cholfact does not throw)

[tkelman]

it doesn't throw, but it gives you an answer that you can't use for many indefinite inputs

[fredrikekre]

That is already what is happening (unless the matrix is annotated with Hermitian/Symmetric). Are you suggesting we should not even try the Cholesky factorization in factorize and do this:

diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index ab2dac2..c39940f 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -768,12 +768,7 @@ function factorize(A::StridedMatrix{T}) where T
             return UpperTriangular(A)
         end
         if herm
-            cf = cholfact(A)
-            if cf.info == 0
-                return cf
-            else
-                return factorize(Hermitian(A))
-            end
+            return factorize(Hermitian(A))
         end
         if sym
             return factorize(Symmetric(A))

[tkelman]

it's a bit of a judgement call on how frequent you expect indefinite input to be and worth benchmarking the cost of trying cholesky before bunch kaufman vs going to bunch kaufman first. The latter has a faster worst case but can be a bit slower on inputs that cholesky works for, which may be more common in some fields

[fredrikekre]

Okay, what I am really after is that my interpretation of the documentation is that this should return a Cholesky factorization since H is positive definite:

julia> A = rand(3, 3); A = A'A; H = Hermitian(A);

julia> isposdef(H)
true

julia> typeof(factorize(H))
Base.LinAlg.BunchKaufman{Float64,Array{Float64,2}}

[fredrikekre]

My original suggestion here would also make it consistent with sparse matrices https://github.com/JuliaLang/julia/blob/master/base/sparse/linalg.jl#L927-L942

[tkelman]

yeah, part of the problem is cholmod can't do sparse bkfact, there are other libraries but it would take some thought to figure out how to make them sanely coexist