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

1. define a logabsdet for SMatrix, fixes #489

2. define methods for det, logdet, logabsdet for LU (cf LinearAlgebra,
which has these methods)

3. make det, logdet, logabsdet fall back to using the LU methods (to
avoid code repetition)

4. include lu.jl first, since det.jl used code from it

5. incidental consistency fixes

Author: tpapp

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,24 @@ end
     -2*rem(s, 2) + 1
 end
 
+det(F::LU) = det(F.U) * _parity(F.p)
+
+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 +73,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])