Merge pull request #66 from JuliaArrays/use-size

Use `Size` more consistently

Author: andyferris

perf/bench9.txt

@@ -0,0 +1,71 @@
+using StaticArrays
+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("=====================================================================")
+        println(" Vectors of length ", N, " and eltype ", T)
+        println("=====================================================================")
+        immutables = [rand(SVector{N,T})]
+        mutables = [rand(T,N), rand(MVector{N,T}), Size(N)(rand(T,N))]
+        instances = vcat(immutables, mutables)
+
+        namelengths = [length(string(typeof(v).name.name)) for v ∈ instances]
+        maxnamelength = maximum(namelengths)
+
+        for v ∈ instances
+            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)) , ")")
+        end
+
+        println()
+
+        for v ∈ mutables
+            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
+
+        println()
+
+        if N > 16
+            continue
+        end
+        println("=====================================================================")
+        println(" Matrices of size ", N, "×", N, " and eltype ", T)
+        println("=====================================================================")
+        immutables = [rand(SMatrix{N,N,T})]
+        mutables = [rand(T,N,N), rand(MMatrix{N,N,T}), Size(N,N)(rand(T,N,N))]
+        instances = vcat(immutables, mutables)
+
+        namelengths = [length(string(typeof(v).name.name)) for v ∈ instances]
+        maxnamelength = maximum(namelengths)
+
+        for m ∈ instances
+            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 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
+
+        println()
+
+    end
+end

perf/benchmark3.jl

@@ -12,21 +12,54 @@ array may be reshaped.
 immutable SizedArray{S,T,N,M} <: StaticArray{T,N}
     data::Array{T,M}
 
-    function SizedArray(a)
+    function SizedArray(a::Array)
         if length(a) != prod(S)
             error("Dimensions $(size(a)) don't match static size $S")
         end
         new(a)
     end
+
+    function SizedArray()
+        new(Array{T,M}(S))
+    end
 end
 
 @inline (::Type{SizedArray{S,T,N}}){S,T,N,M}(a::Array{T,M}) = SizedArray{S,T,N,M}(a)
 @inline (::Type{SizedArray{S,T}}){S,T,M}(a::Array{T,M}) = SizedArray{S,T,_ndims(S),M}(a)
 @inline (::Type{SizedArray{S}}){S,T,M}(a::Array{T,M}) = SizedArray{S,T,_ndims(S),M}(a)
 
+@inline (::Type{SizedArray{S,T,N}}){S,T,N}() = SizedArray{S,T,N,N}()
+@inline (::Type{SizedArray{S,T}}){S,T}() = SizedArray{S,T,_ndims(S),_ndims(S)}()
+
+@generated function (::Type{SizedArray{S,T,N,M}}){S,T,N,M,L}(x::NTuple{L})
+    if L != prod(S)
+        error("Dimension mismatch")
+    end
+    exprs = [:(a[$i] = x[$i]) for i = 1:L]
+    return quote
+        $(Expr(:meta, :inline))
+        a = SizedArray{S,T,N,M}()
+        @inbounds $(Expr(:block, exprs...))
+        return a
+    end
+end
+
+@inline (::Type{SizedArray{S,T,N}}){S,T,N}(x::Tuple) = SizedArray{S,T,N,N}(x)
+@inline (::Type{SizedArray{S,T}}){S,T}(x::Tuple) = SizedArray{S,T,_dims(S),_dims(S)}(x)
+@inline (::Type{SizedArray{S}}){S,T,L}(x::NTuple{L,T}) = SizedArray{S,T,_dims(S),_dims(S)}(x)
+
 # 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
@@ -38,10 +71,18 @@ end
 @propagate_inbounds setindex!(a::SizedArray, v, i::Int) = setindex!(a.data, v, i)
 
 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/SizedArray.jl

@@ -41,6 +41,9 @@ end
 @pure function similar_type{SA<:StaticArray,T}(::Union{SA,Type{SA}}, ::Type{T}, size::Int)
     similar_type(similar_type(SA, T), size)
 end
+@pure function similar_type{SA<:StaticArray,T,S}(::Union{SA,Type{SA}}, ::Type{T}, size::Size{S})
+    similar_type(similar_type(SA, T), size)
+end
 @generated function similar_type{SA<:StaticArray,T}(::Union{SA,Type{SA}}, ::Type{T})
     # This function has a strange error (on tests) regarding double-inference, if it is marked @pure
     if T == eltype(SA)
@@ -115,6 +118,25 @@ end
 
 @pure similar_type{SA<:StaticArray,N}(::Union{SA,Type{SA}}, sizes::Tuple{Vararg{Int,N}}) = SArray{sizes, eltype(SA), N, prod(sizes)}
 
+@generated function similar_type{SA <: StaticArray,S}(::Union{SA,Type{SA}}, ::Size{S})
+    if length(S) == 1
+        return quote
+            $(Expr(:meta, :inline))
+            SVector{$(S[1]), $(eltype(SA))}
+        end
+    elseif length(S) == 2
+        return quote
+            $(Expr(:meta, :inline))
+            SMatrix{$(S[1]), $(S[2]), $(eltype(SA))}
+        end
+    else
+        return quote
+            $(Expr(:meta, :inline))
+            SArray{S, $(eltype(SA)), $(length(S)), $(prod(S))}
+        end
+    end
+end
+
 # Some specializations for the mutable case
 @pure similar_type{MA<:Union{MVector,MMatrix,MArray,SizedArray}}(::Union{MA,Type{MA}}, size::Int) = MVector{size, eltype(MA)}
 @pure similar_type{MA<:Union{MVector,MMatrix,MArray,SizedArray}}(::Union{MA,Type{MA}}, sizes::Tuple{Int}) = MVector{sizes[1], eltype(MA)}
@@ -123,23 +145,73 @@ end
 
 @pure similar_type{MA<:Union{MVector,MMatrix,MArray,SizedArray},N}(::Union{MA,Type{MA}}, sizes::Tuple{Vararg{Int,N}}) = MArray{sizes, eltype(MA), N, prod(sizes)}
 
+@generated function similar_type{MA<:Union{MVector,MMatrix,MArray,SizedArray},S}(::Union{MA,Type{MA}}, ::Size{S})
+    if length(S) == 1
+        return quote
+            $(Expr(:meta, :inline))
+            MVector{$(S[1]), $(eltype(MA))}
+        end
+    elseif length(S) == 2
+        return quote
+            $(Expr(:meta, :inline))
+            MMatrix{$(S[1]), $(S[2]), $(eltype(MA))}
+        end
+    else
+        return quote
+            $(Expr(:meta, :inline))
+            MArray{S, $(eltype(MA)), $(length(S)), $(prod(S))}
+        end
+    end
+end
+
 # And also similar() returning mutable StaticArrays
 @inline similar{SV <: StaticVector}(::SV) = MVector{length(SV),eltype(SV)}()
 @inline similar{SV <: StaticVector, T}(::SV, ::Type{T}) = MVector{length(SV),T}()
-@inline similar{SA <: StaticArray}(::SA, sizes::Tuple{Int}) = MVector{sizes[1], eltype(SA)}()
-@inline similar{SA <: StaticArray}(::SA, size::Int) = MVector{size, eltype(SA)}()
-@inline similar{T}(::StaticArray, ::Type{T}, sizes::Tuple{Int}) = MVector{sizes[1],T}()
-@inline similar{T}(::StaticArray, ::Type{T}, size::Int) = MVector{size,T}()
 
 @inline similar{SM <: StaticMatrix}(m::SM) = MMatrix{size(SM,1),size(SM,2),eltype(SM),length(SM)}()
 @inline similar{SM <: StaticMatrix, T}(::SM, ::Type{T}) = MMatrix{size(SM,1),size(SM,2),T,length(SM)}()
-@inline similar{SA <: StaticArray}(::SA, sizes::Tuple{Int,Int}) = MMatrix{sizes[1], sizes[2], eltype(SA), sizes[1]*sizes[2]}()
-@inline similar(a::StaticArray, T::Type, sizes::Tuple{Int,Int}) = MMatrix{sizes[1], sizes[2], T, sizes[1]*sizes[2]}()
 
 @inline similar{SA <: StaticArray}(m::SA) = MArray{size(SA),eltype(SA),ndims(SA),length(SA)}()
 @inline similar{SA <: StaticArray,T}(m::SA, ::Type{T}) = MArray{size(SA),T,ndims(SA),length(SA)}()
-@inline similar{SA <: StaticArray,N}(m::SA, sizes::NTuple{N, Int}) = MArray{sizes,eltype(SA),N,prod(sizes)}()
-@inline similar{SA <: StaticArray,N,T}(m::SA, ::Type{T}, sizes::NTuple{N, Int}) = MArray{sizes,T,N,prod(sizes)}()
+
+@generated function similar{SA <: StaticArray,S}(::SA, ::Size{S})
+    if length(S) == 1
+        return quote
+            $(Expr(:meta, :inline))
+            MVector{$(S[1]), $(eltype(SA))}()
+        end
+    elseif length(S) == 2
+        return quote
+            $(Expr(:meta, :inline))
+            MMatrix{$(S[1]), $(S[2]), $(eltype(SA))}()
+        end
+    else
+        return quote
+            $(Expr(:meta, :inline))
+            MArray{S, $(eltype(SA))}()
+        end
+    end
+end
+
+@generated function similar{SA <: StaticArray, T, S}(::SA, ::Type{T}, ::Size{S})
+    if length(S) == 1
+        return quote
+            $(Expr(:meta, :inline))
+            MVector{$(S[1]), T}()
+        end
+    elseif length(S) == 2
+        return quote
+            $(Expr(:meta, :inline))
+            MMatrix{$(S[1]), $(S[2]), T}()
+        end
+    else
+        return quote
+            $(Expr(:meta, :inline))
+            MArray{S, T}()
+        end
+    end
+end
+
 
 # This is used in Base.LinAlg quite a lot, and it impacts type stability
 # since some functions like expm() branch on a check for Hermitian or Symmetric

src/abstractarray.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)))))