Tests for inv() and hacky workarounds for Unsigned

Partially address #29

Author: Andy Ferris

src/inv.jl

@@ -1,47 +1,49 @@
+@inline inv(A::StaticMatrix) = _inv(Size(A), A)
 
-@generated function inv{T}(A::StaticMatrix{T})
+@inline _inv(::Size{(1,1)}, A) = similar_type(typeof(A), typeof(inv(one(eltype(A)))))(inv(A[1]))
+
+@inline function _inv(::Size{(2,2)}, A)
+    T = eltype(A)
+    S = typeof((one(T)*zero(T) + zero(T))/one(T))
+    newtype = similar_type(A, S)
+
+    d = det(A)
+    @inbounds return newtype((A[4]/d, -(A[2]/d), -(A[3]/d), A[1]/d))
+end
+
+@inline function _inv(::Size{(3,3)}, A)
+    T = eltype(A)
     S = typeof((one(T)*zero(T) + zero(T))/one(T))
     newtype = similar_type(A, S)
 
-    if size(A) == (1,1)
-        return quote
-            $(Expr(:meta, :inline))
-            @inbounds return $newtype((inv(A[1]),))
-        end
-    elseif size(A) == (2,2)
-        return quote
-            $(Expr(:meta, :inline))
-            d = A[1]*A[4] - A[3]*A[2]
-            @inbounds return $newtype((A[4]/d, -A[2]/d, -A[3]/d, A[1]/d))
-        end
-    elseif size(A) == (3,3)
-        return quote
-            $(Expr(:meta, :inline))
-            @inbounds x0 = SVector(A[1], A[2], A[3])
-            @inbounds x1 = SVector(A[4], A[5], A[6])
-            @inbounds x2 = SVector(A[7], A[8], A[9])
-
-            y0 = cross(x1,x2)
-            d  = vecdot(x0, y0)
-            x0 = x0 / d
-            y0 = y0 / d
-            y1 = cross(x2,x0)
-            y2 = cross(x0,x1)
-
-            @inbounds return $newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3]))
-        end
+    @inbounds x0 = SVector(A[1], A[2], A[3])
+    @inbounds x1 = SVector(A[4], A[5], A[6])
+    @inbounds x2 = SVector(A[7], A[8], A[9])
+
+    y0 = cross(x1,x2)
+    d  = vecdot(x0, y0)
+    x0 = x0 / d
+    y0 = y0 / d
+    y1 = cross(x2,x0)
+    y2 = cross(x0,x1)
+
+    @inbounds return newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3]))
+end
+
+@inline function _inv(::Size, A)
+    T = eltype(A)
+    S = typeof((one(T)*zero(T) + zero(T))/one(T))
+    AA = convert(Array{S}, A) # lufact() doesn't work with StaticArrays at the moment... and the branches below must be type-stable
+    if istriu(A)
+        Ai_ut = inv(UpperTriangular(AA))
+        # TODO double check these routines leave the parent in a clean (upper triangular) state
+        return Size(A)(parent(Ai_ut)) # Return a `SizedArray`
+    elseif istril(A)
+        Ai_lt = inv(LowerTriangular(AA))
+        # TODO double check these routines leave the parent in a clean (lower triangular) state
+        return Size(A)(parent(Ai_lt)) # Return a `SizedArray`
     else
-        return quote # Implementation from Base
-            AA = convert(Array{$S}, A) # lufact() doesn't work with StaticArrays at the moment... and the branches below must be type-stable
-            if istriu(A)
-                Ai = inv(UpperTriangular(AA))
-            elseif istril(A)
-                Ai = inv(LowerTriangular(AA))
-            else
-                Ai = inv(lufact(AA))
-            end
-            # Return a `SizedArray`
-            return $(Size(A))(convert(typeof(parent(Ai)), Ai))
-        end
+        Ai_lu = inv(lufact(AA))
+        return Size(A)(Ai_lu) # Return a `SizedArray`
     end
 end

src/linalg.jl

@@ -282,10 +282,18 @@ end
 end
 
 @generated function cross(a::StaticVector, b::StaticVector)
+    S = typeof(zero(eltype(a))*zero(eltype(b)))
     if length(a) === 3 && length(b) === 3
-        return quote
-            $(Expr(:meta, :inline))
-            similar_type(a, promote_type(eltype(a), eltype(b)))((a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]))
+        if S <: Unsigned
+            return quote
+                $(Expr(:meta, :inline))
+                similar_type(a, $(typeof(Signed(zero(eltype(a))*zero(eltype(b))))))((Signed(a[2]*b[3])-Signed(a[3]*b[2]), Signed(a[3]*b[1])-Signed(a[1]*b[3]), Signed(a[1]*b[2])-Signed(a[2]*b[1])))
+            end
+        else
+            return quote
+                $(Expr(:meta, :inline))
+                similar_type(a, $S)((a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1]))
+            end
         end
     else
         error("Cross product only defined for 3-vectors")

test/inv.jl

@@ -0,0 +1,10 @@
+@testset "Matrix inverse" begin
+    @test inv(@SMatrix [2])::SMatrix ≈ @SMatrix [0.5]
+    @test inv(@SMatrix [1 2; 2 1])::SMatrix ≈ [-1/3 2/3; 2/3 -1/3]
+    @test inv(@SMatrix [1 2 0; 2 1 0; 0 0 1])::SMatrix ≈ [-1/3 2/3 0; 2/3 -1/3 0; 0 0 1]
+    m = randn(Float64, 4,4) + eye(4) # well conditioned
+    @test inv(SMatrix{4,4}(m))::StaticMatrix ≈ inv(m)
+
+    @test inv(@SMatrix [0x01 0x02; 0x02 0x01])::SMatrix ≈ [-1/3 2/3; 2/3 -1/3]
+    @test inv(@SMatrix [0x01 0x02 0x00; 0x02 0x01 0x00; 0x00 0x00 0x01])::SMatrix ≈ [-1/3 2/3 0; 2/3 -1/3 0; 0 0 1]
+end

test/runtests.jl

@@ -19,6 +19,7 @@ using Base.Test
     include("linalg.jl")
     include("matrix_multiply.jl")
     include("det.jl")
+    include("inv.jl")
     include("solve.jl")
     include("eigen.jl")
     include("deque.jl")