Merge pull request #312 from mschauer/bilineardot

Fix 3x3 det and inv for complex entries

Author: c42f

src/det.jl

@@ -15,7 +15,7 @@ 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])
-    return vecdot(x0, cross(x1, x2))
+    return bilinear_vecdot(x0, cross(x1, x2))
 end
 
 @inline function _det(::Size{(4,4)}, A::StaticMatrix)

src/inv.jl

@@ -21,7 +21,7 @@ end
     @inbounds x2 = SVector{3}(A[7], A[8], A[9])
 
     y0 = cross(x1,x2)
-    d  = vecdot(x0, y0)
+    d  = bilinear_vecdot(x0, y0)
     x0 = x0 / d
     y0 = y0 / d
     y1 = cross(x2,x0)

src/linalg.jl

@@ -236,6 +236,23 @@ end
     end
 end
 
+@inline bilinear_vecdot(a::StaticArray, b::StaticArray) = _bilinear_vecdot(same_size(a, b), a, b)
+@generated function _bilinear_vecdot(::Size{S}, a::StaticArray, b::StaticArray) where {S}
+    if prod(S) == 0
+        return :(zero(promote_op(*, eltype(a), eltype(b))))
+    end
+
+    expr = :(a[1] * b[1])
+    for j = 2:prod(S)
+        expr = :($expr + a[$j] * b[$j])
+    end
+
+    return quote
+        @_inline_meta
+        @inbounds return $expr
+    end
+end
+
 @inline norm(v::StaticVector) = vecnorm(v)
 @inline norm(v::StaticVector, p::Real) = vecnorm(v, p)
 

test/det.jl

@@ -3,7 +3,8 @@
     @test logdet(@SMatrix [1]) == 0.0
     @test det(@SMatrix [0 1; 1 0]) == -1
     @test logdet(@SMatrix Complex{Float64}[0 1; 1 0]) == log(det(@SMatrix Complex{Float64}[0 1; 1 0]))
-
+    @test det(eye(SMatrix{3,3})*im) == det(eye(3,3)*im)
+    
     @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) == -1
     m = randn(Float64, 4,4)
     @test det(SMatrix{4,4}(m)) ≈ det(m)

test/inv.jl

@@ -22,7 +22,10 @@ end
     @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]
-
+    @test inv(@SMatrix [1+im 2 0; 1 2-im 1; 1 1 2+im])::SMatrix ≈ [2-2im  -3+1im   1-1im
+                                                         -1+0im   2+1im  -1+0im
+                                                          0+1im   0-1im   1-0.0im]/2
+   
     m = randn(Float64, 10,10) + eye(10) # well conditioned
     @test inv(SMatrix{10,10}(m))::StaticMatrix ≈ inv(m)
 

test/linalg.jl

@@ -86,6 +86,16 @@ using StaticArrays, Base.Test
         @test @inferred(cross(SVector(UInt(1),UInt(2),UInt(3)), SVector(UInt(4),UInt(5),UInt(6)))) === SVector(-3, 6, -3)
     end
 
+    @testset "inner products" begin
+        v1 = @SVector [2,4,6,8]
+        v2 = @SVector [4,3,2,1]
+
+        @test @inferred(dot(v1, v2)) === 40
+        @test @inferred(dot(v1, -v2)) === -40
+        @test @inferred(dot(v1*im, v2*im)) === 40*im*conj(im)
+        @test @inferred(StaticArrays.bilinear_vecdot(v1*im, v2*im)) === 40*im*im
+    end
+
     @testset "transpose() and conj()" begin
         @test @inferred(conj(SVector(1+im, 2+im))) === SVector(1-im, 2-im)
 
@@ -150,6 +160,7 @@ using StaticArrays, Base.Test
 
     @testset "size zero" begin
         @test vecdot(SVector{0, Float64}(()), SVector{0, Float64}(())) === 0.
+        @test StaticArrays.bilinear_vecdot(SVector{0, Float64}(()), SVector{0, Float64}(())) === 0.
         @test vecnorm(SVector{0, Float64}(())) === 0.
         @test vecnorm(SVector{0, Float64}(()), 1) === 0.
         @test trace(SMatrix{0,0,Float64}(())) === 0.