Merge pull request #770 from tpapp/tp/det-cleanup

Define logabsdet(::SMatrix), det & friends for LU

Author: andyferris

src/StaticArrays.jl

@@ -18,7 +18,7 @@ using LinearAlgebra
 import LinearAlgebra: transpose, adjoint, dot, eigvals, eigen, lyap, tr,
                       kron, diag, norm, dot, diagm, lu, svd, svdvals,
                       factorize, ishermitian, issymmetric, isposdef, issuccess, normalize,
-                      normalize!, Eigen, det, logdet, cross, diff, qr, \
+                      normalize!, Eigen, det, logdet, logabsdet, cross, diff, qr, \
 using LinearAlgebra: checksquare
 
 export SOneTo
@@ -128,6 +128,7 @@ include("arraymath.jl")
 include("linalg.jl")
 include("matrix_multiply_add.jl")
 include("matrix_multiply.jl")
+include("lu.jl")
 include("det.jl")
 include("inv.jl")
 include("solve.jl")
@@ -138,7 +139,6 @@ include("lyap.jl")
 include("triangular.jl")
 include("cholesky.jl")
 include("svd.jl")
-include("lu.jl")
 include("qr.jl")
 include("deque.jl")
 include("flatten.jl")

src/det.jl

@@ -42,13 +42,31 @@ end
     -2*rem(s, 2) + 1
 end
 
+det(F::LU) = det(F.U) * _parity(F.p)
+
+function logabsdet(A::Union{LowerTriangular{<:Any,<:StaticMatrix},
+                            UpperTriangular{<:Any,<:StaticMatrix}})
+    checksquare(A)
+    mapreduce(x -> (log(abs(x)), sign(x)), ((l1, s1), (l2, s2)) -> (l1 + l2, s1 * s2),
+              diag(A))
+end
+
+function logdet(F::LU)
+    d, s = logabsdet(F.U)
+    d + log(s * _parity(F.p))
+end
+
+function logabsdet(F::LU)
+    d, s = logabsdet(F.U)
+    d, s * _parity(F.p)
+end
+
 @generated function _det(::Size{S}, A::StaticMatrix) where S
     checksquare(A)
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            _, U, p = _lu(A, Val(true), false)
-            det(UpperTriangular(U))*_parity(p)
+            det(lu(A, Val(true); check = false))
         end
     else
         :(@_inline_meta; det(Matrix(A)))
@@ -62,11 +80,22 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            LUp = lu(A; check = false)
-            d, s = logabsdet(LUp.U)
-            d + log(s*_parity(LUp.p))
+            logdet(lu(A, Val(true); check = false))
         end
     else
         :(@_inline_meta; logdet(drop_sdims(A)))
     end
 end
+
+@inline logabsdet(A::StaticMatrix) = _logabsdet(Size(A), A)
+@generated function _logabsdet(::Size{S}, A::StaticMatrix) where S
+    checksquare(A)
+    if prod(S) ≤ 14*14
+        quote
+            @_inline_meta
+            logabsdet(lu(A, Val(true); check = false))
+        end
+    else
+        :(@_inline_meta; logabsdet(drop_sdims(A)))
+    end
+end

src/lu.jl

@@ -84,7 +84,7 @@ issuccess(F::LU) = _first_zero_on_diagonal(F.U) == 0
             L = similar_type(A, T2, Size($M, $(min(M,N))))(f.L)
             U = similar_type(A, T2, Size($(min(M,N)), $N))(f.U)
             p = similar_type(A, Int, Size($M))(f.p)
-            (L,U,p)
+            L, U, p
         end
     end
 end

test/det.jl

@@ -29,12 +29,17 @@ using StaticArrays, Test, LinearAlgebra
     for sz in (5, 8, 15), typ in (Float64, Complex{Float64})
         A = rand(typ, sz, sz)
         SA = SMatrix{sz,sz,typ}(A)
-        @test det(A) ≈ det(SA)
+        @test det(A) ≈ det(SA) == det(lu(SA))
         if typ == Float64 && det(A) < 0
             A[:,1], A[:,2] = A[:,2], A[:,1]
             SA = SMatrix{sz,sz,typ}(A)
         end
-        @test logdet(A) ≈ logdet(SA)
+        @test logdet(A) ≈ logdet(SA) == logdet(lu(SA))
+        dA, sA = logabsdet(A)
+        dSA, sSA = logabsdet(SA)
+        dLU, sLU = logabsdet(lu(SA))
+        @test dA ≈ dSA == dLU
+        @test sA ≈ sSA == sLU
     end
 
     @test_throws DimensionMismatch det(@SMatrix [0; 1])