LU fixes

Author: fredrikekre

src/det.jl

@@ -47,8 +47,8 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            L, U, p = lu(A)
-            det(U)*_parity(p)
+            LUp = lu(A)
+            det(LUp.U)*_parity(LUp.p)
         end
     else
         :(@_inline_meta; det(Matrix(A)))
@@ -62,9 +62,9 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            L, U, p = lu(A)
-            d, s = logabsdet(U)
-            d + log(s*_parity(p))
+            LUp = lu(A)
+            d, s = logabsdet(LUp.U)
+            d + log(s*_parity(LUp.p))
         end
     else
         :(@_inline_meta; logdet(drop_sdims(A)))

src/inv.jl

@@ -60,8 +60,8 @@ end
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
-            L, U, p = lu(A)
-            U \ (L \ eye(A)[p,:])
+            LUp = lu(A)
+            LUp.U \ (LUp.L \ eye(A)[LUp.p,:])
         end
     else
         :(@_inline_meta; similar_type(A)(inv(Matrix(A))))

src/lu.jl

@@ -1,14 +1,27 @@
+# define our own LU type, since LinearAlgebra.LU requires p::Vector
+struct LU{L,U,p}
+    L::L
+    U::U
+    p::p
+end
+
+# iteration for destructuring into components
+Base.iterate(S::LU) = (S.L, Val(:U))
+Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p))
+Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done))
+Base.iterate(S::LU, ::Val{:done}) = nothing
+
 # LU decomposition
 function lu(A::StaticMatrix, pivot::Union{Val{false},Val{true}}=Val(true))
-    L,U,p = _lu(A, pivot)
-    (L,U,p)
+    L, U, p = _lu(A, pivot)
+    LU(L, U, p)
 end
 
 # For the square version, return explicit lower and upper triangular matrices.
 # We would do this for the rectangular case too, but Base doesn't support that.
 function lu(A::StaticMatrix{N,N}, pivot::Union{Val{false},Val{true}}=Val(true)) where {N}
-    L,U,p = _lu(A, pivot)
-    (LowerTriangular(L), UpperTriangular(U), p)
+    L, U, p = _lu(A, pivot)
+    LU(LowerTriangular(L), UpperTriangular(U), p)
 end
 
 @generated function _lu(A::StaticMatrix{M,N,T}, pivot) where {M,N,T}

src/solve.jl

@@ -34,8 +34,8 @@ end
     if prod(Sa) ≤ 14*14
         quote
             @_inline_meta
-            L, U, p = lu(a)
-            U \ (L \ $(length(Sb) > 1 ? :(b[p,:]) : :(b[p])))
+            LUp = lu(a)
+            LUp.U \ (LUp.L \ $(length(Sb) > 1 ? :(b[LUp.p,:]) : :(b[LUp.p])))
         end
     else
         quote

test/lu.jl

@@ -1,6 +1,5 @@
 using StaticArrays, Test, LinearAlgebra
 
-@testset "LU decomposition" begin
 @testset "LU decomposition ($m×$n, pivot=$pivot)" for pivot in (true, false), m in [0:4..., 15], n in [0:4..., 15]
     a = SMatrix{m,n,Int}(1:(m*n))
     l, u, p = @inferred(lu(a, Val{pivot}()))
@@ -38,4 +37,3 @@ using StaticArrays, Test, LinearAlgebra
     # decomposition is correct
     @test l*u ≈ a[p,:]
 end
-end