Merge pull request #208

Author: c42f

src/matrix_multiply.jl

@@ -1,13 +1,10 @@
 import Base: *,        Ac_mul_B,  A_mul_Bc,  Ac_mul_Bc,  At_mul_B,  A_mul_Bt,  At_mul_Bt
 import Base: A_mul_B!, Ac_mul_B!, A_mul_Bc!, Ac_mul_Bc!, At_mul_B!, A_mul_Bt!, At_mul_Bt!
 
-import Base.LinAlg: BlasFloat
+import Base.LinAlg: BlasFloat, matprod
 
 const StaticVecOrMat{T} = Union{StaticVector{<:Any, T}, StaticMatrix{<:Any, <:Any, T}}
 
-# Idea inspired by https://github.com/JuliaLang/julia/pull/18218
-promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero(T1)*zero(T2))
-
 # TODO Potentially a loop version for rather large arrays? Or try and figure out inference problems?
 
 # Deal with A_mul_Bc, etc...
@@ -60,7 +57,7 @@ promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero
         if length(b) != sa[2]
             throw(DimensionMismatch("Tried to multiply arrays of size $sa and $(size(b))"))
         end
-        T = promote_matprod(Ta, Tb)
+        T = promote_op(matprod,Ta,Tb)
         @inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...)))
     end
 end
@@ -78,7 +75,7 @@ end
 
     return quote
         @_inline_meta
-        T = promote_matprod(Ta, Tb)
+        T = promote_op(matprod,Ta,Tb)
         @inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...)))
     end
 end
@@ -171,7 +168,7 @@ end
 
     return quote
         @_inline_meta
-        T = promote_matprod(Ta, Tb)
+        T = promote_op(matprod,Ta,Tb)
         @inbounds return similar_type(a, T, $S)(tuple($(exprs...)))
     end
 end
@@ -190,7 +187,7 @@ end
 
     return quote
         @_inline_meta
-        T = promote_matprod(Ta, Tb)
+        T = promote_op(matprod,Ta,Tb)
 
         @inbounds $(Expr(:block, exprs_init...))
         for j = 2:$(sa[2])
@@ -218,7 +215,7 @@ end
 
     return quote
         @_inline_meta
-        T = promote_matprod(Ta, Tb)
+        T = promote_op(matprod,Ta,Tb)
         $(Expr(:block,
             vect_exprs...,
             :(@inbounds return similar_type(a, T, $S)(tuple($(exprs...))))
@@ -234,7 +231,7 @@ end
     if sa[2] != 0
         exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(a[$(sub2ind(sa, k, j))]*b[$j]) for j = 1:sa[2]]) for k = 1:sa[1]]
     else
-        exprs = [:(zero(promote_matprod(Ta,Tb))) for k = 1:sa[1]]
+        exprs = [:(zero(promote_op(matprod,Ta,Tb))) for k = 1:sa[1]]
     end
 
     return quote

test/matrix_multiply.jl

@@ -11,6 +11,11 @@
         @test isa(x, SVector{2,CartesianIndex{2}})
         @test x == @SVector [CartesianIndex((7,5)), CartesianIndex((15,13))]
 
+        # block matrices
+        bm = @SMatrix [m m; m m]
+        bv = @SVector [v,v]
+        @test (bm*bv)::SVector{2,SVector{2,Int}} == @SVector [[10,22],[10,22]]
+
         # inner product
         @test @inferred(v'*v) === 5
 
@@ -49,6 +54,13 @@
         v = @SVector [1, 2]
         @test @inferred(v*m) === @SMatrix [1 2 3 4; 2 4 6 8]
 
+        # block matrices
+        m = @SMatrix [1 2; 3 4]
+        bm = @SMatrix [m m; m m]
+        bv = @SVector [v,v]
+        # Broken only because output turns into a normal array:
+        @test_broken (bv'*bm)'::SVector{2,SVector{2,Int}} == @SVector [[14,20],[14,20]]
+
         # Outer product
         v2 = SVector(1, 2)
         v3 = SVector(3, 4)
@@ -87,6 +99,11 @@
         n = @MArray [2 3; 4 5]
         @test (m*n) == @SMatrix [10 13; 22 29]
 
+        # block matrices
+        bm = @SMatrix [m m; m m]
+        bm2 = @SMatrix [14 20; 30 44]
+        @test (bm*bm)::SMatrix{2,2,SMatrix{2,2,Int,4}} == @SMatrix [bm2 bm2; bm2 bm2]
+
         # Alternative methods used between 8 < n <= 14 and n > 14
         m_array = rand(1:10, 10, 10)
         n_array = rand(1:10, 10, 10)