More tests for LU decomposition

Author: martinholters

test/lu.jl

@@ -1,21 +1,41 @@
 using StaticArrays, Base.Test
 
-@testset "LU decomposition" begin
-    # Square case
-    m22 = @SMatrix [1 2; 3 4]
-    @test @inferred(lu(m22)) isa Tuple{LowerTriangular{Float64,SMatrix{2,2,Float64,4}}, UpperTriangular{Float64,SMatrix{2,2,Float64,4}}, SVector{2,Int}}
-    @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]
+@testset "LU decomposition (pivot=$pivot)" for pivot in (true, false)
+    @testset "$m×$n" for m in 0:4, n in 0:4
+        a = SMatrix{m,n,Int}(1:(m*n))
+        l, u, p = @inferred(lu(a, Val{pivot}))
 
-    # Rectangular case
-    m23 = @SMatrix Float64[3 9 4; 6 6 2]
-    @test @inferred(lu(m23)) isa Tuple{SMatrix{2,2,Float64,4}, SMatrix{2,3,Float64,6}, SVector{2,Int}}
-    @test lu(m23)[1] ≊ lu(Matrix(m23))[1]
-    @test lu(m23)[2] ≊ lu(Matrix(m23))[2]
-    @test lu(m23)[3] ≊ lu(Matrix(m23))[3]
+        # expected types
+        @test p isa SVector{m,Int}
+        if m==n
+            @test l isa LowerTriangular{<:Any,<:SMatrix{m,n}}
+            @test u isa UpperTriangular{<:Any,<:SMatrix{m,n}}
+        else
+            @test l isa SMatrix{m,min(m,n)}
+            @test u isa SMatrix{min(m,n),n}
+        end
 
-    @test lu(m23')[1] ≊ lu(Matrix(m23'))[1]
-    @test lu(m23')[2] ≊ lu(Matrix(m23'))[2]
-    @test lu(m23')[3] ≊ lu(Matrix(m23'))[3]
+        if pivot
+            # p is a permutation
+            @test sort(p) == collect(1:m)
+        else
+            @test p == collect(1:m)
+        end
+
+        # l is unit lower triangular
+        for i=1:m, j=(i+1):size(l,2)
+            @test iszero(l[i,j])
+        end
+        for i=1:size(l,2)
+            @test l[i,i] == 1
+        end
+
+        # u is upper triangular
+        for i=1:size(u,1), j=1:i-1
+            @test iszero(u[i,j])
+        end
+
+        # decomposition is correct
+        @test l*u ≈ a[p,:]
+    end
 end