Use Size in eig, chol, inv...

Author: Andy Ferris

perf/benchmark3.jl

@@ -3,6 +3,13 @@ using BenchmarkTools
 
 import BenchmarkTools: prettytime, prettymemory
 
+@noinline plus(a,b) = a+b
+@noinline plus!(c,a,b) = broadcast!(+, c, a, b)
+
+@noinline mul(a,b) = a*b
+@noinline mul!(c,a,b) = A_mul_B!(c, a, b)
+
+
 for T ∈ [Int64, Float64]
     for N ∈ [1,2,4,8,16,32,64,128,256]
         println("=====================================================================")
@@ -16,15 +23,15 @@ for T ∈ [Int64, Float64]
         maxnamelength = maximum(namelengths)
 
         for v ∈ instances
-            result = mean(@benchmark +($(copy(v)), $(copy(v))))
+            result = mean(@benchmark plus($(copy(v)), $(copy(v))))
             padding = maxnamelength - length(string(typeof(v).name.name))
-            println(typeof(v).name.name, ":", " " ^ padding, " v3 = v1 + v2 takes ", prettytime(time(result)), ", ", prettymemory(memory(result)), "(GC ", prettytime(gctime(result)) , ")")
+            println(typeof(v).name.name, ":", " " ^ padding, " v3 = v1 + v2 takes ", prettytime(time(result)), ", ", prettymemory(memory(result)), " (GC ", prettytime(gctime(result)) , ")")
         end
 
         println()
 
         for v ∈ mutables
-            result = mean(@benchmark broadcast!(+, $(copy(v)), $(copy(v)), $(copy(v))))
+            result = mean(@benchmark plus!($(copy(v)), $(copy(v)), $(copy(v))))
             padding = maxnamelength - length(string(typeof(v).name.name))
             println(typeof(v).name.name, ":", " " ^ padding, " v3 .= +.(v1, v2) takes ", prettytime(time(result)), ", ", prettymemory(memory(result)), " (GC ", prettytime(gctime(result)) , ")")
         end
@@ -45,15 +52,15 @@ for T ∈ [Int64, Float64]
         maxnamelength = maximum(namelengths)
 
         for m ∈ instances
-            result = mean(@benchmark *($(copy(m)), $(copy(m))))
+            result = mean(@benchmark mul($(copy(m)), $(copy(m))))
             padding = maxnamelength - length(string(typeof(m).name.name))
             println(typeof(m).name.name, ":", " " ^ padding, " m3 = m1 * m2 takes ", prettytime(time(result)), ", ", prettymemory(memory(result)), " (GC ", prettytime(gctime(result)) , ")")
         end
 
         println()
 
         for m ∈ mutables
-            result = mean(@benchmark A_mul_B!($(copy(m)), $(copy(m)), $(copy(m))))
+            result = mean(@benchmark mul!($(copy(m)), $(copy(m)), $(copy(m))))
             padding = maxnamelength - length(string(typeof(m).name.name))
             println(typeof(m).name.name, ":", " " ^ padding, " A_mul_B!(m3, m1, m2) takes ", prettytime(time(result)), ", ", prettymemory(memory(result)), " (GC ", prettytime(gctime(result)) , ")")
         end

src/SizedArray.jl

@@ -51,6 +51,15 @@ end
 # Overide some problematic default behaviour
 @inline convert{SA<:SizedArray}(::Type{SA}, sa::SizedArray) = SA(sa.data)
 
+# Back to Array (unfortunately need both convert and construct to overide other methods)
+@inline (::Type{Array})(sa::SizedArray) = sa.data
+@inline (::Type{Array{T}}){T,S}(sa::SizedArray{S,T}) = sa.data
+@inline (::Type{Array{T,N}}){T,S,N}(sa::SizedArray{S,T,N}) = sa.data
+
+@inline convert(::Type{Array}, sa::SizedArray) = sa.data
+@inline convert{T,S}(::Type{Array{T}}, sa::SizedArray{S,T}) = sa.data
+@inline convert{T,S,N}(::Type{Array{T,N}}, sa::SizedArray{S,T,N}) = sa.data
+
 @pure _ndims{N}(::NTuple{N,Int}) = N
 
 @pure size{S}(::Type{SizedArray{S}}) = S
@@ -65,11 +74,15 @@ typealias SizedVector{S,T,M} SizedArray{S,T,1,M}
 @pure size{S}(::Type{SizedVector{S}}) = S
 @inline (::Type{SizedVector{S}}){S,T,M}(a::Array{T,M}) = SizedArray{S,T,1,M}(a)
 @inline (::Type{SizedVector{S}}){S,T,L}(x::NTuple{L,T}) = SizedArray{S,T,1,1}(x)
+@inline (::Type{Vector})(sa::SizedVector) = sa.data
+@inline convert(::Type{Vector}, sa::SizedVector) = sa.data
 
 typealias SizedMatrix{S,T,M} SizedArray{S,T,2,M}
 @pure size{S}(::Type{SizedMatrix{S}}) = S
 @inline (::Type{SizedMatrix{S}}){S,T,M}(a::Array{T,M}) = SizedArray{S,T,2,M}(a)
 @inline (::Type{SizedMatrix{S}}){S,T,L}(x::NTuple{L,T}) = SizedArray{S,T,2,2}(x)
+@inline (::Type{Matrix})(sa::SizedMatrix) = sa.data
+@inline convert(::Type{Matrix}, sa::SizedMatrix) = sa.data
 
 
 """

src/cholesky.jl

@@ -1,21 +1,20 @@
 # Generic Cholesky decomposition for fixed-size matrices, mostly unrolled
+@inline function Base.chol(A::StaticMatrix)
+    ishermitian(A) || Base.LinAlg.non_hermitian_error("chol")
+    _chol(Size(A), A)
+end
 
-# Currently all sanity checks are disabled!
-@generated function Base.chol(A::StaticMatrix)
-    if size(A) === (1, 1)
-        return :(_chol1(A))
-    elseif size(A) === (2, 2)
-        #ishermitian(A) || Base.LinAlg.non_hermitian_error("chol")
-        return :(_chol2(A))
-    elseif size(A) === (3, 3)
-        #ishermitian(A) || Base.LinAlg.non_hermitian_error("chol")
-        return :(_chol3(A))
-    else
-        return :(chol(Array(A)))
-    end
+@inline function Base.chol{T<:Real, SM <: StaticMatrix}(A::Base.LinAlg.RealHermSymComplexHerm{T,SM})
+    ishermitian(A) || Base.LinAlg.non_hermitian_error("chol")
+    _chol(Size(A), A)
+end
+
+@inline function Base.chol{SM<:StaticMatrix}(A::Symmetric{SM})
+    eltype(A) <: Real && (ishermitian(A) || Base.LinAlg.non_hermitian_error("chol"))
+    _chol(Size(A), A)
 end
 
-@generated function _chol1(A::StaticMatrix)
+@generated function _chol(::Size{(1,1)}, A::StaticMatrix)
     @assert size(A) == (1,1)
     T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
     newtype = similar_type(A,T)
@@ -26,8 +25,7 @@ end
     end
 end
 
-
-@generated function _chol2(A::StaticMatrix)
+@generated function _chol(::Size{(2,2)}, A::StaticMatrix)
     @assert size(A) == (2,2)
     T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
     newtype = similar_type(A,T)
@@ -41,7 +39,7 @@ end
     end
 end
 
-@generated function _chol3(A::StaticMatrix)
+@generated function _chol(::Size{(3,3)}, A::StaticMatrix)
     @assert size(A) == (3,3)
     T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
     newtype = similar_type(A,T)
@@ -57,3 +55,6 @@ end
         ($newtype)((a11, $(zero(T)), $(zero(T)), a12, a22, $(zero(T)), a13, a23, a33))
     end
 end
+
+# Otherwise default algorithm returning wrapped SizedArray
+@inline _chol(s::Size, A::StaticArray) = s(full(chol(Hermitian(Array(A)))))