LU decomposition via Base

For now, a simple implementation which calls Base.lu, nothing fancy -
just preserve the static sizes correctly in the output.

Author: c42f

src/StaticArrays.jl

@@ -9,7 +9,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
              fill!, det, inv, eig, eigvals, expm, sqrtm, trace, vecnorm, norm, dot, diagm, diag,
-             svd, svdvals!, svdfact,
+             lu, svd, svdvals!, svdfact,
              sum, diff, prod, count, any, all, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,
              randexp!, normalize, normalize!, read, read!, write
@@ -94,9 +94,10 @@ include("eigen.jl")
 include("expm.jl")
 include("sqrtm.jl")
 include("cholesky.jl")
+include("svd.jl")
+include("lu.jl")
 include("deque.jl")
 include("io.jl")
-include("svd.jl")
 
 include("FixedSizeArrays.jl")
 include("ImmutableArrays.jl")

src/linalg.jl

@@ -324,19 +324,3 @@ end
 @inline Base.LinAlg.Symmetric(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
 @inline Base.LinAlg.Hermitian(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))
 
-
-#-------------------------------------------------------------------------------
-# Wrappers for various matrix decompositions which correctly propagate the
-# StaticMatrix type in the returned factorization.
-
-
-#--------------------------------------------------
-# lu
-# TODO
-
-
-#--------------------------------------------------
-# qr
-# TODO
-
-#--------------------------------------------------

src/lu.jl

@@ -0,0 +1,28 @@
+# LU decomposition
+function lu(A::StaticMatrix, pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{true})
+    L,U,p = _lu(A, pivot)
+    (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{Type{Val{false}},Type{Val{true}}}=Val{true}) where {N}
+    L,U,p = _lu(A, pivot)
+    (LowerTriangular(L), UpperTriangular(U), p)
+end
+
+
+@inline function _lu(A::StaticMatrix, pivot)
+    # For now, just call through to Base.
+    # TODO: statically sized LU without allocations!
+    f = lufact(Matrix(A), pivot)
+    T = eltype(A)
+    T2 = typeof((one(T)*zero(T) + zero(T))/one(T))
+    L = similar_type(A, T2, Size(Size(A)[1], diagsize(A)))(f[:L])
+    U = similar_type(A, T2, Size(diagsize(A), Size(A)[2]))(f[:U])
+    p = similar_type(A, Size(Size(A)[1]))(f[:p])
+    (L,U,p)
+end
+
+# Base.lufact() interface is fairly inherently type unstable.  Punt on
+# implementing that, for now...

test/lu.jl

@@ -0,0 +1,21 @@
+using StaticArrays, Base.Test
+
+@testset "LU decomposition" begin
+    # Square case
+    m22 = @SMatrix [1 2; 3 4]
+    @inferred(lu(m22))
+    @test lu(m22)[1]::LowerTriangular{<:Any,<:StaticMatrix} ≊ lu(Matrix(m22))[1]
+    @test lu(m22)[2]::UpperTriangular{<:Any,<:StaticMatrix} ≊ lu(Matrix(m22))[2]
+    @test lu(m22)[3]::StaticVector ≊ lu(Matrix(m22))[3]
+
+    # Rectangular case
+    m23 = @SMatrix Float64[3 9 4; 6 6 2]
+    @inferred lu(m23)
+    @test lu(m23)[1] ≊ lu(Matrix(m23))[1]
+    @test lu(m23)[2] ≊ lu(Matrix(m23))[2]
+    @test lu(m23)[3] ≊ lu(Matrix(m23))[3]
+
+    @test lu(m23')[1] ≊ lu(Matrix(m23'))[1]
+    @test lu(m23')[2] ≊ lu(Matrix(m23'))[2]
+    @test lu(m23')[3] ≊ lu(Matrix(m23'))[3]
+end