Merge pull request #13166 from JuliaLang/anj/symdet

Add det(Symmetric/Hermitian) based on bkfact.

Author: andreasnoack

base/linalg/bunchkaufman.jl

@@ -50,6 +50,31 @@ function A_ldiv_B!{T<:BlasComplex}(B::BunchKaufman{T}, R::StridedVecOrMat{T})
     (issym(B) ? LAPACK.sytrs! : LAPACK.hetrs!)(B.uplo, B.LD, B.ipiv, R)
 end
 
+function det(F::BunchKaufman)
+    M = F.LD
+    p = F.ipiv
+    n = size(F.LD, 1)
+
+    i = 1
+    d = one(eltype(F))
+    while i <= n
+        if p[i] > 0
+            # 1x1 pivot case
+            d *= M[i,i]
+            i += 1
+        else
+            # 2x2 pivot case. Make sure not to square before the subtraction by scaling with the off-diagonal element. This is safe because the off diagonal is always large for 2x2 pivots.
+            if F.uplo == 'U'
+                d *= M[i, i + 1]*(M[i,i]/M[i, i + 1]*M[i + 1, i + 1] - (issym(F) ? M[i, i + 1] : conj(M[i, i + 1])))
+            else
+                d *= M[i + 1,i]*(M[i, i]/M[i + 1, i]*M[i + 1, i + 1] - (issym(F) ? M[i + 1, i] : conj(M[i + 1, i])))
+            end
+            i += 2
+        end
+    end
+    return d
+end
+
 ## reconstruct the original matrix
 ## TODO: understand the procedure described at
 ## http://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf

base/linalg/symmetric.jl

@@ -117,8 +117,16 @@ A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{
 *(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
 *(A::StridedMatrix, B::HermOrSym) = A*full(B)
 
-factorize(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A))
+bkfact(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A))
+factorize(A::HermOrSym) = bkfact(A)
+
+# Is just RealHermSymComplexHerm, but type alias seems to be broken
+det{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}) = real(det(bkfact(A)))
+det{T<:Real}(A::Symmetric{T}) = det(bkfact(A))
+det(A::Symmetric) = det(bkfact(A))
+
 \{T,S<:StridedMatrix}(A::HermOrSym{T,S}, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B)
+
 inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfact(A.data, symbol(A.uplo))), A.uplo)
 inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, symbol(A.uplo), true)), A.uplo)
 

test/linalg/symmetric.jl

@@ -101,6 +101,16 @@ let n=10
         # cond
         @test cond(Hermitian(asym)) ≈ cond(asym)
 
+        # det
+        @test det(asym) ≈ det(Hermitian(asym, :U))
+        @test det(asym) ≈ det(Hermitian(asym, :L))
+        if eltya <: Real
+            @test det(asym) ≈ det(Symmetric(asym, :U))
+            @test det(asym) ≈ det(Symmetric(asym, :L))
+        end
+        @test det(a + a.') ≈ det(Symmetric(a + a.', :U))
+        @test det(a + a.') ≈ det(Symmetric(a + a.', :L))
+
         # rank
         let A = a[:,1:5]*a[:,1:5]'
             @test rank(A) == rank(Hermitian(A))