From bff627a60369a2a3e0be6aa7d9718ff07236fc72 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Thu, 23 Feb 2017 22:42:06 +1000 Subject: [PATCH 1/9] WIP v0.6 support --- src/StaticArrays.jl | 12 +++++-- src/abstractarray.jl | 23 +++--------- src/core.jl | 85 +++++++++++++++++++++++--------------------- src/generated.jl | 15 ++++++++ src/indexing.jl | 21 +++++------ src/traits.jl | 47 ++++++++++++++++++++---- 6 files changed, 124 insertions(+), 79 deletions(-) create mode 100644 src/generated.jl diff --git a/src/StaticArrays.jl b/src/StaticArrays.jl index bfe53906..4f17b40b 100644 --- a/src/StaticArrays.jl +++ b/src/StaticArrays.jl @@ -2,7 +2,9 @@ __precompile__() module StaticArrays -import Base: @pure, @propagate_inbounds, getindex, setindex!, size, similar, +import Base: @_inline_meta, @_propagate_inbounds_meta, @_pure_meta, @propagate_inbounds, @pure + +import Base: getindex, setindex!, size, similar, length, convert, promote_op, map, map!, reduce, reducedim, mapreduce, broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill, @@ -18,7 +20,7 @@ export MArray, MVector, MMatrix export FieldVector, MutableFieldVector export SizedArray, SizedVector, SizedMatrix -export Size +export Size, Length export @SVector, @SMatrix, @SArray export @MVector, @MMatrix, @MArray @@ -27,6 +29,7 @@ export similar_type export push, pop, shift, unshift, insert, deleteat, setindex include("util.jl") +include("generated.jl") include("core.jl") include("traits.jl") @@ -40,8 +43,9 @@ include("MMatrix.jl") include("MArray.jl") include("SizedArray.jl") + include("abstractarray.jl") -include("indexing.jl") +include("indexing.jl") #= include("mapreduce.jl") include("arraymath.jl") include("linalg.jl") @@ -52,6 +56,8 @@ include("det.jl") include("inv.jl") include("eigen.jl") include("cholesky.jl") +=# + if VERSION < v"0.6.0-dev.1671" include("FixedSizeArrays.jl") diff --git a/src/abstractarray.jl b/src/abstractarray.jl index 8d589f23..c4b4b062 100644 --- a/src/abstractarray.jl +++ b/src/abstractarray.jl @@ -1,7 +1,8 @@ -@compat StaticScalar{T} = StaticArray{T,0} +length(a::T) where {T <: StaticArray} = prod(Size(T)) +length(a::Type{T}) where {T<:StaticArray} = prod(Size(T)) -length{T<:StaticArray}(a::T) = prod(Size(T)) -length{T<:StaticArray}(a::Type{T}) = prod(Size(T)) +length_val(a::T) where {T <: StaticArray} = length_val(Size(T)) +length_val(a::Type{T}) where {T<:StaticArray} = length_val(Size(T)) size{T<:StaticArray}(::T) = get(Size(T)) size{T<:StaticArray}(::Type{T}) = get(Size(T)) @@ -9,20 +10,6 @@ size{T<:StaticArray}(::Type{T}) = get(Size(T)) size{T<:StaticArray}(::T, d::Integer) = Size(T)[d] size{T<:StaticArray}(::Type{T}, d::Integer) = Size(T)[d] -# This has to be defined after length and size because it is generated -@generated function convert{SA<:StaticArray}(::Type{SA}, a::AbstractArray) - L = length(SA) - exprs = [:(a[$i]) for i = 1:L] - - return quote - $(Expr(:meta, :inline)) - if length(a) != $L - L = $L - error("Dimension mismatch. Expected input array of length $L, got length $(length(a))") - end - @inbounds return SA($(Expr(:tuple, exprs...))) - end -end # This seems to confuse Julia a bit in certain circumstances (specifically for trailing 1's) @inline function Base.isassigned(a::StaticArray, i::Int...) @@ -30,7 +17,7 @@ end 1 <= ii <= length(a) ? true : false end -@compat Base.IndexStyle{T<:StaticArray}(::Type{T}) = IndexLinear() +Base.IndexStyle{T<:StaticArray}(::Type{T}) = IndexLinear() # Default type search for similar_type """ diff --git a/src/core.jl b/src/core.jl index 4f38d852..5bd371a0 100644 --- a/src/core.jl +++ b/src/core.jl @@ -1,5 +1,6 @@ """ - abstract StaticArray{T, N} <: AbstractArray{T, N} + abstract type StaticArray{T, N} <: AbstractArray{T, N} end + StaticScalar{T} = StaticArray{T, 0} StaticVector{T} = StaticArray{T, 1} StaticMatrix{T} = StaticArray{T, 2} @@ -27,42 +28,53 @@ For mutable containers you may also need to define the following: (see also `SVector`, `SMatrix`, `SArray`, `MVector`, `MMatrix`, `MArray`, `SizedArray` and `FieldVector`) """ -@compat abstract type StaticArray{T, N} <: AbstractArray{T, N} end +abstract type StaticArray{T, N} <: AbstractArray{T, N} end -@compat StaticVector{T} = StaticArray{T, 1} -@compat StaticMatrix{T} = StaticArray{T, 2} +StaticScalar{T} = StaticArray{T, 1} +StaticVector{T} = StaticArray{T, 1} +StaticMatrix{T} = StaticArray{T, 2} -# People might not want to use Tuple for everything (TODO: check this with FieldVector...) -# Generic case, with least 2 inputs -@inline (::Type{SA}){SA<:StaticArray}(x1,x2,xs...) = SA((x1,x2,xs...)) - -@inline convert{SA<:StaticArray}(::Type{SA}, x::Tuple) = error("No precise constructor found. Length of input was $(length(x)) while length of $SA is $(length(SA)).") - -# Avoiding splatting penalties. Being here, implementations of StaticArray will not have to deal with these. TODO check these are necessary or not -#@inline (::Type{SA}){SA<:StaticArray}(x1) = SA((x1,)) # see convert below (lesser precedence than other constructors?) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2) = SA((x1,x2)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3) = SA((x1,x2,x3)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4) = SA((x1,x2,x3,x4)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5) = SA((x1,x2,x3,x4,x5)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6) = SA((x1,x2,x3,x4,x5,x6)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7) = SA((x1,x2,x3,x4,x5,x6,x7)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8) = SA((x1,x2,x3,x4,x5,x6,x7,x8)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15)) -@inline convert{SA<:StaticArray}(::Type{SA},x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16) = SA((x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16)) - -@inline convert{SA<:StaticArray}(::Type{SA}, x1) = SA((x1,)) +(::Type{SA})(x::Tuple) where {SA <: StaticArray} = error("No precise constructor for $SA found. Length of input was $(length(x)).") + +@inline convert(::Type{SA}, x...) where {SA <: StaticArray} = SA(x) # this covers most conversions and "statically-sized reshapes" -@inline convert{SA<:StaticArray}(::Type{SA}, sa::StaticArray) = SA(Tuple(sa)) +@inline convert(::Type{SA}, sa::StaticArray) where {SA<:StaticArray} = SA(Tuple(sa)) +@inline convert(::Type{SA}, sa::SA) where {SA<:StaticArray} = sa + +# A general way of going back to a tuple +@inline function convert(::Type{Tuple}, a::StaticArray) + unroll_tuple((i -> @inbounds return a[i]), length_val(a)) +end + +@inline function convert(::Type{SA}, a::AbstractArray) where {SA <: StaticArray} + if length(a) != length(SA) + error("Dimension mismatch. Expected input array of length $(length(SA)), got length $(length(a))") + end + + return SA(unroll_tuple((i -> @inbounds return a[i]), length_val(SA))) +end -@inline convert{SA<:StaticArray}(::Type{SA}, sa::SA) = sa +#= +@generated function convert(::Type{Tuple}, a::StaticArray) + n = length(a) + exprs = [:(a[$j]) for j = 1:n] + quote + $(Expr(:meta, :inline)) + @inbounds return $(Expr(:tuple, exprs...)) + end +end +=# + + + +# People might not want to use Tuple for everything (TODO: check this with FieldVector...) +# Generic case, with least 2 inputs +#@inline (::Type{SA}){SA<:StaticArray}(x1,x2,xs...) = SA((x1,x2,xs...)) + + +#= function convert{T,N}(::Type{Array}, sa::StaticArray{T,N}) out = Array{T,N}(size(sa)) @inbounds for i = 1:length(sa) @@ -102,13 +114,4 @@ function convert{T}(::Type{Vector}, sa::StaticVector{T}) end return out end - -# A general way of going back to a tuple, etc -@generated function convert(::Type{Tuple}, a::StaticArray) - n = length(a) - exprs = [:(a[$j]) for j = 1:n] - quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:tuple, exprs...)) - end -end +=# diff --git a/src/generated.jl b/src/generated.jl new file mode 100644 index 00000000..36f65ab6 --- /dev/null +++ b/src/generated.jl @@ -0,0 +1,15 @@ +#= +""" + unroll_tuple(f, Val{N}) + +Like the `Base.ntuple` function, but generates fast code beyond N = 16. +""" +=# +@generated function unroll_tuple(f, ::Type{Val{N}}) where {N} + @assert(N isa Int) + exprs = [:(f($i)) for i = 1:N] + quote + @_propagate_inbounds_meta + return $(Expr(:tuple, exprs...)) + end +end diff --git a/src/indexing.jl b/src/indexing.jl index 36653eab..43c51918 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -2,18 +2,15 @@ ## Linear Indexing ## ###################### -# What to do about @boundscheck and @inbounds? It's worse sometimes than @inline, for tuples... -@generated function getindex{SA<:StaticArray, S}(a::SA, inds::NTuple{S,Integer}) - newtype = similar_type(SA, Size(S)) - exprs = [:(a[inds[$i]]) for i = 1:S] - - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end +# Static Vector indexing into AbstractArrays +@propagate_inbounds function getindex(a::AbstractArray{T}, inds::StaticVector{<:Integer}) where {T} + similar_type(a, T, Size(inds))(unroll_tuple(Length(inds)) do i + @_propagate_inbounds_meta + a[inds[i]] + end) end -# Static Vector indexing into AbstractArrays +#= @generated function getindex{T, I <: Integer}( a::AbstractArray{T}, inds::StaticVector{I} ) @@ -25,6 +22,8 @@ end return $(Expr(:call, newtype, Expr(:tuple, exprs...))) end end + + # Convert to StaticArrays using tuples # TODO think about bounds checks here. @generated function getindex{T, I <: Integer}( @@ -434,6 +433,8 @@ end @inbounds return a[$ind_expr] = val end end +=# + # TODO higher dimensional versions with tuples and : diff --git a/src/traits.jl b/src/traits.jl index d712e7d1..79712fa6 100644 --- a/src/traits.jl +++ b/src/traits.jl @@ -21,19 +21,20 @@ _det(::Size{(2,2)}, x::StaticMatrix) = x[1,1]*x[2,2] - x[1,2]*x[2,1] ``` """ immutable Size{S} - function (::Type{Size{S}}){S}() + function Size{S}() where S check_size(S) new{S}() end end + check_size(S::Tuple{Vararg{Int}}) = nothing check_size(S) = error("Size was expected to be a tuple of `Int`s") @pure Size(s::Tuple{Vararg{Int}}) = Size{s}() @pure Size(s::Int...) = Size{s}() -@inline Size(a::StaticArray) = Size(typeof(a)) +Size(a::StaticArray) = Size(typeof(a)) Base.show{S}(io::IO, ::Size{S}) = print(io, "Size", S) @@ -53,13 +54,32 @@ function Size{SA <: StaticArray}(::Type{SA}) """) end -# Some @pure convenience functions. -@pure get{S}(::Size{S}) = S -@pure getindex{S}(::Size{S}, i::Int) = i <= length(S) ? S[i] : 1 +struct Length{L} + function Length{L}() where L + check_length(L) + new{L}() + end +end + +check_length(L::Int) = nothing +check_length(L) = error("Length was expected to be an `Int`") + +Base.show(io::IO, ::Length{L}) where {L} = print(io, "Length(", L, ")") + +Length(a::StaticArray) = Length(Size(a)) +Length(::Type{SA}) where {SA <: StaticArray} = Length(Size(SA)) +@pure Length(::Size{S}) where {S} = Length{prod(S)}() +@pure Length(L::Int) = Length{L}() + -@pure length{S}(::Size{S}) = length(S) -@generated length_val{S}(::Size{S}) = Val{length(S)} +# Some @pure convenience functions for `Size` +@pure get(::Size{S}) where {S} = S + +@pure getindex(::Size{S}, i::Int) where {S} = i <= length(S) ? S[i] : 1 + +@pure length(::Size{S}) where {S} = length(S) +@pure length_val{S}(::Size{S}) = Val{length(S)} @pure Base.:(==){S}(::Size{S}, s::Tuple{Vararg{Int}}) = S == s @pure Base.:(==){S}(s::Tuple{Vararg{Int}}, ::Size{S}) = s == S @@ -70,3 +90,16 @@ end @pure Base.prod{S}(::Size{S}) = prod(S) @pure @inline Base.sub2ind{S}(::Size{S}, x::Int...) = sub2ind(S, x...) + +# Some @pure convenience functions for `Length` +@pure get(::Length{L}) where {L} = L + +@pure Base.:(==)(::Length{L}, l::Int) where {L} = L == l +@pure Base.:(==)(l::Int, ::Length{L}) where {L} = l == L + +@pure Base.:(!=)(::Length{L}, l::Int) where {L} = L != l +@pure Base.:(!=)(l::Int, ::Length{L}) where {L} = l != L + + +# The generated functions work with length, etc... +@propagate_inbounds unroll_tuple(f, ::Length{L}) where {L} = unroll_tuple(f, Val{L}) From a1eab9240c0063f82c339c41f31b635130545c67 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Sun, 26 Feb 2017 22:17:50 +1000 Subject: [PATCH 2/9] More WIP for Julia v0.6 --- REQUIRE | 1 - src/FieldVector.jl | 6 +- src/ImmutableArrays.jl | 51 +- src/MArray.jl | 14 +- src/MMatrix.jl | 12 +- src/MVector.jl | 6 +- src/SArray.jl | 10 +- src/SMatrix.jl | 10 +- src/SUnitRange.jl | 12 + src/SVector.jl | 2 +- src/Scalar.jl | 10 +- src/SizedArray.jl | 4 +- src/StaticArrays.jl | 77 ++- src/abstractarray.jl | 149 +++--- src/arraymath.jl | 239 +++------ src/broadcast.jl | 147 ++++++ src/{core.jl => convert.jl} | 53 +- src/generated.jl | 15 - src/indexing.jl | 937 ++++++++---------------------------- src/linalg.jl | 324 +++++-------- src/mapreduce.jl | 506 ++++++------------- src/matrix_multiply.jl | 753 +++++++---------------------- src/util.jl | 8 +- test/abstractarray.jl | 20 +- test/arraymath.jl | 2 - test/indexing.jl | 36 +- test/linalg.jl | 12 +- test/mapreduce.jl | 4 +- test/matrix_multiply.jl | 25 +- test/runtests.jl | 7 +- 30 files changed, 1136 insertions(+), 2316 deletions(-) create mode 100644 src/SUnitRange.jl create mode 100644 src/broadcast.jl rename src/{core.jl => convert.jl} (55%) delete mode 100644 src/generated.jl diff --git a/REQUIRE b/REQUIRE index b6cc633c..94237c0f 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1 @@ julia 0.5 -Compat 0.19.0 diff --git a/src/FieldVector.jl b/src/FieldVector.jl index b7beee75..289798ba 100644 --- a/src/FieldVector.jl +++ b/src/FieldVector.jl @@ -15,15 +15,15 @@ For example: z::Float64 end """ -@compat abstract type FieldVector{T} <: StaticVector{T} end +abstract type FieldVector{T} <: StaticVector{T} end # Is this a good idea?? Should people just define constructors that accept tuples? @inline (::Type{FV}){FV<:FieldVector}(x::Tuple) = FV(x...) @pure Size{FV<:FieldVector}(::Type{FV}) = Size(nfields(FV)) -@inline getindex(v::FieldVector, i::Integer) = getfield(v, i) -@inline setindex!(v::FieldVector, x, i::Integer) = setfield!(v, i, x) +@inline getindex(v::FieldVector, i::Int) = getfield(v, i) +@inline setindex!(v::FieldVector, x, i::Int) = setfield!(v, i, x) # See #53 Base.cconvert{T}(::Type{Ptr{T}}, v::FieldVector) = Ref(v) diff --git a/src/ImmutableArrays.jl b/src/ImmutableArrays.jl index f3497677..402e3630 100644 --- a/src/ImmutableArrays.jl +++ b/src/ImmutableArrays.jl @@ -1,32 +1,31 @@ module ImmutableArrays using ..StaticArrays -using Compat - -@compat Vector1{T} = SVector{1,T} -@compat Vector2{T} = SVector{2,T} -@compat Vector3{T} = SVector{3,T} -@compat Vector4{T} = SVector{4,T} - -@compat Matrix1x1{T} = SMatrix{1,1,T,1} -@compat Matrix1x2{T} = SMatrix{1,2,T,2} -@compat Matrix1x3{T} = SMatrix{1,3,T,3} -@compat Matrix1x4{T} = SMatrix{1,4,T,4} - -@compat Matrix2x1{T} = SMatrix{2,1,T,2} -@compat Matrix2x2{T} = SMatrix{2,2,T,4} -@compat Matrix2x3{T} = SMatrix{2,3,T,6} -@compat Matrix2x4{T} = SMatrix{2,4,T,8} - -@compat Matrix3x1{T} = SMatrix{3,1,T,3} -@compat Matrix3x2{T} = SMatrix{3,2,T,6} -@compat Matrix3x3{T} = SMatrix{3,3,T,9} -@compat Matrix3x4{T} = SMatrix{3,4,T,12} - -@compat Matrix4x1{T} = SMatrix{4,1,T,4} -@compat Matrix4x2{T} = SMatrix{4,2,T,8} -@compat Matrix4x3{T} = SMatrix{4,3,T,12} -@compat Matrix4x4{T} = SMatrix{4,4,T,16} + +Vector1{T} = SVector{1,T} +Vector2{T} = SVector{2,T} +Vector3{T} = SVector{3,T} +Vector4{T} = SVector{4,T} + +Matrix1x1{T} = SMatrix{1,1,T,1} +Matrix1x2{T} = SMatrix{1,2,T,2} +Matrix1x3{T} = SMatrix{1,3,T,3} +Matrix1x4{T} = SMatrix{1,4,T,4} + +Matrix2x1{T} = SMatrix{2,1,T,2} +Matrix2x2{T} = SMatrix{2,2,T,4} +Matrix2x3{T} = SMatrix{2,3,T,6} +Matrix2x4{T} = SMatrix{2,4,T,8} + +Matrix3x1{T} = SMatrix{3,1,T,3} +Matrix3x2{T} = SMatrix{3,2,T,6} +Matrix3x3{T} = SMatrix{3,3,T,9} +Matrix3x4{T} = SMatrix{3,4,T,12} + +Matrix4x1{T} = SMatrix{4,1,T,4} +Matrix4x2{T} = SMatrix{4,2,T,8} +Matrix4x3{T} = SMatrix{4,3,T,12} +Matrix4x4{T} = SMatrix{4,4,T,16} export Vector1, Vector2, Vector3, Vector4, Matrix1x1, Matrix1x2, Matrix1x3, Matrix1x4, diff --git a/src/MArray.jl b/src/MArray.jl index c8754aee..8f223982 100644 --- a/src/MArray.jl +++ b/src/MArray.jl @@ -85,10 +85,10 @@ end @inline MArray(a::StaticArray) = MArray{size(typeof(a))}(Tuple(a)) # Some more advanced constructor-like functions -@inline eye{Size}(::Type{MArray{Size}}) = eye(MArray{Size,Float64}) -@inline zeros{Size}(::Type{MArray{Size}}) = zeros(MArray{Size,Float64}) -@inline ones{Size}(::Type{MArray{Size}}) = ones(MArray{Size,Float64}) - +@inline one(::Type{MArray{S}}) where {S} = one(MArray{S,Float64,length(S)}) +@inline eye(::Type{MArray{S}}) where {S} = eye(MArray{S,Float64,length(S)}) +@inline one(::Type{MArray{S,T}}) where {S,T} = one(MArray{S,T,length(S)}) +@inline eye(::Type{MArray{S,T}}) where {S,T} = eye(MArray{S,T,length(S)}) #################### ## MArray methods ## @@ -99,13 +99,13 @@ end @pure Size{S,T,N}(::Type{MArray{S,T,N}}) = Size(S) @pure Size{S,T,N,L}(::Type{MArray{S,T,N,L}}) = Size(S) -function getindex(v::MArray, i::Integer) +function getindex(v::MArray, i::Int) Base.@_inline_meta v.data[i] end -@propagate_inbounds setindex!{S,T}(v::MArray{S,T}, val, i::Integer) = setindex!(v, convert(T, val), i) -@inline function setindex!{S,T}(v::MArray{S,T}, val::T, i::Integer) +@propagate_inbounds setindex!{S,T}(v::MArray{S,T}, val, i::Int) = setindex!(v, convert(T, val), i) +@inline function setindex!{S,T}(v::MArray{S,T}, val::T, i::Int) @boundscheck if i < 1 || i > length(v) throw(BoundsError()) end diff --git a/src/MMatrix.jl b/src/MMatrix.jl index cce2fa3e..6e34f0ae 100644 --- a/src/MMatrix.jl +++ b/src/MMatrix.jl @@ -94,9 +94,9 @@ end # Some more advanced constructor-like functions @inline one{N}(::Type{MMatrix{N}}) = one(MMatrix{N,N}) @inline eye{N}(::Type{MMatrix{N}}) = eye(MMatrix{N,N}) -@inline eye{N,M}(::Type{MMatrix{N,M}}) = eye(MMatrix{N,M,Float64}) -@inline zeros{N,M}(::Type{MMatrix{N,M}}) = zeros(MMatrix{N,M,Float64}) -@inline ones{N,M}(::Type{MMatrix{N,M}}) = ones(MMatrix{N,M,Float64}) +#@inline eye{N,M}(::Type{MMatrix{N,M}}) = eye(MMatrix{N,M,Float64}) +#@inline zeros{N,M}(::Type{MMatrix{N,M}}) = zeros(MMatrix{N,M,Float64}) +#@inline ones{N,M}(::Type{MMatrix{N,M}}) = ones(MMatrix{N,M,Float64}) ##################### ## MMatrix methods ## @@ -106,7 +106,7 @@ end @pure Size{S1,S2,T}(::Type{MMatrix{S1,S2,T}}) = Size(S1, S2) @pure Size{S1,S2,T,L}(::Type{MMatrix{S1,S2,T,L}}) = Size(S1, S2) -@propagate_inbounds function getindex{S1,S2,T}(m::MMatrix{S1,S2,T}, i::Integer) +@propagate_inbounds function getindex{S1,S2,T}(m::MMatrix{S1,S2,T}, i::Int) #@boundscheck if i < 1 || i > length(m) # throw(BoundsError(m,i)) #end @@ -121,8 +121,8 @@ end end end -@propagate_inbounds setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val, i::Integer) = setindex!(m, convert(T, val), i) -@propagate_inbounds function setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val::T, i::Integer) +@propagate_inbounds setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val, i::Int) = setindex!(m, convert(T, val), i) +@propagate_inbounds function setindex!{S1,S2,T}(m::MMatrix{S1,S2,T}, val::T, i::Int) #@boundscheck if i < 1 || i > length(m) # throw(BoundsError(m,i)) #end diff --git a/src/MVector.jl b/src/MVector.jl index a65e2ff1..704e9f3e 100644 --- a/src/MVector.jl +++ b/src/MVector.jl @@ -49,13 +49,13 @@ end @pure Size{S}(::Type{MVector{S}}) = Size(S) @pure Size{S,T}(::Type{MVector{S,T}}) = Size(S) -@propagate_inbounds function getindex(v::MVector, i::Integer) +@propagate_inbounds function getindex(v::MVector, i::Int) v.data[i] end # Mutating setindex! -@propagate_inbounds setindex!{S,T}(v::MVector{S,T}, val, i::Integer) = setindex!(v, convert(T, val), i) -@inline function setindex!{S,T}(v::MVector{S,T}, val::T, i::Integer) +@propagate_inbounds setindex!{S,T}(v::MVector{S,T}, val, i::Int) = setindex!(v, convert(T, val), i) +@inline function setindex!{S,T}(v::MVector{S,T}, val::T, i::Int) @boundscheck if i < 1 || i > length(v) throw(BoundsError()) end diff --git a/src/SArray.jl b/src/SArray.jl index 89bf989b..7ba59368 100644 --- a/src/SArray.jl +++ b/src/SArray.jl @@ -65,10 +65,10 @@ end @inline SArray(a::StaticArray) = SArray{size(typeof(a))}(Tuple(a)) # Some more advanced constructor-like functions -@inline eye{Size}(::Type{SArray{Size}}) = eye(SArray{Size,Float64}) -@inline zeros{Size}(::Type{SArray{Size}}) = zeros(SArray{Size,Float64}) -@inline ones{Size}(::Type{SArray{Size}}) = ones(SArray{Size,Float64}) - +@inline one(::Type{SArray{S}}) where {S} = one(SArray{S,Float64,length(S)}) +@inline eye(::Type{SArray{S}}) where {S} = eye(SArray{S,Float64,length(S)}) +@inline one(::Type{SArray{S,T}}) where {S,T} = one(SArray{S,T,length(S)}) +@inline eye(::Type{SArray{S,T}}) where {S,T} = eye(SArray{S,T,length(S)}) #################### ## SArray methods ## @@ -79,7 +79,7 @@ end @pure Size{S,T,N}(::Type{SArray{S,T,N}}) = Size(S) @pure Size{S,T,N,L}(::Type{SArray{S,T,N,L}}) = Size(S) -function getindex(v::SArray, i::Integer) +function getindex(v::SArray, i::Int) Base.@_inline_meta v.data[i] end diff --git a/src/SMatrix.jl b/src/SMatrix.jl index 9a630cfd..3c64d74f 100644 --- a/src/SMatrix.jl +++ b/src/SMatrix.jl @@ -66,7 +66,7 @@ end SMatrix{S1, S2, $T, L}(x) end end -@compat SMatrixNoType{S1, S2, L, T} = SMatrix{S1, S2, T, L} +SMatrixNoType{S1, S2, L, T} = SMatrix{S1, S2, T, L} @generated function (::Type{SMatrixNoType{S1, S2, L}}){S1,S2,L}(x::NTuple{L,Any}) T = promote_tuple_eltype(x) return quote @@ -108,9 +108,9 @@ end # Some more advanced constructor-like functions @inline one{N}(::Type{SMatrix{N}}) = one(SMatrix{N,N}) @inline eye{N}(::Type{SMatrix{N}}) = eye(SMatrix{N,N}) -@inline eye{N,M}(::Type{SMatrix{N,M}}) = eye(SMatrix{N,M,Float64}) -@inline zeros{N,M}(::Type{SMatrix{N,M}}) = zeros(SMatrix{N,M,Float64}) -@inline ones{N,M}(::Type{SMatrix{N,M}}) = ones(SMatrix{N,M,Float64}) +#@inline eye{N,M}(::Type{SMatrix{N,M}}) = eye(SMatrix{N,M,Float64}) +#@inline zeros{N,M}(::Type{SMatrix{N,M}}) = zeros(SMatrix{N,M,Float64}) +#@inline ones{N,M}(::Type{SMatrix{N,M}}) = ones(SMatrix{N,M,Float64}) ##################### ## SMatrix methods ## @@ -120,7 +120,7 @@ end @pure Size{S1,S2,T}(::Type{SMatrix{S1,S2,T}}) = Size(S1, S2) @pure Size{S1,S2,T,L}(::Type{SMatrix{S1,S2,T,L}}) = Size(S1, S2) -function getindex(v::SMatrix, i::Integer) +function getindex(v::SMatrix, i::Int) Base.@_inline_meta v.data[i] end diff --git a/src/SUnitRange.jl b/src/SUnitRange.jl new file mode 100644 index 00000000..42fb7c26 --- /dev/null +++ b/src/SUnitRange.jl @@ -0,0 +1,12 @@ +# `Start` and `End` should be `Int` +struct SUnitRange{Start,End} <: StaticVector{Int} +end + +@pure Size(::SUnitRange{Start, End}) where {Start,End} = Size(End-Start+1) + +@pure @propagate_inbounds function getindex(x::SUnitRange{Start,End}, i::Int) where {Start, End} + @boundscheck if i < Start || i > End + throw(BoundsError(x, i)) + end + return i +end diff --git a/src/SVector.jl b/src/SVector.jl index 51882a10..a83e12e9 100644 --- a/src/SVector.jl +++ b/src/SVector.jl @@ -43,7 +43,7 @@ end @pure Size{S}(::Type{SVector{S}}) = Size(S) @pure Size{S,T}(::Type{SVector{S,T}}) = Size(S) -@propagate_inbounds function getindex(v::SVector, i::Integer) +@propagate_inbounds function getindex(v::SVector, i::Int) v.data[i] end diff --git a/src/Scalar.jl b/src/Scalar.jl index 9a02b051..6a5a0a9c 100644 --- a/src/Scalar.jl +++ b/src/Scalar.jl @@ -6,9 +6,15 @@ Construct a statically-sized 0-dimensional array that contains a single element, """ immutable Scalar{T} <: StaticArray{T,0} data::T + + Scalar{T}(x::AbstractArray) where {T} = new{T}(convert(T,x)) + Scalar{T}(x::Tuple{T2}) where {T, T2} = new{T}(convert(T,x[1])) + Scalar{T}(x) where {T} = new{T}(convert(T, x)) end -@inline (::Type{Scalar{T}}){T}(x::Tuple{T}) = Scalar{T}(x[1]) +@inline Scalar(x::Tuple{T}) where {T} = Scalar{T}(x[1]) +@inline Scalar(a::AbstractArray) = Scalar{typeof(a)}(a) +@inline Scalar(a::AbstractScalar) = Scalar{eltype(a)}(a[]) # Do we want this to convert or wrap? @pure Size(::Type{Scalar}) = Size() @pure Size{T}(::Type{Scalar{T}}) = Size() @@ -24,4 +30,4 @@ end @inline Tuple(v::Scalar) = (v.data,) # A lot more compact than the default array show -Base.show{T}(io::IO, ::MIME"text/plain", x::Scalar{T}) = print(io, "Scalar{$T}(", x.data, ")") +Base.show(io::IO, ::MIME"text/plain", x::Scalar{T}) where {T} = print(io, "Scalar{$T}(", x.data, ")") diff --git a/src/SizedArray.jl b/src/SizedArray.jl index 3c5a9d76..60817d38 100644 --- a/src/SizedArray.jl +++ b/src/SizedArray.jl @@ -71,14 +71,14 @@ end @propagate_inbounds getindex(a::SizedArray, i::Int) = getindex(a.data, i) @propagate_inbounds setindex!(a::SizedArray, v, i::Int) = setindex!(a.data, v, i) -@compat SizedVector{S,T,M} = SizedArray{S,T,1,M} +SizedVector{S,T,M} = SizedArray{S,T,1,M} @pure Size{S}(::Type{SizedVector{S}}) = Size(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 -@compat SizedMatrix{S,T,M} = SizedArray{S,T,2,M} +SizedMatrix{S,T,M} = SizedArray{S,T,2,M} @pure Size{S}(::Type{SizedMatrix{S}}) = Size(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) diff --git a/src/StaticArrays.jl b/src/StaticArrays.jl index 4f17b40b..5cb7acf3 100644 --- a/src/StaticArrays.jl +++ b/src/StaticArrays.jl @@ -5,19 +5,18 @@ module StaticArrays import Base: @_inline_meta, @_propagate_inbounds_meta, @_pure_meta, @propagate_inbounds, @pure import Base: getindex, setindex!, size, similar, - length, convert, promote_op, map, map!, reduce, reducedim, + length, convert, promote_op, map, map!, reduce, reducedim, mapreducedim, mapreduce, broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill, fill!, det, inv, eig, eigvals, trace, vecnorm, norm, dot, diagm, sum, diff, prod, count, any, all, sumabs, sumabs2, minimum, - maximum, extrema, mean, copy - -using Compat + maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!, + randexp!, normalize, normalize! export StaticScalar, StaticArray, StaticVector, StaticMatrix export Scalar, SArray, SVector, SMatrix export MArray, MVector, MMatrix -export FieldVector, MutableFieldVector +export FieldVector export SizedArray, SizedVector, SizedMatrix export Size, Length @@ -28,11 +27,48 @@ export @MVector, @MMatrix, @MArray export similar_type export push, pop, shift, unshift, insert, deleteat, setindex -include("util.jl") -include("generated.jl") +""" + abstract type StaticArray{T, N} <: AbstractArray{T, N} end + StaticScalar{T} = StaticArray{T, 0} + StaticVector{T} = StaticArray{T, 1} + StaticMatrix{T} = StaticArray{T, 2} + +`StaticArray`s are Julia arrays with fixed, known size. + +## Dev docs + +They must define the following methods: + - Constructors that accept a flat tuple of data. + - `Size()` on the *type*, returning an *instance* of `Size{(dim1, dim2, ...)}` (preferably `@pure`). + - `getindex()` with an integer (linear indexing) (preferably `@inline` with `@boundscheck`). + - `Tuple()`, returning the data in a flat Tuple. + +It may be useful to implement: + +- `similar_type(::Type{MyStaticArray}, ::Type{NewElType}, ::Size{NewSize})`, returning a + type (or type constructor) that accepts a flat tuple of data. + +For mutable containers you may also need to define the following: -include("core.jl") + - `setindex!` for a single element (linear indexing). + - `similar(::Type{MyStaticArray}, ::Type{NewElType}, ::Size{NewSize})`. + - In some cases, a zero-parameter constructor, `MyStaticArray{...}()` for unintialized data + is assumed to exist. + +(see also `SVector`, `SMatrix`, `SArray`, `MVector`, `MMatrix`, `MArray`, `SizedArray` and `FieldVector`) +""" +abstract type StaticArray{T, N} <: AbstractArray{T, N} end +const StaticScalar{T} = StaticArray{T, 0} +const StaticVector{T} = StaticArray{T, 1} +const StaticMatrix{T} = StaticArray{T, 2} + +const AbstractScalar{T} = AbstractArray{T, 0} # not exported, but useful none-the-less + +include("util.jl") include("traits.jl") +include("convert.jl") + +include("SUnitRange.jl") include("Scalar.jl") include("SVector.jl") include("FieldVector.jl") @@ -43,25 +79,22 @@ include("MMatrix.jl") include("MArray.jl") include("SizedArray.jl") - include("abstractarray.jl") -include("indexing.jl") #= +include("indexing.jl") +include("broadcast.jl") include("mapreduce.jl") include("arraymath.jl") include("linalg.jl") include("matrix_multiply.jl") -include("solve.jl") -include("deque.jl") -include("det.jl") -include("inv.jl") -include("eigen.jl") -include("cholesky.jl") -=# - - -if VERSION < v"0.6.0-dev.1671" - include("FixedSizeArrays.jl") -end +#include("solve.jl") +#include("deque.jl") +#include("det.jl") +#include("inv.jl") +#include("eigen.jl") +#include("cholesky.jl") + + +#include("FixedSizeArrays.jl") # Currently defunct include("ImmutableArrays.jl") diff --git a/src/abstractarray.jl b/src/abstractarray.jl index c4b4b062..2e618549 100644 --- a/src/abstractarray.jl +++ b/src/abstractarray.jl @@ -1,14 +1,11 @@ length(a::T) where {T <: StaticArray} = prod(Size(T)) length(a::Type{T}) where {T<:StaticArray} = prod(Size(T)) -length_val(a::T) where {T <: StaticArray} = length_val(Size(T)) -length_val(a::Type{T}) where {T<:StaticArray} = length_val(Size(T)) - size{T<:StaticArray}(::T) = get(Size(T)) size{T<:StaticArray}(::Type{T}) = get(Size(T)) -size{T<:StaticArray}(::T, d::Integer) = Size(T)[d] -size{T<:StaticArray}(::Type{T}, d::Integer) = Size(T)[d] +size{T<:StaticArray}(::T, d::Int) = Size(T)[d] +size{T<:StaticArray}(::Type{T}, d::Int) = Size(T)[d] # This seems to confuse Julia a bit in certain circumstances (specifically for trailing 1's) @@ -58,25 +55,25 @@ similar_type{A<:AbstractArray,T,S}(::A,::Type{T},s::Size{S}) = similar_type(A,T, similar_type{A<:AbstractArray,T,S}(::Type{A},::Type{T},s::Size{S}) = default_similar_type(T,s,length_val(s)) default_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{0}}) = Scalar{T} -@generated default_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{1}}) = SVector{S[1],T} -@generated default_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{2}}) = SMatrix{S[1],S[2],T,prod(S)} -@generated default_similar_type{T,S,D}(::Type{T}, s::Size{S}, ::Type{Val{D}}) = SArray{S,T,D,prod(S)} +default_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{1}}) = SVector{S[1],T} +default_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{2}}) = SMatrix{S[1],S[2],T,prod(s)} +default_similar_type{T,S,D}(::Type{T}, s::Size{S}, ::Type{Val{D}}) = SArray{S,T,D,prod(s)} -# mutable things stay mutable -similar_type{SA<:Union{MVector,MMatrix,MArray},T,S}(::Type{SA},::Type{T},s::Size{S}) = mutable_similar_type(T,s,length_val(s)) +# should mutable things stay mutable? +#similar_type{SA<:Union{MVector,MMatrix,MArray},T,S}(::Type{SA},::Type{T},s::Size{S}) = mutable_similar_type(T,s,length_val(s)) mutable_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{0}}) = SizedArray{(),T,0,0} -@generated mutable_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{1}}) = MVector{S[1],T} -@generated mutable_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{2}}) = MMatrix{S[1],S[2],T,prod(S)} -@generated mutable_similar_type{T,S,D}(::Type{T}, s::Size{S}, ::Type{Val{D}}) = MArray{S,T,D,prod(S)} +mutable_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{1}}) = MVector{S[1],T} +mutable_similar_type{T,S}(::Type{T}, s::Size{S}, ::Type{Val{2}}) = MMatrix{S[1],S[2],T,prod(s)} +mutable_similar_type{T,S,D}(::Type{T}, s::Size{S}, ::Type{Val{D}}) = MArray{S,T,D,prod(s)} -# `SizedArray` stays the same, and also takes over an `Array`. -similar_type{SA<:SizedArray,T,S}(::Type{SA},::Type{T},s::Size{S}) = sizedarray_similar_type(T,s,length_val(s)) -similar_type{A<:Array,T,S}(::Type{A},::Type{T},s::Size{S}) = sizedarray_similar_type(T,s,length_val(s)) +# Should `SizedArray` stay the same, and also take over an `Array`? +#similar_type{SA<:SizedArray,T,S}(::Type{SA},::Type{T},s::Size{S}) = sizedarray_similar_type(T,s,length_val(s)) +#similar_type{A<:Array,T,S}(::Type{A},::Type{T},s::Size{S}) = sizedarray_similar_type(T,s,length_val(s)) sizedarray_similar_type{T,S,D}(::Type{T},s::Size{S},::Type{Val{D}}) = SizedArray{S,T,D,D} -# Field vectors are user controlled, and default to SVector, etc +# Field vectors are user controlled, and currently default to SVector, etc """ similar(static_array) @@ -107,103 +104,65 @@ similar{SA<:SizedArray,T,S}(::Type{SA},::Type{T},s::Size{S}) = sizedarray_simila similar{A<:Array,T,S}(::Type{A},::Type{T},s::Size{S}) = sizedarray_similar_type(T,s,length_val(s))() -# 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 -# TODO much more work on type stability. Base functions are using similar() with -# size, which poses difficulties!! -@generated function Base.full{T,SM<:StaticMatrix}(sym::Symmetric{T,SM}) - exprs_up = [i <= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(m[$(sub2ind(size(SM), j, i))]) for i = 1:size(SM,1), j=1:size(SM,2)] - exprs_lo = [i >= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(m[$(sub2ind(size(SM), j, i))]) for i = 1:size(SM,1), j=1:size(SM,2)] - - return quote - $(Expr(:meta, :inline)) - m = sym.data - if sym.uplo == 'U' - @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_up...)))) - else - @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_lo...)))) - end +@inline reshape(a::StaticArray, s::Size) = similar_type(a, s)(Tuple(a)) +@generated function reshape(a::AbstractArray, s::Size{S}) where {S} + if IndexStyle(a) == IndexLinear() # TODO this isn't "hyperpure" + exprs = [:(a[$i]) for i = 1:prod(S)] + else + exprs = [:(a[$(inds...)]) for inds ∈ CartesianRange(S)] end -end - -@generated function Base.full{T,SM<:StaticMatrix}(sym::Hermitian{T,SM}) - exprs_up = [i <= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(conj(m[$(sub2ind(size(SM), j, i))])) for i = 1:size(SM,1), j=1:size(SM,2)] - exprs_lo = [i >= j ? :(m[$(sub2ind(size(SM), i, j))]) : :(conj(m[$(sub2ind(size(SM), j, i))])) for i = 1:size(SM,1), j=1:size(SM,2)] return quote - $(Expr(:meta, :inline)) - m = sym.data - if sym.uplo == 'U' - @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_up...)))) - else - @inbounds return SM($(Expr(:call, SM, Expr(:tuple, exprs_lo...)))) + @_inline_meta + if length(a) != prod(s) + throw(DimensionMismatch("Tried to resize dynamic object of size $(size(a)) to $s")) end + return similar_type(a, s)(tuple($(exprs...))) end end -# Reshape used types to specify size (and also conveniently, the output type) -@generated function reshape{SA<:StaticArray}(a::StaticArray, ::Type{SA}) - if !(SA <: StaticVector) && length(a) != length(SA) - error("Static array of size $(size(a)) cannot be reshaped to size $(size(SA))") - end - - Base.depwarn("Use reshape(array, Size(dims...)) rather than reshape(array, StaticArrayType)", :reshape) - - return quote - $(Expr(:meta, :inline)) - return SA(Tuple(a)) - end -end - -function reshape{SA<:StaticArray}(a::AbstractArray, ::Type{SA}) - if !(SA <: StaticVector) && length(a) != length(SA) - error("Static array of size $(size(a)) cannot be reshaped to size $(size(SA))") - end +reshape(a::Array, s::Size{S}) where {S} = s(a) - Base.depwarn("Use reshape(array, Size(dims...)) rather than reshape(array, StaticArrayType)", :reshape) +@inline copy(a::StaticArray) = similar_type(a)(Tuple(a)) - return SA((a...)) -end +# TODO permutedims? +# TODO perhaps could move `Symmetric`, etc into seperate files. -# Versions using Size{} -@generated function reshape{S}(a::StaticArray, ::Size{S}) - if length(a) != prod(S) - error("Static array of size $(size(a)) cannot be reshaped to size $S") - 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 +# TODO much more work on type stability. Base functions are using similar() with +# size, which poses difficulties!! +@inline Base.full(sym::Symmetric{T,SM}) where {T,SM <: StaticMatrix} = _full(Size(SM), sym) - newtype = similar_type(a, Size(S)) +@generated function _full(::Size{S}, sym::Symmetric{T,SM}) where {S, T, SM <: StaticMatrix} + exprs_up = [i <= j ? :(m[$(sub2ind(S, i, j))]) : :(m[$(sub2ind(S, j, i))]) for i = 1:S[1], j=1:S[2]] + exprs_lo = [i >= j ? :(m[$(sub2ind(S, i, j))]) : :(m[$(sub2ind(S, j, i))]) for i = 1:S[1], j=1:S[2]] return quote - $(Expr(:meta, :inline)) - return $newtype(a) + @_inline_meta + m = sym.data + if sym.uplo == 'U' + @inbounds return SM(tuple($(exprs_up...))) + else + @inbounds return SM(tuple($(exprs_lo...))) + end end end -@generated function reshape{S}(a::Array, ::Size{S}) - newtype = SizedArray{S, eltype(a), length(S)} +@inline Base.full(herm::Hermitian{T,SM}) where {T,SM <: StaticMatrix} = _full(Size(SM), herm) - return quote - $(Expr(:meta, :inline)) - return $newtype(a) - end -end +@generated function _full(::Size{S}, herm::Hermitian{T,SM}) where {S, T, SM <: StaticMatrix} + exprs_up = [i <= j ? :(m[$(sub2ind(S, i, j))]) : :(conj(m[$(sub2ind(S, j, i))])) for i = 1:S[1], j=1:S[2]] + exprs_lo = [i >= j ? :(m[$(sub2ind(S, i, j))]) : :(conj(m[$(sub2ind(S, j, i))])) for i = 1:S[1], j=1:S[2]] -# Clever ruse to determine if a type is "mutable" -# Definitely *not* a deep copy. -@generated function copy(a::StaticArray) - try - out = a() - return quote - $(Expr(:meta, :inline)) - out = $(a)() - out .= a - return out - end - catch - return quote - $(Expr(:meta, :inline)) - a + return quote + @_inline_meta + m = herm.data + if herm.uplo == 'U' + @inbounds return SM(tuple($(exprs_up...))) + else + @inbounds return SM(tuple($(exprs_lo...))) end end end diff --git a/src/arraymath.jl b/src/arraymath.jl index c1d4e12f..93d28ec9 100644 --- a/src/arraymath.jl +++ b/src/arraymath.jl @@ -1,232 +1,149 @@ -import Base: .+, .-, .*, ./, .//, .\, .%, .^, .<<, .>>, .==, .<, .>, .<=, .>=, !, &, |, $ - # Support for elementwise ops on AbstractArray{S<:StaticArray} with Number -Base.promote_op{Op,A<:StaticArray,T<:Number}(op::Op, ::Type{A}, ::Type{T}) = similar_type(A, promote_op(op, eltype(A), T)) -Base.promote_op{Op,T<:Number,A<:StaticArray}(op::Op, ::Type{T}, ::Type{A}) = similar_type(A, promote_op(op, T, eltype(A))) - -@static if VERSION < v"0.6.0-dev.1632" - @inline .-(a1::StaticArray) = broadcast(-, a1) - - @inline .+(a1::StaticArray, a2::StaticArray) = broadcast(+, a1, a2) - @inline .-(a1::StaticArray, a2::StaticArray) = broadcast(-, a1, a2) - @inline .*(a1::StaticArray, a2::StaticArray) = broadcast(*, a1, a2) - @inline ./(a1::StaticArray, a2::StaticArray) = broadcast(/, a1, a2) - @inline .//(a1::StaticArray, a2::StaticArray) = broadcast(//, a1, a2) - @inline .\(a1::StaticArray, a2::StaticArray) = broadcast(\, a1, a2) - @inline .%(a1::StaticArray, a2::StaticArray) = broadcast(%, a1, a2) - @inline .^(a1::StaticArray, a2::StaticArray) = broadcast(^, a1, a2) - @inline .<<(a1::StaticArray, a2::StaticArray) = broadcast(<<, a1, a2) - @inline .>>(a1::StaticArray, a2::StaticArray) = broadcast(>>, a1, a2) - @inline .==(a1::StaticArray, a2::StaticArray) = broadcast(==, a1, a2) - @inline .<(a1::StaticArray, a2::StaticArray) = broadcast(<, a1, a2) - @inline .>(a1::StaticArray, a2::StaticArray) = broadcast(>, a1, a2) - @inline .<=(a1::StaticArray, a2::StaticArray) = broadcast(<=, a1, a2) - @inline .>=(a1::StaticArray, a2::StaticArray) = broadcast(>=, a1, a2) - - @inline .+(a1::StaticArray, a2::Number) = broadcast(+, a1, a2) - @inline .-(a1::StaticArray, a2::Number) = broadcast(-, a1, a2) - @inline .*(a1::StaticArray, a2::Number) = broadcast(*, a1, a2) - @inline ./(a1::StaticArray, a2::Number) = broadcast(/, a1, a2) - @inline .//(a1::StaticArray, a2::Number) = broadcast(//, a1, a2) - @inline .\(a1::StaticArray, a2::Number) = broadcast(\, a1, a2) - @inline .%(a1::StaticArray, a2::Number) = broadcast(%, a1, a2) - @inline .^(a1::StaticArray, a2::Number) = broadcast(^, a1, a2) - @inline .<<(a1::StaticArray, a2::Number) = broadcast(<<, a1, a2) - @inline .>>(a1::StaticArray, a2::Number) = broadcast(>>, a1, a2) - @inline .==(a1::StaticArray, a2) = broadcast(==, a1, a2) - @inline .<(a1::StaticArray, a2) = broadcast(<, a1, a2) - @inline .>(a1::StaticArray, a2) = broadcast(>, a1, a2) - @inline .<=(a1::StaticArray, a2) = broadcast(<=, a1, a2) - @inline .>=(a1::StaticArray, a2) = broadcast(>=, a1, a2) - - @inline .+(a1::Number, a2::StaticArray) = broadcast(+, a1, a2) - @inline .-(a1::Number, a2::StaticArray) = broadcast(-, a1, a2) - @inline .*(a1::Number, a2::StaticArray) = broadcast(*, a1, a2) - @inline ./(a1::Number, a2::StaticArray) = broadcast(/, a1, a2) - @inline .//(a1::Number, a2::StaticArray) = broadcast(//, a1, a2) - @inline .\(a1::Number, a2::StaticArray) = broadcast(\, a1, a2) - @inline .%(a1::Number, a2::StaticArray) = broadcast(%, a1, a2) - @inline .^(a1::Number, a2::StaticArray) = broadcast(^, a1, a2) - @inline .<<(a1::Number, a2::StaticArray) = broadcast(<<, a1, a2) - @inline .>>(a1::Number, a2::StaticArray) = broadcast(>>, a1, a2) - @inline .==(a1, a2::StaticArray) = broadcast(==, a1, a2) - @inline .<(a1, a2::StaticArray) = broadcast(<, a1, a2) - @inline .>(a1, a2::StaticArray) = broadcast(>, a1, a2) - @inline .<=(a1, a2::StaticArray) = broadcast(<=, a1, a2) - @inline .>=(a1, a2::StaticArray) = broadcast(>=, a1, a2) -end - -# The remaining auto-vectorized boolean operators -# TODO Julia v0.6 changes this -@inline !(a::StaticArray{Bool}) = broadcast(!, a) - -@inline (&){T1,T2}(a1::StaticArray{T1}, a2::StaticArray{T2}) = broadcast(&, a1, a2) -@inline (|){T1,T2}(a1::StaticArray{T1}, a2::StaticArray{T2}) = broadcast(|, a1, a2) -@inline ($){T1,T2}(a1::StaticArray{T1}, a2::StaticArray{T2}) = broadcast($, a1, a2) - -@inline (&){T}(a1::StaticArray{T}, a2::Number) = broadcast(&, a1, a2) -@inline (|){T}(a1::StaticArray{T}, a2::Number) = broadcast(|, a1, a2) -@inline ($){T}(a1::StaticArray{T}, a2::Number) = broadcast($, a1, a2) +#Base.promote_op{Op,A<:StaticArray,T<:Number}(op::Op, ::Type{A}, ::Type{T}) = similar_type(A, promote_op(op, eltype(A), T)) +#Base.promote_op{Op,T<:Number,A<:StaticArray}(op::Op, ::Type{T}, ::Type{A}) = similar_type(A, promote_op(op, T, eltype(A))) -@inline (&){T}(a1::Number, a2::StaticArray{T}) = broadcast(&, a1, a2) -@inline (|){T}(a1::Number, a2::StaticArray{T}) = broadcast(|, a1, a2) -@inline ($){T}(a1::Number, a2::StaticArray{T}) = broadcast($, a1, a2) - -@inline zeros{SA <: StaticArray}(::SA) = zeros(SA) -@generated function Base.zeros{SA <: StaticArray}(::Type{SA}) - s = size(SA) +@inline zeros(::SA) where {SA <: StaticArray} = zeros(SA) +@generated function zeros(::Type{SA}) where {SA <: StaticArray} T = eltype(SA) - if T == Any - T = Float64 - end - v = zeros(T, s...) - return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + if T === Any + return quote + @_inline_meta + _fill(zero(Float64), Size(SA), SA) + end + else + return quote + @_inline_meta + _fill(zero($T), Size(SA), SA) + end end end -@inline Base.ones{SA <: StaticArray}(::SA) = ones(SA) -@generated function Base.ones{SA <: StaticArray}(::Type{SA}) - s = size(SA) +@inline ones(::SA) where {SA <: StaticArray} = ones(SA) +@generated function ones(::Type{SA}) where {SA <: StaticArray} T = eltype(SA) - if T == Any - T = Float64 - end - v = ones(T, s...) - return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + if T === Any + return quote + @_inline_meta + _fill(one(Float64), Size(SA), SA) + end + else + return quote + @_inline_meta + _fill(one($T), Size(SA), SA) + end end end -@inline Base.fill{SA <: StaticArray}(val, ::SA) = fill(val, SA) -@generated function Base.fill{SA <: StaticArray}(val, ::Type{SA}) - l = length(SA) - T = eltype(SA) - expr = [:valT for i = 1:l] +@inline fill(val, ::SA) where {SA <: StaticArray} = ones(val, SA) +@inline fill(val, ::Type{SA}) where {SA <: StaticArray} = _fill(val, Size(SA), SA) +@generated function _fill(val, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} + v = [:val for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - valT = convert($T, val) - $(Expr(:call, SA, Expr(:tuple, expr...))) + @_inline_meta + $SA(tuple($(v...))) end end # Also consider randcycle, randperm? Also faster rand!(staticarray, collection) -@generated function Base.rand{SA <: StaticArray}(rng::AbstractRNG, ::Type{SA}) - s = size(SA) + +@inline rand(rng::AbstractRNG, ::SA) where {SA <: StaticArray} = rand(rng, SA) +@inline rand(rng::AbstractRNG, ::Type{SA}) where {SA <: StaticArray} = _rand(rng, Size(SA), SA) +@generated function _rand(rng::AbstractRNG, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} T = eltype(SA) if T == Any T = Float64 end v = [:(rand(rng, $T)) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + @_inline_meta + $SA(tuple($(v...))) end end -@generated function Base.rand{SA <: StaticArray}(rng::AbstractRNG, range::AbstractArray, ::Type{SA}) - s = size(SA) - T = eltype(range) +@inline rand(rng::AbstractRNG, range::AbstractArray, ::SA) where {SA <: StaticArray} = rand(rng, range, SA) +@inline rand(rng::AbstractRNG, range::AbstractArray, ::Type{SA}) where {SA <: StaticArray} = _rand(rng, range, Size(SA), SA) +@generated function _rand(rng::AbstractRNG, range::AbstractArray, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} v = [:(rand(rng, range)) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) - end -end - -@generated function Base.rand{SA <: StaticArray}(range::AbstractArray, ::Type{SA}) - s = size(SA) - T = eltype(range) - v = [:(rand(range)) for i = 1:prod(s)] - return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + @_inline_meta + $SA(tuple($(v...))) end end -@generated function Base.randn{SA <: StaticArray}(rng::AbstractRNG, ::Type{SA}) - s = size(SA) +@inline randn(rng::AbstractRNG, ::SA) where {SA <: StaticArray} = randn(rng, SA) +@inline randn(rng::AbstractRNG, ::Type{SA}) where {SA <: StaticArray} = _randn(rng, Size(SA), SA) +@generated function _randn(rng::AbstractRNG, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} T = eltype(SA) if T == Any T = Float64 end v = [:(randn(rng, $T)) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + @_inline_meta + $SA(tuple($(v...))) end end -@generated function Base.randexp{SA <: StaticArray}(rng::AbstractRNG, ::Type{SA}) - s = size(SA) +@inline randexp(rng::AbstractRNG, ::SA) where {SA <: StaticArray} = randexp(rng, SA) +@inline randexp(rng::AbstractRNG, ::Type{SA}) where {SA <: StaticArray} = _randexp(rng, Size(SA), SA) +@generated function _randexp(rng::AbstractRNG, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} T = eltype(SA) if T == Any T = Float64 end v = [:(randexp(rng, $T)) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:call, SA, Expr(:tuple, v...))) + @_inline_meta + $SA(tuple($(v...))) end end +# Mutable versions + # Why don't these two exist in Base? # @generated function Base.zeros!{SA <: StaticArray}(a::SA) # @generated function Base.ones!{SA <: StaticArray}(a::SA) -@generated function Base.fill!{SA <: StaticArray}(a::SA, val) - l = length(SA) - T = eltype(SA) - exprs = [:(@inbounds a[$i] = valT) for i = 1:l] +@inline fill!(a::SA, val) where {SA <: StaticArray} = _fill!(Size(SA), a, val) +@generated function _fill!(::Size{s}, a::SA, val) where {s, SA <: StaticArray} + exprs = [:(a[$i] = valT) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - valT = convert($T, val) - $(Expr(:block, exprs...)) + @_inline_meta + valT = convert(eltype(SA), val) + @inbounds $(Expr(:block, exprs...)) return a end end -@generated function Base.rand!{SA <: StaticArray}(rng::AbstractRNG, a::SA) - l = length(SA) - T = eltype(SA) - exprs = [:(@inbounds a[$i] = rand(rng, $T)) for i = 1:l] +@inline rand!(rng::AbstractRNG, a::SA) where {SA <: StaticArray} = _rand!(rng, Size(SA), a) +@generated function _rand!(rng::AbstractRNG, ::Size{s}, a::SA) where {s, SA <: StaticArray} + exprs = [:(a[$i] = rand(rng, eltype(SA))) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:block, exprs...)) + @_inline_meta + @inbounds $(Expr(:block, exprs...)) return a end end -@generated function Base.rand!{N}(rng::MersenneTwister, a::StaticArray{Float64,N}) # ambiguity with AbstractRNG and non-Float64... possibly an optimized form in Base? - l = length(a) - exprs = [:(@inbounds a[$i] = rand(rng, Float64)) for i = 1:l] - return quote - $(Expr(:meta, :inline)) - $(Expr(:block, exprs...)) - return a - end -end +# ambiguity with AbstractRNG and non-Float64... possibly an optimized form in Base? +@inline rand!(rng::MersenneTwister, a::SA) where {SA <: StaticArray{Float64}} = _rand!(rng, Size(SA), a) -@generated function Base.randn!{SA <: StaticArray}(rng::AbstractRNG, a::SA) - l = length(SA) - T = eltype(SA) - exprs = [:(@inbounds a[$i] = randn(rng, $T)) for i = 1:l] +@inline randn!(rng::AbstractRNG, a::SA) where {SA <: StaticArray} = _randn!(rng, Size(SA), a) +@generated function _randn!(rng::AbstractRNG, ::Size{s}, a::SA) where {s, SA <: StaticArray} + exprs = [:(a[$i] = randn(rng, eltype(SA))) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:block, exprs...)) + @_inline_meta + @inbounds $(Expr(:block, exprs...)) return a end end -@generated function Base.randexp!{SA <: StaticArray}(rng::AbstractRNG, a::SA) - l = length(SA) - T = eltype(SA) - exprs = [:(@inbounds a[$i] = randexp(rng, $T)) for i = 1:l] +@inline randexp!(rng::AbstractRNG, a::SA) where {SA <: StaticArray} = _randexp!(rng, Size(SA), a) +@generated function _randexp!(rng::AbstractRNG, ::Size{s}, a::SA) where {s, SA <: StaticArray} + exprs = [:(a[$i] = randexp(rng, eltype(SA))) for i = 1:prod(s)] return quote - $(Expr(:meta, :inline)) - $(Expr(:block, exprs...)) + @_inline_meta + @inbounds $(Expr(:block, exprs...)) return a end end diff --git a/src/broadcast.jl b/src/broadcast.jl new file mode 100644 index 00000000..22696818 --- /dev/null +++ b/src/broadcast.jl @@ -0,0 +1,147 @@ +################ +## broadcast! ## +################ + +# TODO: bad codegen for `broadcast(-, SVector(1,2,3))` + +@propagate_inbounds function broadcast(f, a::Union{Number, StaticArray}...) + _broadcast(f, broadcast_sizes(a...), a...) +end + +@inline broadcast_sizes(a...) = _broadcast_sizes((), a...) +@inline _broadcast_sizes(t::Tuple) = t +@inline _broadcast_sizes(t::Tuple, a::StaticArray, as...) = _broadcast_sizes((t..., Size(a)), as...) +@inline _broadcast_sizes(t::Tuple, a::Number, as...) = _broadcast_sizes((t..., Size()), as...) + +function broadcasted_index(oldsize, newindex) + index = ones(Int, length(oldsize)) + for i = 1:length(oldsize) + if oldsize[i] != 1 + index[i] = newindex[i] + end + end + return sub2ind(oldsize, index...) +end + +@generated function _broadcast(f, s::Tuple{Vararg{Size}}, a::Union{Number, StaticArray}...) + first_staticarray = 0 + for i = 1:length(a) + if a[i] <: StaticArray + first_staticarray = a[i] + break + end + end + + sizes = [sz.parameters[1] for sz ∈ s.parameters] + + ndims = 0 + for i = 1:length(sizes) + ndims = max(ndims, length(sizes[i])) + end + + newsize = ones(Int, ndims) + for i = 1:length(sizes) + s = sizes[i] + for j = 1:length(s) + if newsize[j] == 1 || newsize[j] == s[j] + newsize[j] = s[j] + else + throw(DimensionMismatch("Tried to broadcast on inputs sized $sizes")) + end + end + end + newsize = tuple(newsize...) + + exprs = Array{Expr}(newsize) + more = true + current_ind = ones(Int, length(newsize)) + + while more + exprs_vals = [(a[i] <: Number ? :(a[$i]) : :(a[$i][$(broadcasted_index(sizes[i], current_ind))])) for i = 1:length(sizes)] + exprs[current_ind...] = :(f($(exprs_vals...))) + + # increment current_ind + current_ind[1] += 1 + for i ∈ 1:length(newsize) + if current_ind[i] > newsize[i] + if i == length(newsize) + more = false + break + else + current_ind[i] = 1 + current_ind[i+1] += 1 + end + else + break + end + end + end + + eltype_exprs = [:(eltype($t)) for t ∈ a] + newtype_expr = :(Core.Inference.return_type(f, Tuple{$(eltype_exprs...)})) + + return quote + @_inline_meta + @inbounds return similar_type($first_staticarray, $newtype_expr, Size($newsize))(tuple($(exprs...))) + end +end + + +################ +## broadcast! ## +################ + +@propagate_inbounds function broadcast!(f, dest::StaticArray, a::Union{Number, StaticArray}...) + _broadcast!(f, Size(dest), dest, broadcast_sizes(a...), a...) +end + + +@generated function _broadcast!(f, ::Size{newsize}, dest::StaticArray, s::Tuple{Vararg{Size}}, a::Union{Number, StaticArray}...) where {newsize} + sizes = [sz.parameters[1] for sz ∈ s.parameters] + sizes = tuple(sizes...) + + ndims = 0 + for i = 1:length(sizes) + ndims = max(ndims, length(sizes[i])) + end + + for i = 1:length(sizes) + s = sizes[i] + for j = 1:length(s) + if s[j] != 1 && s[j] != (j <= length(newsize) ? newsize[j] : 1) + throw(DimensionMismatch("Tried to broadcast to destination sized $newsize from inputs sized $sizes")) + end + end + end + + exprs = Array{Expr}(newsize) + j = 1 + more = true + current_ind = ones(Int, max(length(newsize), length.(sizes)...)) + while more + exprs_vals = [(a[i] <: Number ? :(a[$i]) : :(a[$i][$(broadcasted_index(sizes[i], current_ind))])) for i = 1:length(sizes)] + exprs[current_ind...] = :(dest[$j] = f($(exprs_vals...))) + + # increment current_ind + current_ind[1] += 1 + for i ∈ 1:length(newsize) + if current_ind[i] > newsize[i] + if i == length(newsize) + more = false + break + else + current_ind[i] = 1 + current_ind[i+1] += 1 + end + else + break + end + end + j += 1 + end + + return quote + @_inline_meta + @inbounds $(Expr(:block, exprs...)) + end +end diff --git a/src/core.jl b/src/convert.jl similarity index 55% rename from src/core.jl rename to src/convert.jl index 5bd371a0..d799aa6c 100644 --- a/src/core.jl +++ b/src/convert.jl @@ -1,42 +1,7 @@ -""" - abstract type StaticArray{T, N} <: AbstractArray{T, N} end - StaticScalar{T} = StaticArray{T, 0} - StaticVector{T} = StaticArray{T, 1} - StaticMatrix{T} = StaticArray{T, 2} - -`StaticArray`s are Julia arrays with fixed, known size. - -## Dev docs - -They must define the following methods: - - Constructors that accept a flat tuple of data. - - `Size()` on the *type*, returning an *instance* of `Size{(dim1, dim2, ...)}` (preferably `@pure`). - - `getindex()` with an integer (linear indexing) (preferably `@inline` with `@boundscheck`). - - `Tuple()`, returning the data in a flat Tuple. - -It may be useful to implement: - -- `similar_type(::Type{MyStaticArray}, ::Type{NewElType}, ::Size{NewSize})`, returning a - type (or type constructor) that accepts a flat tuple of data. - -For mutable containers you may also need to define the following: - - - `setindex!` for a single element (linear indexing). - - `similar(::Type{MyStaticArray}, ::Type{NewElType}, ::Size{NewSize})`. - - In some cases, a zero-parameter constructor, `MyStaticArray{...}()` for unintialized data - is assumed to exist. - -(see also `SVector`, `SMatrix`, `SArray`, `MVector`, `MMatrix`, `MArray`, `SizedArray` and `FieldVector`) -""" -abstract type StaticArray{T, N} <: AbstractArray{T, N} end - -StaticScalar{T} = StaticArray{T, 1} -StaticVector{T} = StaticArray{T, 1} -StaticMatrix{T} = StaticArray{T, 2} - (::Type{SA})(x::Tuple) where {SA <: StaticArray} = error("No precise constructor for $SA found. Length of input was $(length(x)).") -@inline convert(::Type{SA}, x...) where {SA <: StaticArray} = SA(x) +@inline (::Type{SA})(x...) where {SA <: StaticArray} = SA(x) +@inline (::Type{SA})(a::AbstractArray) where {SA <: StaticArray} = convert(SA, a) # Is this a good idea? # this covers most conversions and "statically-sized reshapes" @inline convert(::Type{SA}, sa::StaticArray) where {SA<:StaticArray} = SA(Tuple(sa)) @@ -44,7 +9,7 @@ StaticMatrix{T} = StaticArray{T, 2} # A general way of going back to a tuple @inline function convert(::Type{Tuple}, a::StaticArray) - unroll_tuple((i -> @inbounds return a[i]), length_val(a)) + unroll_tuple(a, Length(a)) end @inline function convert(::Type{SA}, a::AbstractArray) where {SA <: StaticArray} @@ -52,9 +17,19 @@ end error("Dimension mismatch. Expected input array of length $(length(SA)), got length $(length(a))") end - return SA(unroll_tuple((i -> @inbounds return a[i]), length_val(SA))) + return SA(unroll_tuple(a, Length(SA))) end +length_val(a::T) where {T <: StaticArray} = length_val(Size(T)) +length_val(a::Type{T}) where {T<:StaticArray} = length_val(Size(T)) + +@generated function unroll_tuple(a::AbstractArray, ::Length{L}) where {L} + exprs = [:(a[$j]) for j = 1:L] + quote + @_inline_meta + @inbounds return $(Expr(:tuple, exprs...)) + end +end #= @generated function convert(::Type{Tuple}, a::StaticArray) diff --git a/src/generated.jl b/src/generated.jl deleted file mode 100644 index 36f65ab6..00000000 --- a/src/generated.jl +++ /dev/null @@ -1,15 +0,0 @@ -#= -""" - unroll_tuple(f, Val{N}) - -Like the `Base.ntuple` function, but generates fast code beyond N = 16. -""" -=# -@generated function unroll_tuple(f, ::Type{Val{N}}) where {N} - @assert(N isa Int) - exprs = [:(f($i)) for i = 1:N] - quote - @_propagate_inbounds_meta - return $(Expr(:tuple, exprs...)) - end -end diff --git a/src/indexing.jl b/src/indexing.jl index 43c51918..8f2cf667 100644 --- a/src/indexing.jl +++ b/src/indexing.jl @@ -1,823 +1,312 @@ +# Default error messages to help users with new types and to avoid subsequent stack overflows +getindex(a::StaticArray, i::Int) = error("getindex(::$typeof(a), ::Int) is not defined.") +setindex!(a::StaticArray, value, i::Int) = error("setindex!(::$(typeof(a)), value, ::Int) is not defined.") + ###################### -## Linear Indexing ## +## Scalar Indexing ## ###################### -# Static Vector indexing into AbstractArrays -@propagate_inbounds function getindex(a::AbstractArray{T}, inds::StaticVector{<:Integer}) where {T} - similar_type(a, T, Size(inds))(unroll_tuple(Length(inds)) do i - @_propagate_inbounds_meta - a[inds[i]] - end) -end - -#= -@generated function getindex{T, I <: Integer}( - a::AbstractArray{T}, inds::StaticVector{I} - ) - S = length(inds) - newtype = similar_type(inds, T, Size(S)) - exprs = [:(a[inds[$i]]) for i = 1:S] - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - +# Note: all indexing behavior defaults to dense, linear indexing -# Convert to StaticArrays using tuples -# TODO think about bounds checks here. -@generated function getindex{T, I <: Integer}( - m::AbstractArray{T}, inds1::StaticVector{I}, i2::Integer - ) - S = length(inds1) - newtype = similar_type(inds1, T, Size(S)) # drop singular dimension like in base - exprs = [:(m[inds1[$j], i2]) for j = 1:S] - return Expr(:call, newtype, Expr(:tuple, exprs...)) +@propagate_inbounds function getindex(a::StaticArray, inds::Int...) + _getindex_scalar(Size(a), a, inds...) end -@generated function getindex{T, I <: Integer}( - m::AbstractArray{T}, i1::Integer, inds2::StaticVector{I} - ) - S = length(inds2) - newtype = similar_type(inds2, T, Size(S)) - exprs = [:(m[i1, inds2[$j]]) for j = 1:S] - return Expr(:call, newtype, Expr(:tuple, exprs...)) -end - -@generated function getindex{I1 <: Integer,I2 <: Integer, T}( - m::AbstractArray{T}, inds1::StaticVector{I1}, inds2::StaticVector{I2} - ) - S1, S2 = length(inds1), lengths(inds2) - exprs = [:(m[inds1[$j1], inds2[$j2]]) for j1 = 1:S1, j2 = 1:S2] - return Expr(:call, SMatrix{S1, S2, T}, Expr(:tuple, exprs...)) # I guess SMatrix should be fine? -end - - -@generated function getindex{SA<:StaticArray}(a::SA, ::Colon) - if SA <: StaticVector && a == similar_type(SA) - return quote - $(Expr(:meta, :inline)) - a - end - else - l = length(SA) - inds = 1:l - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - $(Expr(:call, :getindex, :a, Expr(:tuple, inds...))) +@generated function _getindex_scalar(::Size{S}, a::StaticArray, inds::Int...) where S + stride = 1 + for i ∈ 1:length(inds) + if i == 1 + ind_expr = :(inds[1]) + else + ind_expr = :($ind_expr + $stride * (inds[$i] - 1)) end + stride *= Size(S)[i] end -end - -# MAYBE: fixed-size indexing with bools? They would have to be Val's... - -# Size-indeterminate linear indexing seems to be provided by AbstractArray, -# returning a `Vector`. - -# We seem to get an error from Base's implementation with UnitRange -#= -function Base.getindex(v::StaticVector, r::UnitRange) - l = length(r) - out = similar(v, (l,)) - for i in r - out[i] = v[i] - end - return out -end -=# - - -# Same for setindex! -@generated function setindex!{SA <: StaticArray, I <: Integer}( - a::SA, vals, inds::Union{Tuple{Vararg{I}}, StaticVector{I}} - ) - S = inds <: Tuple ? length(inds.parameters) : length(inds) - exprs = [:(a[inds[$i]] = vals[$i]) for i = 1:S] - if vals <: StaticArray - if length(vals) != S # We should be able to annotate that the RHS of the above exprs is @inbounds, but not the LHS... - error("Dimension mismatch") - end - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - $(Expr(:block, exprs...)) - return :vals - end - else - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - @boundscheck if length(vals) != $S - error("Dimension mismatch") - end - $(Expr(:block, exprs...)) - return :vals - end + return quote + @_propagate_inbounds_meta + a[$ind_expr] end end -# this one for consistency -@generated function setindex!{I <: Integer}( - a::Array, vals::AbstractArray, inds::StaticVector{I} - ) - S = inds <: Tuple ? length(inds.parameters) : length(inds) - exprs = [:(a[inds[$i]] = vals[$i]) for i = 1:S] - if vals <: StaticArray - if length(vals) != S - error("Dimension mismatch") - end - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - $(Expr(:block, exprs...)) - return :vals - end - else - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - @boundscheck if length(vals) != $S - error("Dimension mismatch") - end - $(Expr(:block, exprs...)) - return :vals - end - end +@propagate_inbounds function setindex!(a::StaticArray, value, inds::Int...) + _setindex!_scalar(Size(a), a, value, inds...) end -@generated function setindex!{SA<:StaticArray}(a::SA, vals, ::Colon) - exprs = [:(a[$i] = vals[$i]) for i = 1:length(SA)] - if vals <: StaticArray - if length(vals) != length(SA) - error("Dimension mismatch") - end - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - $(Expr(:block, exprs...)) - return :vals - end - else - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - @boundscheck if length(vals) != $(length(SA)) - error("Dimension mismatch") - end - $(Expr(:block, exprs...)) - return :vals +@generated function _setindex!_scalar(::Size{S}, a::StaticArray, value, inds::Int...) where S + stride = 1 + for i ∈ 1:length(inds) + if i == 1 + ind_expr = :(inds[1]) + else + ind_expr = :($ind_expr + $stride * (inds[$i] - 1)) end + stride *= S[i] end -end - -############################### -## Two-dimensional Indexing ## -############################### -# Special 2D case to begin with (and possibly good for avoiding splatting penalties?) -# Furthermore, avoids stupidity regarding two-dimensional indexing on 3+ dimensional arrays! -@generated function getindex{SM<:StaticMatrix}(m::SM, i1::Integer, i2::Integer) return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - @boundscheck if (i1 < 1 || i1 > $(size(SM,1)) || i2 < 1 || i2 > $(size(SM,2))) - throw(BoundsError(m, (i1,i2))) - end - - @inbounds return m[i1 + $(size(SM,1))*(i2-1)] + @_propagate_inbounds_meta + a[$ind_expr] = value end end -# TODO put bounds checks here, as they should have less overhead here -@generated function getindex{SM<:StaticMatrix, S1, S2}(m::SM, inds1::NTuple{S1,Integer}, inds2::NTuple{S2,Integer}) - newtype = similar_type(SM, Size(S1, S2)) - exprs = [:(m[inds1[$i1], inds2[$i2]]) for i1 = 1:S1, i2 = 1:S2] - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function getindex{SM<:StaticMatrix, S2}(m::SM, i1::Integer, inds2::NTuple{S2,Integer}) - newtype = similar_type(SM, Size(S2)) - exprs = [:(m[i1, inds2[$i2]]) for i2 = 1:S2] - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end +######################### +## Indexing utilities ## +######################### -@generated function getindex{SM<:StaticMatrix, S1}(m::SM, inds1::NTuple{S1,Integer}, i2::Integer) - newtype = similar_type(SM, Size(S1)) - exprs = [:(m[inds1[$i1], i2]) for i1 = 1:S1] +@pure increment(::Type{Val{N}}) where {N} = Val{N+1} - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end +@inline index_sizes(s::Size, inds...) = _index_sizes(s, Val{1}, (), inds...) +@inline _index_sizes(s::Size, ::Type{Val{N}}, x::Tuple) where {N} = x +@inline _index_sizes(s::Size, v::Type{Val{N}}, x::Tuple, ::Int, inds...) where {N} = _index_sizes(s, increment(v), (x..., Size()), inds...) +@inline _index_sizes(s::Size, v::Type{Val{N}}, x::Tuple, a::StaticArray, inds...) where {N} = _index_sizes(s, increment(v), (x..., Size(a)), inds...) +@inline _index_sizes(s::Size, v::Type{Val{N}}, x::Tuple, a::Colon, inds...) where {N} = _index_sizes(s, increment(v), (x..., Size(s[N])), inds...) +@inline index_sizes(inds...) = _index_sizes(Val{1}, (), inds...) +@inline _index_sizes(::Type{Val{N}}, x::Tuple) where {N} = x +@inline _index_sizes(v::Type{Val{N}}, x::Tuple, ::Int, inds...) where {N} = _index_sizes(increment(v), (x..., Size()), inds...) +@inline _index_sizes(v::Type{Val{N}}, x::Tuple, a::StaticArray, inds...) where {N} = _index_sizes(increment(v), (x..., Size(a)), inds...) +out_index_size(ind_sizes::Type{<:Size}...) = Size(_out_index_size((), ind_sizes...)) +@inline _out_index_size(t::Tuple) = t +@inline _out_index_size(t::Tuple, ::Type{Size{S}}, ind_sizes...) where {S} = _out_index_size((t..., S...), ind_sizes...) +linear_index_size(ind_sizes::Type{<:Size}...) = _linear_index_size((), ind_sizes...) +@inline _linear_index_size(t::Tuple) = t +@inline _linear_index_size(t::Tuple, ::Type{Size{S}}, ind_sizes...) where {S} = _linear_index_size((t..., prod(S)), ind_sizes...) +_ind(i::Int, ::Int, ::Type{Int}) = :(inds[$i]) +_ind(i::Int, j::Int, ::Type{<:StaticArray}) = :(inds[$i][$j]) +_ind(i::Int, j::Int, ::Type{Colon}) = j -# TODO: bounds checking? -@generated function setindex!{SM<:StaticMatrix, S1, S2}(m::SM, val, inds1::NTuple{S1,Integer}, inds2::NTuple{S2,Integer}) - exprs = [:(m[inds1[$i1], inds2[$i2]] = val[$i1,$i2]) for i1 = 1:S1, i2 = 1:S2] +##################### +## Array Indexing ## +##################### - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - $(Expr(:block, exprs...)) - return val - end -end -# TODO put bounds checks here, as they should have less overhead here -@generated function getindex{SM<:StaticMatrix}(m::SM, ::Colon, inds2::Union{Integer, Tuple{Vararg{Integer}}}) - inds1 = ntuple(identity, size(SM,1)) - quote - $(Expr(:meta, :inline, :propagate_inbounds)) - m[$inds1, inds2] - end +@propagate_inbounds function getindex(a::StaticArray, inds::Union{Int, StaticArray{Int}, Colon}...) + _getindex(a, index_sizes(Size(a), inds...), inds) end -# TODO put bounds checks here, as they should have less overhead here -@generated function getindex{SM<:StaticMatrix}(m::SM, inds1::Union{Integer, Tuple{Vararg{Integer}}}, ::Colon) - inds2 = ntuple(identity, size(SM,2)) - quote - $(Expr(:meta, :inline, :propagate_inbounds)) - m[inds1, $inds2] - end +# Hard to describe "Union{Int, StaticArray{Int}} with at least one StaticArray{Int}" +# Here we require the first StaticArray{Int} to be within the first four dimensions +@propagate_inbounds function getindex(a::AbstractArray, i1::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, inds...), (i1, inds...)) end - -@generated function getindex{SM<:StaticMatrix}(m::SM, ::Colon, ::Colon) - inds1 = ntuple(identity, size(SM,1)) - inds2 = ntuple(identity, size(SM,2)) - quote - $(Expr(:meta, :inline)) - @inbounds return m[$inds1, $inds2] - end +@propagate_inbounds function getindex(a::AbstractArray, i1::Int, i2::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, inds...), (i1, i2, inds...)) end - -@generated function setindex!{SM<:StaticMatrix, S2}(m::SM, val, i1::Integer, inds2::NTuple{S2,Integer}) - newtype = similar_type(SM, Size(S2)) - exprs = [:(m[i1, inds2[$i2]] = val[$i2]) for i2 = 1:S2] - - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:block, exprs...)) - end +@propagate_inbounds function getindex(a::AbstractArray, i1::Int, i2::Int, i3::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, i3, inds...), (i1, i2, i3, inds...)) end -# TODO expand out tuples to vectors in size-indeterminate cases -# 2D setindex! - -@generated function setindex!{SM<:StaticMatrix}(m::SM, val, i1::Integer, i2::Integer) - return quote - $(Expr(:meta, :inline)) - @boundscheck if (i1 < 1 || i1 > $(size(SM,1)) || i2 < 1 || i2 > $(size(SM,2))) - throw(BoundsError(m, (i1,i2))) - end - - @inbounds return m[i1 + $(size(SM,1))*(i2-1)] = val - end +@propagate_inbounds function getindex(a::AbstractArray, i1::Int, i2::Int, i3::Int, i4::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, i3, i4, inds...), (i1, i2, i3, i4, inds...)) end - - -@generated function setindex!{SM<:StaticMatrix, S1}(m::SM, val, inds1::NTuple{S1,Integer}, i2::Integer) - newtype = similar_type(SM, Size(S1)) - exprs = [:(m[inds1[$i1], i2] = val[$i1]) for i1 = 1:S1] - - return quote - $(Expr(:meta, :inline, :propagate_inbounds)) - return $(Expr(:block, exprs...)) - end +# Disambuguity methods for the above +@propagate_inbounds function getindex(a::StaticArray, i1::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, inds...), (i1, inds...)) end -# TODO put bounds checks here, as they should have less overhead here -@generated function setindex!{SM<:StaticMatrix}(m::SM, val, ::Colon, inds2::Union{Integer, Tuple{Vararg{Integer}}}) - inds1 = ntuple(identity, size(SM,1)) - quote - $(Expr(:meta, :inline, :propagate_inbounds)) - m[$inds1, inds2] = val - end +@propagate_inbounds function getindex(a::StaticArray, i1::Int, i2::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, inds...), (i1, i2, inds...)) end -# TODO put bounds checks here, as they should have less overhead here -@generated function setindex!{SM<:StaticMatrix}(m::SM, val, inds1::Union{Integer, Tuple{Vararg{Integer}}}, ::Colon) - inds2 = ntuple(identity, size(SM,2)) - quote - $(Expr(:meta, :inline, :propagate_inbounds)) - m[inds1, $inds2] = val - end +@propagate_inbounds function getindex(a::StaticArray, i1::Int, i2::Int, i3::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, i3, inds...), (i1, i2, i3, inds...)) end -# Check bounds for val here? -@generated function setindex!{SM<:StaticMatrix}(m::SM, val, ::Colon, ::Colon) - inds1 = ntuple(identity, size(SM,1)) - inds2 = ntuple(identity, size(SM,2)) - quote - $(Expr(:meta, :inline)) - return m[$inds1, $inds2] = val - end +@propagate_inbounds function getindex(a::StaticArray, i1::Int, i2::Int, i3::Int, i4::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _getindex(a, index_sizes(i1, i2, i3, i4, inds...), (i1, i2, i3, i4, inds...)) end -################################# -## Multi-dimensional Indexing ## -################################# - -# Scalar indexing: -@generated function getindex(a::StaticArray, i1::Integer, i2::Integer, i3::Integer) - @assert ndims(a) <= 3 +@generated function _getindex(a::AbstractArray, ind_sizes::Tuple{Vararg{Size}}, inds) + newsize = out_index_size(ind_sizes.parameters...) + linearsizes = linear_index_size(ind_sizes.parameters...) + exprs = Array{Expr}(linearsizes) - strides = [1, size(a,1), size(a,1)*size(a,2)] + # Iterate over input indices + ind_types = inds.parameters + current_ind = ones(Int,length(linearsizes)) + more = true + while more + exprs_tmp = [_ind(i, current_ind[i], ind_types[i]) for i = 1:length(linearsizes)] + exprs[current_ind...] = :(getindex(a, $(exprs_tmp...))) - N = 3 - ind_exprs = [:i1, :i2, :i3] - exprs = [j == 1 ? ind_exprs[1] : :( $(strides[j]) * ($(ind_exprs[j])-1) ) for j = 1:N] - ind_expr = Expr(:call, :+, exprs...) - - bounds_expr = :(i1 < 1 || i1 > $(size(a,1)) || i2 < 1 || i2 > $(size(a,2)) || i3 < 1 || i3 > $(size(a,3))) - - return quote - $(Expr(:meta, :inline)) - @boundscheck if $bounds_expr - throw(BoundsError(a, (i1,i2,i3))) - end - - @inbounds return a[$ind_expr] - end -end - -# TODO check speed (splatting penalty) -@generated function getindex(a::StaticArray, i1::Integer, i2::Integer, i_n::Integer...) - @assert ndims(a) <= 2 + length(i_n) - - N = 2 + length(i_n) - - strides = ones(Int, N) - for j = 2:N - for k = 1:j - 1 - strides[j] *= size(a, k) + # increment current_ind + current_ind[1] += 1 + for i ∈ 1:length(linearsizes) + if current_ind[i] > linearsizes[i] + if i == length(linearsizes) + more = false + break + else + current_ind[i] = 1 + current_ind[i+1] += 1 + end + else + break + end end end - ind_exprs = [:i1, :i2, [:(i_n[$j]) for j = 1:length(i_n)]...] - exprs = [j == 1 ? ind_exprs[1] : :( $(strides[j]) * ($(ind_exprs[j])-1) ) for j = 1:N] - ind_expr = Expr(:call, :+, exprs...) - - bounds_expr = :(i1 < 1 || i1 > $(size(a,1)) || i2 < 1 || i2 > $(size(a,2))) - for j = 1:N-2 - bound_expr = :($bounds_expr || i_n[$j] < 1 || i_n[$j] > $(size(a,2+j))) - end - - return quote - $(Expr(:meta, :inline)) - @boundscheck if $bounds_expr - throw(BoundsError(a, (i1,i2,i_n...))) - end - - @inbounds return a[$ind_expr] + quote + @_propagate_inbounds_meta + similar_type(a, $newsize)(tuple($(exprs...))) end end -@generated function setindex!(a::StaticArray, val, i1::Integer, i2::Integer, i3::Integer) - @assert ndims(a) <= 3 - - strides = [1, size(a,1), size(a,1)*size(a,2)] - N = 3 - ind_exprs = [:i1, :i2, :i3] - exprs = [j == 1 ? ind_exprs[1] : :($(strides[j]) * ($(ind_exprs[j])-1) ) for j = 1:N] - ind_expr = Expr(:call, :+, exprs...) - bounds_expr = :(i1 < 1 || i1 > $(size(a,1)) || i2 < 1 || i2 > $(size(a,2)) || i3 < 1 || i3 > $(size(a,3))) - - return quote - $(Expr(:meta, :inline)) - @boundscheck if $bounds_expr - throw(BoundsError(a, (i1,i2,3))) - end - - @inbounds return a[$ind_expr] = val - end +@propagate_inbounds function setindex!(a::StaticArray, value, inds::Union{Int, StaticArray{Int}, Colon}...) + _setindex!(a, value, index_sizes(Size(a), inds...), inds) end -# TODO check speed (splatting penalty) -@generated function setindex!(a::StaticArray, val, i1::Integer, i2::Integer, i_n::Integer...) - @assert ndims(a) <= 2 + length(i_n) - - N = 2 + length(i_n) - - strides = ones(Int, N) - for j = 2:N - for k = 1:j - 1 - strides[j] *= size(a, k) - end - end - - ind_exprs = [:i1, :i2, [:(i_n[$j]) for j = 1:length(i_n)]...] - exprs = [j == 1 ? ind_exprs[1] : :( $(strides[j]) * ($(ind_exprs[j])-1) ) for j = 1:N] - ind_expr = Expr(:call, :+, exprs...) - - bounds_expr = :(i1 < 1 || i1 > $(size(a,1)) || i2 < 1 || i2 > $(size(a,2))) - for j = 1:N-2 - bound_expr = :($bounds_expr || i_n[$j] < 1 || i_n[$j] > $(size(a,2+j))) - end - - return quote - $(Expr(:meta, :inline)) - @boundscheck if $bounds_expr - throw(BoundsError(a, (i1,i2,i_n...))) - end - - @inbounds return a[$ind_expr] = val - end +# Hard to describe "Union{Int, StaticArray{Int}} with at least one StaticArray{Int}" +# Here we require the first StaticArray{Int} to be within the first four dimensions +@propagate_inbounds function setindex!(a::AbstractArray, value, i1::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, inds...), (i1, inds...)) end -=# - - -# TODO higher dimensional versions with tuples and : - -#= -############## -## Indexing ## -############## -# Indexing with no components -Base.getindex(a::StaticArray) = a.data[1] - -# Can index linearly with a scalar, a tuple, or colon (overspecified to ovoid ambiguity problems in julia 0.5) -Base.getindex{Sizes,T,N,D}(a::SArray{Sizes,T,N,D}, i::Int) = a.data[i] -Base.getindex{Sizes,T,N,D}(a::MArray{Sizes,T,N,D}, i::Int) = a.data[i] -@generated function Base.getindex{N}(a::SArray, i::NTuple{N,Int}) - newtype = similar_type(a, Val{(N,)}) - exprs = ntuple(n -> :(a[i[$n]]), N) - return :($newtype($(Expr(:tuple, exprs...)))) -end -@generated function Base.getindex{N}(a::MArray, i::NTuple{N,Int}) - newtype = similar_type(a, Val{(N,)}) - exprs = ntuple(n -> :(a[i[$n]]), N) - return :($newtype($(Expr(:tuple, exprs...)))) +@propagate_inbounds function setindex!(a::AbstractArray, value, i1::Int, i2::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, inds...), (i1, i2, inds...)) end -Base.getindex(a::SArray, ::Colon) = SVector(a.data) -Base.getindex(a::MArray, ::Colon) = MVector(a.data) # Makes a copy? - -# Multidimensional index generalizes the above -@generated function Base.getindex{Sizes,T,N}(a::SArray{Sizes,T,N}, i...) - # Striding lengths - strides = [1, cumprod(collect(Sizes)[1:end-1])...] - - # Get the parameters of the new matrix - apl_slicing = VERSION >= v"0.5-" - NewN = 0 - NewSizes = Vector{Int}() - OldSizes = Vector{Int}() # Same as NewSizes but includes singleton 1's - at_end = true - is_singleton = trues(N) - for j = length(i):-1:1 - if i[j] == Int - unshift!(OldSizes,1) - if !apl_slicing - if !at_end - NewN += 1 - unshift!(NewSizes,1) - is_singleton[j] = false - end - end - elseif i[j] <: TupleN{Int} - NewN += 1 - unshift!(NewSizes,length(i[j].parameters)) - unshift!(OldSizes,length(i[j].parameters)) - is_singleton[j] = false - at_end = false - elseif i[j] == Colon - NewN += 1 - unshift!(NewSizes,Sizes[j]) - unshift!(OldSizes,Sizes[j]) - is_singleton[j] = false - at_end = false - else - str = "Cannot index dimension $j of $a with a $(i[j]). Use an Int, NTuple{N,Int} or Colon." - return :(error($str)) - end - end - NewSizes = (NewSizes...) - NewM = prod(NewSizes) - # Bail early if possible - if NewN == 0 - return :(a.data[sub2ind(Sizes,i...)]) - end - - NewType = MArray{NewSizes,T,NewN,NTuple{NewM,T}} - - # Now we build an expression for each new element - exprs = Vector{Expr}() - inds_old = ones(Int,N) - for j = 1:NewM - sum_exprs = ntuple(p -> :($(strides[p]) * (i[$p][$(inds_old[p])] - 1)), N) - push!(exprs, :(a.data[$(Expr(:call, :+, 1, sum_exprs...))])) - - if j < NewM - inds_old[1] += 1 - end - for k = 1:N - if inds_old[k] > OldSizes[k] - inds_old[k] = 1 - inds_old[k+1] += 1 - else - break - end - end - end - - return :(SArray{$NewSizes,$T,$NewN,NTuple{$NewM,$T}}($(Expr(:tuple, exprs...)))) +@propagate_inbounds function setindex!(a::AbstractArray, value, i1::Int, i2::Int, i3::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, inds...), (i1, i2, i3, inds...)) end -@generated function Base.getindex{Sizes,T,N}(a::MArray{Sizes,T,N}, i...) - # Striding lengths - strides = [1, cumprod(collect(Sizes)[1:end-1])...] - - # Get the parameters of the new matrix - apl_slicing = VERSION >= v"0.5-" - NewN = 0 - NewSizes = Vector{Int}() - OldSizes = Vector{Int}() # Same as NewSizes but includes singleton 1's - at_end = true - is_singleton = trues(N) - for j = length(i):-1:1 - if i[j] == Int - unshift!(OldSizes,1) - if !apl_slicing - if !at_end - NewN += 1 - unshift!(NewSizes,1) - is_singleton[j] = false - end - end - elseif i[j] <: TupleN{Int} - NewN += 1 - unshift!(NewSizes,length(i[j].parameters)) - unshift!(OldSizes,length(i[j].parameters)) - is_singleton[j] = false - at_end = false - elseif i[j] == Colon - NewN += 1 - unshift!(NewSizes,Sizes[j]) - unshift!(OldSizes,Sizes[j]) - is_singleton[j] = false - at_end = false - else - str = "Cannot index dimension $j of $a with a $(i[j]). Use an Int, NTuple{N,Int} or Colon." - return :(error($str)) - end - end - NewSizes = (NewSizes...) - NewM = prod(NewSizes) - - # Bail early if possible - if NewN == 0 - return :(a.data[sub2ind(Sizes,i...)]) - end - - NewType = MArray{NewSizes,T,NewN,NTuple{NewM,T}} - - # Now we build an expression for each new element - exprs = Vector{Expr}() - inds_old = ones(Int,N) - for j = 1:NewM - sum_exprs = ntuple(p -> :($(strides[p]) * (i[$p][$(inds_old[p])] - 1)), N) - push!(exprs, :(a.data[$(Expr(:call, :+, 1, sum_exprs...))])) - - if j < NewM - inds_old[1] += 1 - end - for k = 1:N - if inds_old[k] > OldSizes[k] - inds_old[k] = 1 - inds_old[k+1] += 1 - else - break - end - end - end - - return :(MArray{$NewSizes,$T,$NewN,NTuple{$NewM,$T}}($(Expr(:tuple, exprs...)))) +@propagate_inbounds function setindex!(a::AbstractArray, value, i1::Int, i2::Int, i3::Int, i4::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, i4, inds...), (i1, i2, i3, i4, inds...)) end -# setindex! (no index) -Base.setindex!(a::MArray, v) = setindex!(a, v, 1) -Base.unsafe_setindex!{Sizes,T}(a::MArray{Sizes,T}, v) = Base.unsafe_setindex!(a, v, 1) - - -# setindex! (linear scalar) -function Base.setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i::Int) - if i < 0 || i > prod(Sizes) - throw(BoundsError(a,i)) - end - - Base.unsafe_setindex!(a, v, i) +# Disambiguity methods for the above +@propagate_inbounds function setindex!(a::StaticArray, value, i1::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, inds...), (i1, inds...)) end -function Base.unsafe_setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i::Int) - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - Base.unsafe_store!(p,convert(T,v),i) +@propagate_inbounds function setindex!(a::StaticArray, value, i1::Int, i2::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, inds...), (i1, i2, inds...)) end - -# setindex! arbitray set of indices -function Base.setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i) - M = prod(Sizes) - for ind in i - if ind < 0 || ind > M - throw(BoundsError(a,i)) - end - end - - Base.unsafe_setindex!(a, v, i) +@propagate_inbounds function setindex!(a::StaticArray, value, i1::Int, i2::Int, i3::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, inds...), (i1, i2, i3, inds...)) end -function Base.unsafe_setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i) - N = length(i) - - # Check if v is OK in size - if length(v) != N - throw(DimensionMismatch("tried to assign $(prod(size(v))) elements to $N destinations")) - end - - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for j in i - Base.unsafe_store!(p,convert(T,v[j]),i[j]) - end +@propagate_inbounds function setindex!(a::StaticArray, value, i1::Int, i2::Int, i3::Int, i4::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, i4, inds...), (i1, i2, i3, i4, inds...)) end - -# setindex! (linear tuple optimization) -@generated function Base.unsafe_setindex!{Sizes,T,N}(a::MArray{Sizes,T}, v::StaticArray{Sizes}, i::NTuple{N,Int}) - # Compile-time check if v is OK in size - if length(v) != N - str = "tried to assign $(length(v)) elements to $N destinations" - return :(throw(DimensionMismatch($str))) - end - - quote - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for j in i - Base.unsafe_store!(p,convert(T,v[j]),i[j]) - end - end +# disambiguities from Base +@propagate_inbounds function setindex!(a::Array, value::AbstractArray, i1::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, inds...), (i1, inds...)) end -# setindex! linear with : (Colon) -function Base.setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i::Colon) - Base.unsafe_setindex!(a, v, i) +@propagate_inbounds function setindex!(a::Array, value::AbstractArray, i1::Int, i2::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, inds...), (i1, i2, inds...)) end -function Base.unsafe_setindex!{Sizes,T}(a::MArray{Sizes,T}, v, ::Colon) - # Check if v is OK in size - if length(v) != length(a) - throw(DimensionMismatch("tried to assign $(prod(size(v))) elements to $(length(a)) destinations")) - end - - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for j = 1:length(a) - Base.unsafe_store!(p,convert(T,v[j]),j) - end +@propagate_inbounds function setindex!(a::Array, value::AbstractArray, i1::Int, i2::Int, i3::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, inds...), (i1, i2, i3, inds...)) end -# setindex(m_array, static_array, :) - optimization for know sizes -@generated function Base.unsafe_setindex!{Sizes1,Sizes2,T}(a::MArray{Sizes1,T}, v::StaticArray{Sizes2}, i::Colon) - # Compile-time check if v is OK in size - if length(v) != length(a) - str = "tried to assign $(length(v)) elements to $(length(a)) destinations" - return :(throw(DimensionMismatch($str))) - end - - quote - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for j = 1:$(length(a)) - Base.unsafe_store!(p,convert(T,v[j]),j) - end - end +@propagate_inbounds function setindex!(a::Array, value::AbstractArray, i1::Int, i2::Int, i3::Int, i4::StaticArray{Int}, inds::Union{Int, StaticArray{Int}}...) + _setindex!(a, value, index_sizes(i1, i2, i3, i4, inds...), (i1, i2, i3, i4, inds...)) end +# setindex! from a scalar +@generated function _setindex!(a::AbstractArray, value, ind_sizes::Tuple{Vararg{Size}}, inds) + linearsizes = linear_index_size(ind_sizes.parameters...) + exprs = Array{Expr}(linearsizes) -# setindex! multi-dimensional general case -function Base.setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i...) - # check that all i > 0 - for ii in 1:length(i) - if isa(i[ii], Colon) - continue - elseif any(j -> j<1 || j>Sizes[ii], i[ii]) - error("Indices must be positive integers and in-range: got $i for $Sizes-dimensional array") - end - end - - Base.unsafe_setindex!(a, v, i...) -end - -function Base.unsafe_setindex!{Sizes,T}(a::MArray{Sizes,T}, v, i...) - # First check v and i have correct sizes (ignoring singleton dimensions, mimicking Base.Array) - i_sizes = ntuple(j -> i[j] == Colon() ? Sizes[j] : length(i[j]), length(i) ) + # Iterate over input indices + ind_types = inds.parameters + current_ind = ones(Int,length(ind_types)) + more = true + while more + exprs_tmp = [_ind(i, current_ind[i], ind_types[i]) for i = 1:length(ind_types)] + exprs[current_ind...] = :(setindex!(a, value, $(exprs_tmp...))) - # "Flatten" out any singleton dimensions on input v and indices i - Sizes2 = size(v) - Sizes2_flattened = Vector{Int}() - for j = 1:length(Sizes2) - if Sizes2[j] > 1 - push!(Sizes2_flattened, Sizes2[j]) + # increment current_ind + current_ind[1] += 1 + for i ∈ 1:length(linearsizes) + if current_ind[i] > linearsizes[i] + if i == length(linearsizes) + more = false + break + else + current_ind[i] = 1 + current_ind[i+1] += 1 + end + else + break + end end end - i_sizes_flattened = Vector{Int}() - for j = 1:length(i_sizes) - if i_sizes[j] > 1 - push!(i_sizes_flattened, i_sizes[j]) + quote + @_propagate_inbounds_meta + $(exprs...) + return value + end +end + + +# setindex! from an array +@generated function _setindex!(a::AbstractArray, v::AbstractArray, ind_sizes::Tuple{Vararg{Size}}, inds) + linearsizes = linear_index_size(ind_sizes.parameters...) + exprs = Array{Expr}(linearsizes) + + # Iterate over input indices + ind_types = inds.parameters + current_ind = ones(Int,length(ind_types)) + more = true + j = 1 + while more + exprs_tmp = [_ind(i, current_ind[i], ind_types[i]) for i = 1:length(ind_types)] + exprs[current_ind...] = :(setindex!(a, v[$j], $(exprs_tmp...))) + + # increment current_ind + current_ind[1] += 1 + for i ∈ 1:length(linearsizes) + if current_ind[i] > linearsizes[i] + if i == length(linearsizes) + more = false + break + else + current_ind[i] = 1 + current_ind[i+1] += 1 + end + else + break + end end + j += 1 end - # If they aren't consistent, then we can't find a sensible way to assign the data - if Sizes2_flattened != i_sizes_flattened - println(Sizes2_flattened) - println(i_sizes_flattened) - throw(DimensionMismatch("tried to assign $(Sizes2)-dimensional data to $(i_sizes)-dimensional destination")) - end - - M = prod(Sizes2) - - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for ind_v = 1:M - sub_i = ind2sub(i_sizes, ind_v) - sub_a = ntuple(k -> i[k][sub_i[k]], length(i_sizes)) - ind_a = sub2ind(Sizes, sub_a...) - Base.unsafe_store!(p, convert(T, v[ind_v]), ind_a) - end -end - -# setindex! multi-dimensional run-time-optimized case -@generated function Base.unsafe_setindex!{Sizes1,Sizes2,T}(a::MArray{Sizes1,T}, v::StaticArray{Sizes2}, i::Union{Int,TupleN{Int},Colon}...) - # First check v and i have correct sizes (ignoring singleton dimensions, mimicking Base.Array) - i_sizes = ntuple(j -> i[j] == Colon ? Sizes1[j] : length(i[j].parameters), length(i) ) - - # "Flatten" out any singleton dimensions on input v and indices i - Sizes2_flattened = Vector{Int}() - for j = 1:length(Sizes2) - if Sizes2[j] > 1 - push!(Sizes2_flattened, Sizes2[j]) + if v <: StaticArray + if Length(v) != prod(linearsizes) + return DimensionMismatch("tried to assign $(length(v))-element array to $newsize destination") end - end - - i_sizes_flattened = Vector{Int}() - for j = 1:length(i_sizes) - if i_sizes[j] > 1 - push!(i_sizes_flattened, i_sizes[j]) + quote + @_propagate_inbounds_meta + $(exprs...) + return v end - end - - # If they aren't consistent, then we can't find a sensible way to assign the data - if Sizes2_flattened != i_sizes_flattened - str = "tried to assign $(Sizes2)-dimensional data to $(i_sizes)-dimensional destination" - return :(throw(DimensionMismatch($str))) - end - - M = prod(Sizes2) - - return quote - # Get a pointer to the object - p = Base.data_pointer_from_objref(a) - p = Base.unsafe_convert(Ptr{T}, p) - - # Store the value - for ind_v = 1:$M - sub_i = ind2sub($(i_sizes), ind_v) - sub_a = ntuple(k -> i[k][sub_i[k]], $(length(i_sizes))) - ind_a = sub2ind($(Sizes1), sub_a...) - Base.unsafe_store!(p, convert(T, v[ind_v]), ind_a) + else + quote + @_propagate_inbounds_meta + if length(v) != $(prod(linearsizes)) + newsize = $linearsizes + throw(DimensionMismatch("tried to assign $(length(v))-element array to $newsize destination")) + end + $(exprs...) + return v end end end -=# diff --git a/src/linalg.jl b/src/linalg.jl index 07f9f869..f78dbcde 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -30,67 +30,20 @@ import Base: +, -, *, /, \ # With UniformScaling -@generated function +(a::StaticMatrix, b::UniformScaling) - n = size(a,1) - newtype = similar_type(a, promote_type(eltype(a), eltype(b))) +@inline +(a::StaticMatrix, b::UniformScaling) = _plus_uniform(Size(a), a, b.λ) +@inline +(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), b, a.λ) +@inline -(a::StaticMatrix, b::UniformScaling) = _plus_uniform(Size(a), a, -b.λ) +@inline -(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), -b, a.λ) - if size(a,2) != n - error("Dimension mismatch") - end - - exprs = [i == j ? :(a[$(sub2ind(size(a), i, j))] + b.λ) : :(a[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n] - - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function +(a::UniformScaling, b::StaticMatrix) - n = size(b,1) - newtype = similar_type(b, promote_type(eltype(a), eltype(b))) - - if size(b,2) != n - error("Dimension mismatch") - end - - exprs = [i == j ? :(a.λ + b[$(sub2ind(size(b), i, j))]) : :(b[$(sub2ind(size(b), i, j))]) for i = 1:n, j = 1:n] - - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function -(a::StaticMatrix, b::UniformScaling) - n = size(a,1) - newtype = similar_type(a, promote_type(eltype(a), eltype(b))) - - if size(a,2) != n - error("Dimension mismatch") +@generated function _plus_uniform(::Size{S}, a::StaticMatrix, λ) where {S} + if S[1] != S[2] + throw(DimensionMismatch("matrix is not square: dimensions are $S")) end - - exprs = [i == j ? :(a[$(sub2ind(size(a), i, j))] - b.λ) : :(a[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n] - + n = S[1] + exprs = [i == j ? :(a[$(sub2ind(size(a), i, j))] + λ) : :(a[$(sub2ind(size(a), i, j))]) for i = 1:n, j = 1:n] return quote $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function -(a::UniformScaling, b::StaticMatrix) - n = size(b,1) - newtype = similar_type(b, promote_type(eltype(a), eltype(b))) - - if size(b,2) != n - error("Dimension mismatch") - end - - exprs = [i == j ? :(a.λ - b[$(sub2ind(size(b), i, j))]) : :(-b[$(sub2ind(size(b), i, j))]) for i = 1:n, j = 1:n] - - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type(a, promote_type(eltype(a), typeof(λ)))(tuple($(exprs...))) end end @@ -104,6 +57,7 @@ end # TODO different methods for v0.5, v0.6 (due to `RowVector`) @inline conj(a::StaticArray) = map(conj, a) +#= @generated function transpose(v::StaticVector) n = length(v) newtype = similar_type(v, Size(1,n)) @@ -125,60 +79,56 @@ end @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) end end +=# -@generated function transpose(m::StaticMatrix) - (s1,s2) = size(m) - if s1 == s2 - newtype = m - else - newtype = similar_type(m, Size(s2,s1)) - end +@inline transpose(m::StaticMatrix) = _transpose(Size(m), m) + +@generated function _transpose(::Size{S}, m::StaticMatrix) where {S} + Snew = (S[2], S[1]) - exprs = [:(m[$j1, $j2]) for j2 = 1:s2, j1 = 1:s1] + exprs = [:(m[$(sub2ind(S, j1, j2))]) for j2 = 1:S[2], j1 = 1:S[1]] return quote $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type($m, Size($Snew))(tuple($(exprs...))) end end -@generated function ctranspose(m::StaticMatrix) - (s1,s2) = size(m) - if s1 == s2 - newtype = m - else - newtype = similar_type(m, Size(s2,s1)) - end +@inline ctranspose(m::StaticMatrix) = _transpose(Size(m), m) + +@generated function _ctranspose(::Size{S}, m::StaticMatrix) where {S} + Snew = (S[2], S[1]) - exprs = [:(conj(m[$j1, $j2])) for j2 = 1:s2, j1 = 1:s1] + exprs = [:(conj(m[$(sub2ind(S, j1, j2))])) for j2 = 1:S[2], j1 = 1:S[1]] return quote $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type($m, Size($Snew))(tuple($(exprs...))) end end - @inline vcat(a::Union{StaticVector,StaticMatrix}) = a -@generated function vcat(a::Union{StaticVector, StaticMatrix}, b::Union{StaticVector,StaticMatrix}) - if size(a,2) != size(b,2) - error("Dimension mismatch") +@inline vcat(a::Union{StaticVector, StaticMatrix}, b::Union{StaticVector,StaticMatrix}) = _vcat(Size(a), Size(b), a, b) +@generated function _vcat(::Size{Sa}, ::Size{Sb}, a::Union{StaticVector, StaticMatrix}, b::Union{StaticVector,StaticMatrix}) where {Sa, Sb} + if Size(Sa)[2] != Size(Sb)[2] + throw(DimensionMismatch("Tried to vcat arrays of size $Sa and $Sb")) end + # TODO cleanup? if a <: StaticVector && b <: StaticVector - newtype = similar_type(a, Size(length(a) + length(b))) - exprs = vcat([:(a[$i]) for i = 1:length(a)], - [:(b[$i]) for i = 1:length(b)]) + Snew = (Sa[1] + Sb[1],) + exprs = vcat([:(a[$i]) for i = 1:Sa[1]], + [:(b[$i]) for i = 1:Sb[1]]) else - newtype = similar_type(a, Size(size(a,1) + size(b,1), size(a,2))) + Snew = (Sa[1] + Sb[1], Size(Sa)[2]) exprs = [((i <= size(a,1)) ? ((a <: StaticVector) ? :(a[$i]) : :(a[$i,$j])) : ((b <: StaticVector) ? :(b[$(i-size(a,1))]) : :(b[$(i-size(a,1)),$j]))) - for i = 1:(size(a,1)+size(b,1)), j = 1:size(a,2)] + for i = 1:(Sa[1]+Sb[1]), j = 1:Size(Sa)[2]] end return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...))) end end # TODO make these more efficient @@ -187,28 +137,24 @@ end @inline vcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}, c::Union{StaticVector,StaticMatrix}...) = vcat(vcat(a,b), c...) -@generated function hcat(a::StaticVector) - newtype = similar_type(a, Size(length(a),1)) - exprs = [:(a[$i]) for i = 1:length(a)] - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end + +@inline hcat(a::StaticVector) = similar_type(a, Size(Size(a)[1],1))(a) @inline hcat(a::StaticMatrix) = a -@generated function hcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}) - if size(a,1) != size(b,1) +@inline hcat(a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}) = _hcat(Size(a), Size(b), a, b) + +@generated function _hcat(::Size{Sa}, ::Size{Sb}, a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}) where {Sa, Sb} + if Sa[1] != Sb[1] error("Dimension mismatch") end - exprs1 = [:(a[$i]) for i = 1:length(a)] - exprs2 = [:(b[$i]) for i = 1:length(b)] + exprs = vcat([:(a[$i]) for i = 1:prod(Sa)], + [:(b[$i]) for i = 1:prod(Sb)]) - newtype = similar_type(a, Size(size(a,1), size(a,2) + size(b,2))) + Snew = (Sa[1], Size(Sa)[2] + Size(Sb)[2]) return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs1..., exprs2...))) + @_inline_meta + @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...))) end end # TODO make these more efficient @@ -220,108 +166,96 @@ end @inline Base.zero{SA <: StaticArray}(a::SA) = zeros(SA) @inline Base.zero{SA <: StaticArray}(a::Type{SA}) = zeros(SA) -@inline one{SM <: StaticMatrix}(::SM) = one(SM) -@generated function one{SM <: StaticArray}(::Type{SM}) - s = size(SM) - if (length(s) != 2) || (s[1] != s[2]) +@inline one(::SM) where {SM <: StaticMatrix} = _one(Size(SM), SM) +@inline one(::Type{SM}) where {SM <: StaticMatrix} = _one(Size(SM), SM) +@generated function _one(::Size{S}, ::Type{SM}) where {S, SM <: StaticArray} + if (length(S) != 2) || (S[1] != S[2]) error("multiplicative identity defined only for square matrices") end - T = eltype(SM) + T = eltype(SM) # should be "hyperpure" if T == Any T = Float64 end - e = eye(T, s...) + exprs = [i == j ? :(one($T)) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]] return quote $(Expr(:meta, :inline)) - $(Expr(:call, SM, Expr(:tuple, e...))) + SM(tuple($(exprs...))) end end -@inline eye{SM <: StaticMatrix}(::SM) = eye(SM) -@generated function eye{SM <: StaticArray}(::Type{SM}) - s = size(SM) - if length(s) != 2 - error("`eye` is only defined for matrices") - end - T = eltype(SM) +@inline eye(::SM) where {SM <: StaticMatrix} = _eye(Size(SM), SM) +@inline eye(::Type{SM}) where {SM <: StaticMatrix} = _eye(Size(SM), SM) +@generated function _eye(::Size{S}, ::Type{SM}) where {S, SM <: StaticArray} + T = eltype(SM) # should be "hyperpure" if T == Any T = Float64 end - e = eye(T, s...) + exprs = [i == j ? :(one($T)) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]] return quote $(Expr(:meta, :inline)) - $(Expr(:call, SM, Expr(:tuple, e...))) + SM(tuple($(exprs...))) end end -@generated function diagm(v::StaticVector) +@inline diagm(v::StaticVector) = _diagm(Size(v), v) +@generated function _diagm(::Size{S}, v::StaticVector) where {S} + Snew = (S[1], S[1]) T = eltype(v) - exprs = [i == j ? :(v[$i]) : zero(T) for i = 1:length(v), j = 1:length(v)] - newtype = similar_type(v, Size(length(v), length(v))) + exprs = [i == j ? :(v[$i]) : zero(T) for i = 1:S[1], j = 1:S[1]] return quote $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type($v, Size($Snew))(tuple($(exprs...))) end end -@generated function cross(a::StaticVector, b::StaticVector) - S = typeof(zero(eltype(a))*zero(eltype(b))) - if length(a) === 3 && length(b) === 3 - 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 - elseif length(a) === 2 && length(b) === 2 - return quote - $(Expr(:meta, :inline)) - promote_type(eltype(a), eltype(b))(a[1]*b[2]-a[2]*b[1]) - end - elseif length(a) === 2 && length(b) === 2 - return quote - $(Expr(:meta, :inline)) - promote_type(eltype(a), eltype(b))(a[1]*b[2]-a[2]*b[1]) - end - else - error("Cross product only defined for 3-vectors") - end +@inline cross(a::StaticVector, b::StaticVector) = _cross(same_size(a, b), a, b) +_cross(::Size{S}, a::StaticVector, b::StaticVector) where {S} = error("Cross product not defined for $(S[1])-vectors") +@inline function _cross(::Size{(2,)}, a::StaticVector, b::StaticVector) + @inbounds return a[1]*b[2] - a[2]*b[1] +end +@inline function _cross(::Size{(3,)}, a::StaticVector, b::StaticVector) + T = typeof(a[2]*b[3]-a[3]*b[2]) + @inbounds return similar_type(a, typeof(a[2]*b[3]-a[3]*b[2]))((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 +@inline function _cross(::Size{(2,)}, a::StaticVector{<:Unsigned}, b::StaticVector{<:Unsigned}) + @inbounds return Signed(a[1]*b[2]) - Signed(a[2]*b[1]) +end +@inline function _cross(::Size{(3,)}, a::StaticVector{<:Unsigned}, b::StaticVector{<:Unsigned}) + T = typeof(a[2]*b[3]-a[3]*b[2]) + @inbounds return similar_type(a, typeof(Signed(a[2]*b[3])-Signed(a[3]*b[2])))(((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 -@generated function dot(a::StaticVector, b::StaticVector) - if length(a) == length(b) - expr = :(conj(a[1]) * b[1]) - for j = 2:length(a) - expr = :($expr + conj(a[$j]) * b[$j]) - end +@inline dot(a::StaticVector, b::StaticVector) = _dot(same_size(a, b), a, b) +@generated function _dot(::Size{S}, a::StaticVector, b::StaticVector) where {S} + if S[1] == 0 + return :(zero(promote_op(*, eltype(a), eltype(b)))) + end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - else - error("dot() expects vectors of same length. Got lengths $(length(a)) and $(length(b)).") + expr = :(conj(a[1]) * b[1]) + for j = 2:S[1] + expr = :($expr + conj(a[$j]) * b[$j]) + end + + return quote + @_inline_meta + @inbounds return $expr end end -@generated function vecdot(a::StaticArray, b::StaticArray) - if length(a) == length(b) - expr = :(conj(a[1]) * b[1]) - for j = 2:length(a) - expr = :($expr + conj(a[$j]) * b[$j]) - end +@inline vecdot(a::StaticArray, b::StaticArray) = _vecdot(same_size(a, b), a, b) +@generated function _vecdot(::Size{S}, a::StaticArray, b::StaticArray) where {S} + if prod(S) == 0 + return :(zero(promote_op(*, eltype(a), eltype(b)))) + end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - else - error("vecdot() expects arrays of same length. Got sizes $(size(a)) and $(size(b)).") + expr = :(conj(a[1]) * b[1]) + for j = 2:prod(S) + expr = :($expr + conj(a[$j]) * b[$j]) + end + + return quote + @_inline_meta + @inbounds return $expr end end @@ -330,13 +264,14 @@ end @inline Base.LinAlg.norm_sqr(v::StaticVector) = mapreduce(abs2, +, zero(real(eltype(v))), v) -@generated function vecnorm(a::StaticArray) - if length(a) == 0 +@inline vecnorm(a::StaticArray) = _vecnorm(Size(a), a) +@generated function _vecnorm(::Size{S}, a::StaticArray) where {S} + if prod(S) == 0 return zero(real(eltype(a))) end expr = :(real(conj(a[1]) * a[1])) - for j = 2:length(a) + for j = 2:prod(S) expr = :($expr + real(conj(a[$j]) * a[$j])) end @@ -348,18 +283,19 @@ end _norm_p0(x) = x == 0 ? zero(x) : one(x) -@generated function vecnorm(a::StaticArray, p::Real) - if length(a) == 0 +@inline vecnorm(a::StaticArray, p::Real) = _vecnorm(Size(a), a, p) +@generated function _vecnorm(::Size{S}, a::StaticArray, p::Real) where {S} + if prod(S) == 0 return zero(real(eltype(a))) end expr = :(abs(a[1])^p) - for j = 2:length(a) + for j = 2:prod(S) expr = :($expr + abs(a[$j])^p) end expr_p1 = :(abs(a[1])) - for j = 2:length(a) + for j = 2:prod(S) expr_p1 = :($expr_p1 + abs(a[$j])) end @@ -379,36 +315,36 @@ _norm_p0(x) = x == 0 ? zero(x) : one(x) end end -@inline Base.normalize(a::StaticVector) = inv(vecnorm(a))*a -@inline Base.normalize(a::StaticVector, p::Real) = inv(vecnorm(a, p))*a +@inline normalize(a::StaticVector) = inv(vecnorm(a))*a +@inline normalize(a::StaticVector, p::Real) = inv(vecnorm(a, p))*a -@inline Base.normalize!(a::StaticVector) = (a .*= inv(vecnorm(a))) -@inline Base.normalize!(a::StaticVector, p::Real) = (a.*= inv(vecnorm(a, p))) +@inline normalize!(a::StaticVector) = (a .*= inv(vecnorm(a))) +@inline normalize!(a::StaticVector, p::Real) = (a .*= inv(vecnorm(a, p))) - -@generated function trace(a::StaticMatrix) - s = size(a) - if s[1] != s[2] - error("matrix is not square") +@inline trace(a::StaticMatrix) = _trace(Size(a), a) +@generated function _trace(::Size{S}, a::StaticMatrix) where {S} + if S[1] != S[2] + throw(DimensionMismatch("matrix is not square")) end - if s[1] == 0 + if S[1] == 0 return zero(eltype(a)) end - exprs = [:(a[$(sub2ind(s, i, i))]) for i = 1:s[1]] + exprs = [:(a[$(sub2ind(S, i, i))]) for i = 1:S[1]] total = reduce((ex1, ex2) -> :($ex1 + $ex2), exprs) return quote - $(Expr(:meta, :inline)) + @_inline_meta @inbounds return $total end end -@inline Size{T,SA<:StaticArray}(::Union{Symmetric{T,SA}, Type{Symmetric{T,SA}}}) = Size(SA) -@inline Size{T,SA<:StaticArray}(::Union{Hermitian{T,SA}, Type{Hermitian{T,SA}}}) = Size(SA) +# TODO same for `RowVector`? +@inline Size(::Union{Symmetric{T,SA}, Type{Symmetric{T,SA}}}) where {T,SA<:StaticArray} = Size(SA) +@inline Size(::Union{Hermitian{T,SA}, Type{Hermitian{T,SA}}}) where {T,SA<:StaticArray} = Size(SA) -# some micro-optimizations +# some micro-optimizations (TODO check these make sense for v0.6) @inline Base.LinAlg.checksquare{SM<:StaticMatrix}(::SM) = _checksquare(Size(SM)) @inline Base.LinAlg.checksquare{SM<:StaticMatrix}(::Type{SM}) = _checksquare(Size(SM)) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index d0ea5e80..67544efa 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -1,157 +1,60 @@ -######### -## map ## -######### - -# Single input -@generated function map{T}(f, a1::StaticArray{T}) - exprs = [:(f(a1[$j])) for j = 1:length(a1)] - return quote - $(Expr(:meta, :inline)) - newtype = similar_type(typeof(a1), promote_op(f, T)) - @inbounds return $(Expr(:call, :newtype, Expr(:tuple, exprs...))) +# Returns the common Size of the inputs (or else throws a DimensionMismatch) +@inline same_size(a1::StaticArray, as::StaticArray...) = _same_size(Size(a1), as...) +@inline _same_size(s::Size) = s +@inline function _same_size(s::Size, a1::StaticArray, as::StaticArray...) + if s === Size(a1) + return _same_size(s, as...) + else + throw(DimensionMismatch("Dimensions must match. Got inputs with $s and $(Size(a1)).")) end end -# Two inputs -@generated function map{T1,T2}(f, a1::StaticArray{T1}, a2::StaticArray{T2}) - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))") - end +@inline _first(a1, as...) = a1 - exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a1)] - return quote - $(Expr(:meta, :inline)) - newtype = similar_type(typeof(a1), promote_op(f, T1, T2)) - @inbounds return $(Expr(:call, :newtype, Expr(:tuple, exprs...))) - end -end - -# TODO these assume linear fast... -@generated function map{T1,T2}(f, a1::StaticArray{T1}, a2::AbstractArray{T2}) - exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a1)] - return quote - $(Expr(:meta, :inline)) - - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))") - end +################ +## map / map! ## +################ - newtype = similar_type(typeof(a1), promote_op(f, T1, T2)) - @inbounds return $(Expr(:call, :newtype, Expr(:tuple, exprs...))) - end +@inline function map(f, a::StaticArray...) + _map(f, same_size(a...), a...) end -@generated function map{T1,T2}(f, a1::AbstractArray{T1}, a2::StaticArray{T2}) - exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a2)] +@generated function _map(f, ::Size{S}, a::StaticArray...) where {S} + exprs = Vector{Expr}(prod(S)) + for i ∈ 1:prod(S) + tmp = [:(a[$j][$i]) for j ∈ 1:length(a)] + exprs[i] = :(f($(tmp...))) + end + eltypes = [eltype(a[j]) for j ∈ 1:length(a)] # presumably, `eltype` is "hyperpure"? + newT = :(Core.Inference.return_type(f, Tuple{$(eltypes...)})) return quote - $(Expr(:meta, :inline)) - - @boundscheck if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))") - end - - newtype = similar_type(typeof(a2), promote_op(f, T1, T2)) - @inbounds return $(Expr(:call, :newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(typeof(_first(a...)), $newT)(tuple($(exprs...))) end end -# TODO General case involving arbitrary many inputs? - -############ -## reduce ## -############ -@generated function reduce(op, a::StaticArray) - if length(a) == 1 - return :(@inbounds return a[1]) - else - expr = :(op(a[1], a[2])) - for j = 3:length(a) - expr = :(op($expr, a[$j])) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - end +@inline function map!(f, dest::StaticArray, a::StaticArray...) + _map!(f, dest, same_size(dest, a...), a...) end -@generated function reduce(op, v0, a::StaticArray) - if length(a) == 0 - return :(v0) - else - expr = :(op(v0, a[1])) - for j = 2:length(a) - expr = :(op($expr, a[$j])) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - end +# Ambiguities with Base: +@inline function map!(f, dest::StaticArray, a::StaticArray) + _map!(f, dest, same_size(dest, a), a) end - -@generated function reducedim{D}(op, a::StaticArray, ::Type{Val{D}}) - S = size(a) - if S[D] == 1 - return :(return a) - else - N = ndims(a) - Snew = ([n==D ? 1 : S[n] for n = 1:N]...) - newtype = similar_type(a, Size(Snew)) - - exprs = Array{Expr}(Snew) - itr = [1:n for n = Snew] - for i = Base.product(itr...) - ik = copy([i...]) - ik[D] = 2 - expr = :(op(a[$(i...)], a[$(ik...)])) - for k = 3:S[D] - ik[D] = k - expr = :(op($expr, a[$(ik...)])) - end - - exprs[i...] = expr - end - - return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end - end +@inline function map!(f, dest::StaticArray, a::StaticArray, b::StaticArray) + _map!(f, dest, same_size(dest, a, b), a, b) end -# These are all similar in Base but not @inline'd -@inline sum{T}(a::StaticArray{T}) = reduce(+, zero(T), a) -@inline prod{T}(a::StaticArray{T}) = reduce(*, one(T), a) -@inline count(a::StaticArray{Bool}) = reduce(+, 0, a) -@inline all(a::StaticArray{Bool}) = reduce(&, true, a) # non-branching versions -@inline any(a::StaticArray{Bool}) = reduce(|, false, a) # (benchmarking needed) -@inline mean(a::StaticArray) = sum(a) / length(a) -@inline sumabs{T}(a::StaticArray{T}) = mapreduce(abs, +, zero(T), a) -@inline sumabs2{T}(a::StaticArray{T}) = mapreduce(abs2, +, zero(T), a) -@inline minimum(a::StaticArray) = reduce(min, a) # base has mapreduce(idenity, scalarmin, a) -@inline maximum(a::StaticArray) = reduce(max, a) # base has mapreduce(idenity, scalarmax, a) -@inline minimum{D}(a::StaticArray, dim::Type{Val{D}}) = reducedim(min, a, dim) -@inline maximum{D}(a::StaticArray, dim::Type{Val{D}}) = reducedim(max, a, dim) - -@generated function diff{D}(a::StaticArray, ::Type{Val{D}}=Val{1}) - S = size(a) - N = ndims(a) - Snew = ([n==D ? S[n]-1 : S[n] for n = 1:N]...) - newtype = similar_type(a, Size(Snew)) - exprs = Array{Expr}(Snew) - itr = [1:n for n = Snew] - - for i1 = Base.product(itr...) - i2 = copy([i1...]) - i2[D] = i1[D] + 1 - exprs[i1...] = :(a[$(i2...)] - a[$(i1...)]) +@generated function _map!(f, dest, ::Size{S}, a::StaticArray...) where {S} + exprs = Vector{Expr}(prod(S)) + for i ∈ 1:prod(S) + tmp = [:(a[$j][$i]) for j ∈ 1:length(a)] + exprs[i] = :(dest[$i] = f($(tmp...))) end - return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds $(Expr(:block, exprs...)) end end @@ -159,272 +62,153 @@ end ## mapreduce ## ############### -# Single array -@generated function mapreduce(f, op, a1::StaticArray) - if length(a1) == 1 - return :(f(a1[1])) - else - expr = :(op(f(a1[1]), f(a1[2]))) - for j = 3:length(a1) - expr = :(op($expr, f(a1[$j]))) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - end +@inline function mapreduce(f, op, a::StaticArray...) + _mapreduce(f, op, same_size(a...), a...) end -@generated function mapreduce(f, op, v0, a1::StaticArray) - if length(a1) == 0 - return :(v0) - else - expr = :(op(v0, f(a1[1]))) - for j = 2:length(a1) - expr = :(op($expr, f(a1[$j]))) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end - end +@inline function mapreduce(f, op, v0, a::StaticArray...) + _mapreduce(f, op, v0, same_size(a...), a...) end -# Two arrays (e.g. dot has f(a,b) = a' * b, op = +) -@generated function mapreduce(f, op, a1::StaticArray, a2::StaticArray) - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))") +@generated function _mapreduce(f, op, ::Size{S}, a::StaticArray...) where {S} + tmp = [:(a[$j][1]) for j ∈ 1:length(a)] + expr = :(f($(tmp...))) + for i ∈ 2:prod(S) + tmp = [:(a[$j][$i]) for j ∈ 1:length(a)] + expr = :(op($expr, f($(tmp...)))) end - - if length(a1) == 1 - return :(f(a1[1], a2[1])) - else - expr = :(op(f(a1[1], a2[1]), f(a1[2], a2[2]))) - for j = 3:length(a1) - expr = :(op($expr, f(a1[$j], a2[$j]))) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end + return quote + @_inline_meta + @inbounds return $expr end end -@generated function mapreduce(f, op, v0, a1::StaticArray, a2::StaticArray) - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))") +@generated function _mapreduce(f, op, v0, ::Size{S}, a::StaticArray...) where {S} + expr = :v0 + for i ∈ 1:prod(S) + tmp = [:(a[$j][$i]) for j ∈ 1:length(a)] + expr = :(op($expr, f($(tmp...)))) end - - if length(a1) == 0 - return :(v0) - else - expr = :(op(v0, f(a1[1], a2[1]))) - for j = 2:length(a1) - expr = :(op($expr, f(a1[$j], a2[$j]))) - end - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end + return quote + @_inline_meta + @inbounds return $expr end end -# TODO General case involving arbitrary many inputs? - -############### -## broadcast ## -############### -# Single input version -@inline broadcast(f, a::StaticArray) = map(f, a) - - -# Two input versions -@generated function broadcast(f, a1::StaticArray, a2::StaticArray) - if size(a1) == size(a2) - return quote - $(Expr(:meta, :inline)) - map(f, a1, a2) - end - else - s1 = size(a1) - s2 = size(a2) - ndims = max(length(s1), length(s2)) - - s = Vector{Int}(ndims) - expands1 = Vector{Bool}(ndims) - expands2 = Vector{Bool}(ndims) - for i = 1:ndims - if length(s1) < i - s[i] = s2[i] - expands1[i] = false - expands2[i] = s2[i] > 1 - elseif length(s2) < i - s[i] = s1[i] - expands1[i] = s1[i] > 1 - expands2[i] = false - else - s[i] = max(s1[i], s1[i]) - @assert s1[i] == 1 || s1[i] == s[i] - @assert s2[i] == 1 || s2[i] == s[i] - expands1[i] = s1[i] > 1 - expands2[i] = s2[i] > 1 - end - end - s = (s...) - L = prod(s) - - if s == s1 - newtype = :( similar_type($a1, promote_op(f, $(eltype(a1)), $(eltype(a2)))) ) - else - newtype = :( similar_type($a1, promote_op(f, $(eltype(a1)), $(eltype(a2))), $(Size(s))) ) - end +################## +## mapreducedim ## +################## - exprs = Vector{Expr}(L) - - i = 1 - ind = ones(Int, ndims) - while i <= L - ind1 = [expands1[j] ? ind[j] : 1 for j = 1:length(s1)] - ind2 = [expands2[j] ? ind[j] : 1 for j = 1:length(s2)] - - exprs[i] = Expr(:call, :f, Expr(:ref, :a1, ind1...), Expr(:ref, :a2, ind2...)) +# I'm not sure why the signature for this from Base precludes multiple arrays? +# (also, why now mutating `mapreducedim!` and `reducedim!`?) +# (similarly, `broadcastreduce` and `broadcastreducedim` sounds useful) +@inline function mapreducedim(f, op, a::StaticArray, ::Type{Val{D}}) where {D} + _mapreducedim(f, op, Size(a), a, Val{D}) +end - i += 1 +@inline function mapreducedim(f, op, a::StaticArray, ::Type{Val{D}}, v0) where {D} + _mapreducedim(f, op, Size(a), a, Val{D}, v0) +end - ind[1] += 1 - j = 1 - while j < length(s) - if ind[j] > s[j] - ind[j] = 1 - ind[j+1] += 1 - else - break - end +@generated function _mapreducedim(f, op, ::Size{S}, a::StaticArray, ::Type{Val{D}}) where {S, D} + N = length(S) + Snew = ([n==D ? 1 : S[n] for n = 1:N]...) - j += 1 - end + exprs = Array{Expr}(Snew) + itr = [1:n for n ∈ Snew] + for i ∈ Base.product(itr...) + expr = :(f(a[$(i...)])) + for k = 2:S[D] + ik = collect(i) + ik[D] = k + expr = :(op($expr, f(a[$(ik...)]))) end - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end + exprs[i...] = expr end -end - -@inline broadcast(f, a::StaticArray, n::Number) = map(x -> f(x, n), a) -@inline broadcast(f, n::Number, a::StaticArray) = map(x -> f(n, x), a) -# Other two-input versions with AbstractArray - -########## -## map! ## -########## - -# Single input -@generated function map!{F}(f::F, out::StaticArray, a1::StaticArray) - exprs = [:(out[$j] = f(a1[$j])) for j = 1:length(a1)] + # TODO element type might change return quote - $(Expr(:meta, :inline)) - @inbounds $(Expr(:block, exprs...)) + @_inline_meta + @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...))) end end -# Two inputs -@generated function map!{F}(f::F, out::StaticArray, a1::StaticArray, a2::StaticArray) - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))") - end +@generated function _mapreducedim(f, op, ::Size{S}, a::StaticArray, ::Type{Val{D}}, v0) where {S, D} + N = ndims(a) + Snew = ([n==D ? 1 : S[n] for n = 1:N]...) - if size(a1) != size(a2) - error("Dimensions must match. Got sizes $(size(out)) and $(size(a1))") + exprs = Array{Expr}(Snew) + itr = [1:n for n = Snew] + for i ∈ Base.product(itr...) + expr = :v0 + for k = 1:S[D] + ik[D] = k + expr = :(op($expr, f(a[$(ik...)]))) + end + + exprs[i...] = expr end - exprs = [:(out[$j] = f(a1[$j], a2[$j])) for j = 1:length(a1)] + # TODO element type might change return quote - #$(Expr(:meta, :inline)) - @inbounds $(Expr(:block, exprs...)) + @_inline_meta + @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...))) end end +############ +## reduce ## +############ -################ -## broadcast! ## -################ - -@inline broadcast!{F}(f::F, out::StaticArray, a::StaticArray) = map!(f, out, a) -@inline broadcast!(f::typeof(identity), out::StaticArray, a::StaticArray) = map!(f, out, a) +@inline reduce(op, a::StaticArray) = mapreduce(identity, op, a) +@inline reduce(op, v0, a::StaticArray) = mapreduce(identity, op, v0, a) -# Two input versions -@generated function broadcast!{F}(f::F, out::StaticArray, a1::StaticArray, a2::StaticArray) - if size(a1) == size(a2) && size(out) == size(a1) - return quote - $(Expr(:meta, :inline)) - @inbounds map!(f, out, a1, a2) - end - else - s1 = size(a1) - s2 = size(a2) - ndims = max(length(s1), length(s2)) - - s = Vector{Int}(ndims) - expands1 = Vector{Bool}(ndims) - expands2 = Vector{Bool}(ndims) - for i = 1:ndims - if length(s1) < i - s[i] = s2[i] - expands1[i] = false - expands2[i] = s2[i] > 1 - elseif length(s2) < i - s[i] = s1[i] - expands1[i] = s1[i] > 1 - expands2[i] = false - else - s[i] = max(s1[i], s1[i]) - @assert s1[i] == 1 || s1[i] == s[i] - @assert s2[i] == 1 || s2[i] == s[i] - expands1[i] = s1[i] > 1 - expands2[i] = s2[i] > 1 - end - end - s = (s...) - L = prod(s) +############### +## reducedim ## +############### - if s != size(out) - error("Dimension mismatch") - end +@inline reducedim(op, a::StaticArray, ::Val{D}) where {D} = mapreducedim(identity, op, a, Val{D}) +@inline reducedim(op, a::StaticArray, ::Val{D}, v0) where {D} = mapreducedim(identity, op, a, Val{D}, v0) - exprs = Vector{Expr}(L) +####################### +## related functions ## +####################### - i = 1 - ind = ones(Int, ndims) - while i <= L - ind1 = [expands1[j] ? ind[j] : 1 for j = 1:length(s1)] - ind2 = [expands2[j] ? ind[j] : 1 for j = 1:length(s2)] - index1 = sub2ind(s1, ind1...) - index2 = sub2ind(s2, ind2...) +# These are all similar in Base but not @inline'd +@inline sum(a::StaticArray{T}) where {T} = reduce(+, zero(T), a) +@inline prod(a::StaticArray{T}) where {T} = reduce(*, one(T), a) +@inline count(a::StaticArray{Bool}) = reduce(+, 0, a) +@inline all(a::StaticArray{Bool}) = reduce(&, true, a) # non-branching versions +@inline any(a::StaticArray{Bool}) = reduce(|, false, a) # (benchmarking needed) +@inline mean(a::StaticArray) = sum(a) / length(a) +@inline sumabs(a::StaticArray{T}) where {T} = mapreduce(abs, +, zero(T), a) +@inline sumabs2(a::StaticArray{T}) where {T} = mapreduce(abs2, +, zero(T), a) +@inline minimum(a::StaticArray) = reduce(min, a) # base has mapreduce(idenity, scalarmin, a) +@inline maximum(a::StaticArray) = reduce(max, a) # base has mapreduce(idenity, scalarmax, a) +@inline minimum(a::StaticArray, dim::Type{Val{D}}) where {D} = reducedim(min, a, dim) +@inline maximum(a::StaticArray, dim::Type{Val{D}}) where {D} = reducedim(max, a, dim) - exprs[i] = :(out[$i] = $(Expr(:call, :f, Expr(:ref, :a1, index1), Expr(:ref, :a2, index2)))) +# Diff is slightly different +@inline diff(a::StaticArray) = diff(a, Val{1}) +@inline diff(a::StaticArray, ::Type{Val{D}}) where {D} = _diff(Size(a), a, Val{D}) - i += 1 +@generated function _diff(::Size{S}, a::StaticArray, ::Type{Val{D}}) where {S, D} + N = length(S) + Snew = ([n==D ? S[n]-1 : S[n] for n = 1:N]...) - ind[1] += 1 - j = 1 - while j < length(s) - if ind[j] > s[j] - ind[j] = 1 - ind[j+1] += 1 - else - break - end + exprs = Array{Expr}(Snew) + itr = [1:n for n = Snew] - j += 1 - end - end + for i1 = Base.product(itr...) + i2 = copy([i1...]) + i2[D] = i1[D] + 1 + exprs[i1...] = :(a[$(i2...)] - a[$(i1...)]) + end - return quote - $(Expr(:meta, :inline)) - @inbounds $(Expr(:block, exprs...)) - end + # TODO element type might change + return quote + @_inline_meta + @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...))) end end diff --git a/src/matrix_multiply.jl b/src/matrix_multiply.jl index e2a5b2f2..da1fef23 100644 --- a/src/matrix_multiply.jl +++ b/src/matrix_multiply.jl @@ -1,676 +1,274 @@ 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 + +const StaticVecOrMat{T} = Union{StaticVector{T}, StaticMatrix{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... # TODO make faster versions of A*_mul_B* -@generated function A_mul_Bc(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A * ctranspose(B) - end -end -@generated function Ac_mul_B(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - ctranspose(A) * B - end -end -@generated function Ac_mul_Bc(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - ctranspose(B*A) # is this always safe? - end -end +@inline A_mul_Bc(A::StaticVecOrMat, B::StaticVecOrMat) = A * ctranspose(B) +@inline Ac_mul_Bc(A::StaticVecOrMat, B::StaticVecOrMat) = ctranspose(A) * ctranspose(B) +@inline Ac_mul_B(A::StaticVecOrMat, B::StaticVecOrMat) = ctranspose(A) * B -@generated function A_mul_Bt(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A * transpose(B) - end -end -@generated function At_mul_B(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - transpose(A) * B - end -end -@generated function At_mul_Bt(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - transpose(B*A) # is this always safe? - end -end +@inline A_mul_Bt(A::StaticVecOrMat, B::StaticVecOrMat) = A * transpose(B) +@inline At_mul_Bt(A::StaticVecOrMat, B::StaticVecOrMat) = transpose(A) * transpose(B) +@inline At_mul_B(A::StaticVecOrMat, B::StaticVecOrMat) = transpose(A) * B -# mutating +@inline A_mul_Bc!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, A, ctranspose(B)) +@inline Ac_mul_Bc!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, ctranspose(A), ctranspose(B)) +@inline Ac_mul_B!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, ctranspose(A), B) -@generated function A_mul_Bc!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, A, ctranspose(B)) - end -end -@generated function Ac_mul_B!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, ctranspose(A), B) - end -end -@generated function Ac_mul_Bc!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, ctranspose(A), ctranspose(B)) - end -end +@inline A_mul_Bt!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, A, transpose(B)) +@inline At_mul_Bt!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, transpose(A), transpose(B)) +@inline At_mul_B!(dest::StaticVecOrMat, A::StaticVecOrMat, B::StaticVecOrMat) = A_mul_B!(dest, transpose(A), B) -@generated function A_mul_Bt!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, A, transpose(B)) - end -end -@generated function At_mul_B!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, transpose(A), B) - end -end -@generated function At_mul_Bt!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector}) - return quote - $(Expr(:meta, :inline)) - A_mul_B!(C, transpose(A), transpose(B)) - end -end +# Manage dispatch of * and A_mul_B! +# TODO RowVector? (Inner product?) +@inline *(A::StaticMatrix, B::StaticVector) = _A_mul_B(Size(A), Size(B), A, B) +@inline *(A::StaticMatrix, B::StaticMatrix) = _A_mul_B(Size(A), Size(B), A, B) +@inline *(A::StaticVector, B::StaticMatrix) = *(reshape(A, Size(Size(A)[1], 1)), B) -@generated function *{TA,Tb}(A::StaticMatrix{TA}, b::StaticVector{Tb}) - sA = size(A) - sb = size(b) +@inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticVector) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B) +@inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticMatrix) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B) +@inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::StaticMatrix) = A_mul_B!(dest, reshape(A, Size(Size(A)[1], 1)), B) - s = Size(sA[1]) - T = promote_matprod(TA, Tb) - #println(T) +#@inline *{TA<:Base.LinAlg.BlasFloat,Tb}(A::StaticMatrix{TA}, b::StaticVector{Tb}) - if sb[1] != sA[2] - error("Dimension mismatch") - end +# Implementations - if s == sb - if T == Tb - newtype = b - else - newtype = similar_type(b, T) - end - else - if T == Tb - newtype = similar_type(b, s) - else - newtype = similar_type(b, T, s) - end +@generated function _A_mul_B(::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticVector{Tb}) where {sa, sb, Ta, Tb} + if sb[1] != sa[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) 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]] + 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(T) for k = 1:sA[1]] + exprs = [:(zero(T)) for k = 1:sa[1]] end return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + T = promote_matprod(Ta, Tb) + @inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...))) end end -# For an ambiguity relating to the below two functions -@generated function *{TA<:Base.LinAlg.BlasFloat,Tb}(A::StaticMatrix{TA}, b::StaticVector{Tb}) - sA = size(A) - sb = size(b) - - s = Size(sA[1]) - T = promote_matprod(TA, Tb) - #println(T) - - if sb[1] != sA[2] - error("Dimension mismatch") - end - - if s == sb - if T == Tb - newtype = b - else - newtype = similar_type(b, T) - end - else - if T == Tb - newtype = similar_type(b, s) - else - newtype = similar_type(b, T, s) - end - 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(T) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -# This happens to be size-inferrable from A -@generated function *{TA,Tb}(A::StaticMatrix{TA}, b::AbstractVector{Tb}) - sA = size(A) - #sb = size(b) - - s = Size(sA[1]) - T = promote_matprod(TA, Tb) - - if T == TA - newtype = similar_type(A, s) - else - newtype = similar_type(A, T, s) - 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(T) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - if length(b) != $(sA[2]) - error("Dimension mismatch") - end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -# Ambiguity with BLAS function -@generated function *{TA <: Base.LinAlg.BlasFloat, Tb}(A::StaticMatrix{TA}, b::StridedVector{Tb}) - sA = size(A) - - s = Size(sA[1]) - T = promote_matprod(TA, Tb) - - if T == TA - newtype = similar_type(A, s) - else - newtype = similar_type(A, T, s) - 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(T) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - if length(b) != $(sA[2]) - error("Dimension mismatch") - end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function *(a::StaticVector, B::StaticMatrix) - Ta = eltype(a) - TB = eltype(B) - sa = size(a) - sB = size(B) - - s = Size(sa[1],sB[2]) - T = promote_matprod(Ta, TB) - - if sB[1] != 1 - error("Dimension mismatch") - end - - if s == sB - if T == TB - newtype = B - else - newtype = similar_type(B, T) - end - else - if T == TB - newtype = similar_type(B, s) - else - newtype = similar_type(B, T, s) - end - end - - exprs = [:(a[$j] * B[$k]) for j = 1:s[1], k = 1:s[2]] - - return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -#@inline *{S1,S2,S3}(A::MMatrix{S1,S2}, B::MMatrix{S2,S3}) = MMatrix{S1,S3}(SMatrix{S1,S2}(A)*SMatrix{S2,S3}(B)) - -@generated function *(A::StaticMatrix, B::StaticMatrix) - TA = eltype(A) - TB = eltype(B) +# TODO: I removed StaticMatrix * AbstractVector. Reinstate? - T = promote_matprod(TA, TB) - - can_mutate = !isbits(A) || !isbits(B) # !isbits implies can get a persistent pointer (to pass to BLAS). Probably will change to !isimmutable in a future version of Julia. - can_blas = T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}} +@generated function _A_mul_B(Sa::Size{sa}, Sb::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} + can_mutate = a.mutable && b.mutable # TODO this probably isn't safe. Maybe a trait?? + can_blas = Ta == Tb && Ta <: BlasFloat if can_mutate - sA = size(A) - sB = size(B) - s = Size(sA[1], sB[2]) + S = Size(sa[1], sb[2]) # Heuristic choice between BLAS and explicit unrolling (or chunk-based unrolling) - if can_blas && size(A,1)*size(A,2)*size(B,2) >= 14*14*14 + if can_blas && sa[1]*sa[2]*sb[2] >= 14*14*14 return quote - $(Expr(:meta, :inline)) - C = similar(A, $T, $(Size(s))) - A_mul_B_blas!(C, A, B) + @_inline_meta + T = promote_matprod(Ta, Tb) + C = similar(a, T, $S) + A_mul_B_blas!($S, C, Sa, Sb, a, b) return C end - elseif size(A,1)*size(A,2)*size(B,2) < 8*8*8 + elseif sa[1]*sa[2]*sb[2] < 8*8*8 return quote - $(Expr(:meta, :inline)) - return A_mul_B_unrolled(A, B) + @_inline_meta + return A_mul_B_unrolled(Sa, Sb, a, b) end - elseif size(A,1) <= 14 && size(A,2) <= 14 && size(B,2) <= 14 + elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 return quote - $(Expr(:meta, :inline)) - return $(similar_type(A, T, s))(A_mul_B_unrolled_chunks(A, B)) + @_inline_meta + T = promote_matprod(Ta, Tb) + return similar_type(a, T, $S)(A_mul_B_unrolled_chunks(Sa, Sb, a, b)) end else return quote - $(Expr(:meta, :inline)) - return A_mul_B_loop(A, B) + @_inline_meta + return A_mul_B_loop(Sa, Sb, a, b) end end else # both are isbits type... # Heuristic choice for amount of codegen - if size(A,1)*size(A,2)*size(B,2) <= 8*8*8 + if sa[1]*sa[2]*sb[2] <= 8*8*8 return quote - $(Expr(:meta, :inline)) - return A_mul_B_unrolled(A, B) + @_inline_meta + return A_mul_B_unrolled(Sa, Sb, a, b) end - elseif size(A,1) <= 14 && size(A,2) <= 14 && size(B,2) <= 14 + elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 return quote - $(Expr(:meta, :inline)) - return A_mul_B_unrolled_chunks(A, B) + @_inline_meta + return A_mul_B_unrolled_chunks(Sa, Sb, a, b) end else return quote - $(Expr(:meta, :inline)) - return A_mul_B_loop(A, B) + @_inline_meta + return A_mul_B_loop(Sa, Sb, a, b) end end end end -@generated function A_mul_B_unrolled(A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - - s = Size(sA[1], sB[2]) - T = promote_matprod(TA, TB) - - if sB[1] != sA[2] - error("Dimension mismatch") +@generated function A_mul_B_unrolled(::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} + if sb[1] != sa[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) end - # TODO think about which to be similar to - if s == sB - if T == TB - newtype = B - else - newtype = similar_type(B, T) - end - else - if T == TB - newtype = similar_type(B, s) - else - newtype = similar_type(B, T, s) - end - end + S = Size(sa[1], sb[2]) - if sA[2] != 0 - exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k1, j))]*B[$(sub2ind(sB, j, k2))]) for j = 1:sA[2]]) for k1 = 1:sA[1], k2 = 1:sB[2]] + if sa[2] != 0 + exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(a[$(sub2ind(sa, k1, j))]*b[$(sub2ind(sb, j, k2))]) for j = 1:sa[2]]) for k1 = 1:sa[1], k2 = 1:sb[2]] else - exprs = [zero(T) for k1 = 1:sA[1], k2 = 1:sB[2]] + exprs = [:(zero(T)) for k1 = 1:sa[1], k2 = 1:sb[2]] end return quote - $(Expr(:meta,:inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + T = promote_matprod(Ta, Tb) + @inbounds return similar_type(a, T, $S)(tuple($(exprs...))) end end -@generated function A_mul_B_loop(A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - - s = Size(sA[1], sB[2]) - T = promote_matprod(TA, TB) - - if sB[1] != sA[2] - error("Dimension mismatch") +@generated function A_mul_B_loop(::Size{sa}, ::Size{sb}, b::StaticMatrix{Ta}, a::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} + if sb[1] != sa[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) end - # TODO think about which to be similar to - if s == sB - if T == TB - newtype = B - else - newtype = similar_type(B, T) - end - else - if T == TB - newtype = similar_type(B, s) - else - newtype = similar_type(B, T, s) - end - end + S = Size(sa[1], sb[2]) - tmps = [Symbol("tmp_$(k1)_$(k2)") for k1 = 1:sA[1], k2 = 1:sB[2]] - exprs_init = [:($(tmps[k1,k2]) = A[$k1] * B[1 + $((k2-1) * sB[1])]) for k1 = 1:sA[1], k2 = 1:sB[2]] - exprs_loop = [:($(tmps[k1,k2]) += A[$(k1-sA[1]) + $(sA[1])*j] * B[j + $((k2-1) * sB[1])]) for k1 = 1:sA[1], k2 = 1:sB[2]] + tmps = [Symbol("tmp_$(k1)_$(k2)") for k1 = 1:sa[1], k2 = 1:sb[2]] + exprs_init = [:($(tmps[k1,k2]) = a[$k1] * b[1 + $((k2-1) * sb[1])]) for k1 = 1:sa[1], k2 = 1:sb[2]] + exprs_loop = [:($(tmps[k1,k2]) += a[$(k1-sa[1]) + $(sa[1])*j] * b[j + $((k2-1) * sb[1])]) for k1 = 1:sa[1], k2 = 1:sb[2]] return quote - $(Expr(:meta,:inline)) + @_inline_meta + T = promote_matprod(Ta, Tb) @inbounds $(Expr(:block, exprs_init...)) - for j = 2:$(sA[2]) + for j = 2:$(sa[2]) @inbounds $(Expr(:block, exprs_loop...)) end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, tmps...))) + @inbounds return similar_type(a, T, $S)(tuple($(tmps...))) end end # Concatenate a series of matrix-vector multiplications # Each function is N^2 not N^3 - aids in compile time. -@generated function A_mul_B_unrolled_chunks(A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - - s = Size(sA[1], sB[2]) - T = promote_matprod(TA, TB) - - if sB[1] != sA[2] - error("Dimension mismatch") +@generated function A_mul_B_unrolled_chunks(::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} + if sb[1] != sa[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) end - # TODO think about which to be similar to - if s == sB - if T == TB - newtype = B - else - newtype = similar_type(B, T) - end - else - if T == TB - newtype = similar_type(B, s) - else - newtype = similar_type(B, T, s) - end - end - - #vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, B[:, $k2])) for k2 = 1:sB[2]] + S = Size(sa[1], sb[2]) - # Do a custom B[:, k2] to return a SVector (an isbits type) rather than (possibly) a mutable type. Avoids allocation == faster - tmp_type_in = SVector{sB[1], T} - tmp_type_out = SVector{sA[1], T} - vect_exprs = [:($(Symbol("tmp_$k2"))::$tmp_type_out = partly_unrolled_multiply(A, $(Expr(:call, tmp_type_in, [Expr(:ref, :B, sub2ind(s, i, k2)) for i = 1:sB[1]]...)))::$tmp_type_out) for k2 = 1:sB[2]] + # Do a custom b[:, k2] to return a SVector (an isbits type) rather than (possibly) a mutable type. Avoids allocation == faster + tmp_type_in = :(SVector{$(sa[1]), T}) + tmp_type_out = :(SVector{$(sb[1]), T}) + vect_exprs = [:($(Symbol("tmp_$k2"))::$tmp_type_out = partly_unrolled_multiply(Size(a), Size($(sa[1])), a, $(Expr(:call, tmp_type_in, [Expr(:ref, :b, sub2ind(S, i, k2)) for i = 1:sb[1]]...)))::$tmp_type_out) for k2 = 1:sb[2]] - exprs = [:($(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sA[1], k2 = 1:sB[2]] - - return Expr(:block, - Expr(:meta,:inline), - vect_exprs..., - :(@inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))) - ) -end - -@generated function partly_unrolled_multiply(A::StaticMatrix, b::StaticVector) - TA = eltype(A) - Tb = eltype(b) - sA = size(A) - sb = size(b) - - s = Size(sA[1]) - T = promote_matprod(TA, Tb) - - if sb[1] != sA[2] - error("Dimension mismatch") - end - - if s == sb - if T == Tb - newtype = b - else - newtype = similar_type(b, T) - end - else - if T == Tb - newtype = similar_type(b, s) - else - newtype = similar_type(b, T, s) - end - 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(T) for k = 1:sA[1]] - end + exprs = [:($(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sa[1], k2 = 1:sb[2]] return quote - $(Expr(:meta,:noinline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + T = promote_matprod(Ta, Tb) + $(Expr(:block, + vect_exprs..., + :(@inbounds return similar_type(a, T, $S)(tuple($(exprs...)))) + )) end end - -# TODO aliasing problems if c === b? -@generated function A_mul_B!{T1,T2,T3}(c::StaticVector{T1}, A::StaticMatrix{T2}, b::StaticVector{T3}) - sA = size(A) - sb = size(b) - s = size(c) - - if sb[1] != sA[2] || s[1] != sA[1] - error("Dimension mismatch") +@generated function partly_unrolled_multiply(::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticVector{Tb}) where {sa, sb, Ta, Tb} + if sa[2] != sb[1] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) end - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] + 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 = [:(c[$k] = $(zero(T1))) for k = 1:sA[1]] + exprs = [:(zero(promote_matprod(Ta,Tb))) for k = 1:sa[1]] end return quote - $(Expr(:meta,:inline)) - @inbounds $(Expr(:block, exprs...)) - end -end - -# These two for ambiguity with a BLAS calling function -@generated function A_mul_B!{T<:Union{Float32, Float64}}(c::StaticVector{T}, A::StaticMatrix{T}, b::StaticVector{T}) - sA = size(A) - sb = size(b) - s = size(c) - - if sb[1] != sA[2] || s[1] != sA[1] - error("Dimension mismatch") - end - - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] - else - exprs = [:(c[$k] = $(zero(T))) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - @inbounds $(Expr(:block, exprs...)) - end -end -@generated function A_mul_B!{T<:Union{Complex{Float32}, Complex{Float64}}}(c::StaticVector{T}, A::StaticMatrix{T}, b::StaticVector{T}) - sA = size(A) - sb = size(b) - s = size(c) - - if sb[1] != sA[2] || s[1] != sA[1] - error("Dimension mismatch") - end - - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] - else - exprs = [:(c[$k] = $(zero(T))) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - @inbounds $(Expr(:block, exprs...)) - end -end - -# The unrolled code is inferrable from the size of A -@generated function A_mul_B!{T1,T2,T3}(c::AbstractVector{T1}, A::StaticMatrix{T2}, b::AbstractVector{T3}) - sA = size(A) - T = eltype(c) - - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] - else - exprs = [:(c[$k] = $(zero(T))) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - if length(b) != $(sA[2]) || length(c) != $(sA[1]) - error("Dimension mismatch") - end - @inbounds $(Expr(:block, exprs...)) + $(Expr(:meta,:noinline)) + @inbounds return SVector(tuple($(exprs...))) end end -# Ambiguity with a BLAS specialized function -# Also possible bug makes this harder to resolve (see https://github.com/JuliaLang/julia/issues/19124) -# (problem being that I can't use T<:BlasFloat) -@generated function A_mul_B!{T<:Union{Float64,Float32}}(c::StridedVector{T}, A::StaticMatrix{T}, b::StridedVector{T}) - sA = size(A) - - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] - else - exprs = [:(c[$k] = $(zero(T))) for k = 1:sA[1]] - end - - return quote - $(Expr(:meta,:inline)) - if length(b) != $(sA[2]) || length(c) != $(sA[1]) - error("Dimension mismatch") - end - @inbounds $(Expr(:block, exprs...)) +# TODO aliasing problems if c === b? +@generated function _A_mul_B!(::Size{sc}, c::StaticVector, ::Size{sa}, ::Size{sb}, a::StaticMatrix, b::StaticVector) where {sa, sb, sc} + if sb[1] != sa[2] || sc[1] != sa[1] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) end -end -@generated function A_mul_B!{T<:Union{Complex{Float64},Complex{Float32}}}(c::StridedVector{T}, A::StaticMatrix{T}, b::StridedVector{T}) - sA = size(A) - if sA[2] != 0 - exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]))) for k = 1:sA[1]] + if sa[2] != 0 + exprs = [:(c[$k] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(a[$(sub2ind(sa, k, j))]*b[$j]) for j = 1:sa[2]]))) for k = 1:sa[1]] else - exprs = [:(c[$k] = $(zero(T))) for k = 1:sA[1]] + exprs = [:(c[$k] = zero(eltype(c))) for k = 1:sa[1]] end return quote - $(Expr(:meta,:inline)) - if length(b) != $(sA[2]) || length(c) != $(sA[1]) - error("Dimension mismatch") - end + @_inline_meta @inbounds $(Expr(:block, exprs...)) end end - -@generated function A_mul_B!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix) - if isbits(C) - error("Cannot mutate $C") - end - - TA = eltype(A) - TB = eltype(B) - T = promote_matprod(TA, TB) - - can_blas = T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}} +@generated function _A_mul_B!(::Size{sc}, c::StaticMatrix{Tc}, ::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, sc, Ta, Tb, Tc} + can_blas = Tc == Ta && Tc == Tb && Tc <: BlasFloat if can_blas - if size(A,1) * size(A,2) * size(A,3) < 4*4*4 + if sa[1] * sa[2] * sb[2] < 4*4*4 return quote - $(Expr(:meta, :inline)) - A_mul_B_unrolled!(C, A, B) - return C + @_inline_meta + A_mul_B_unrolled!(c, a, b) + return c end - elseif size(A,1) * size(A,2) * size(A,3) < 14*14*14 # Something seems broken for this one with large matrices (becomes allocating) + elseif sa[1] * sa[2] * sb[2] < 14*14*14 # Something seems broken for this one with large matrices (becomes allocating) return quote - $(Expr(:meta, :inline)) - A_mul_B_unrolled_chunks!(C, A, B) - return C + @_inline_meta + A_mul_B_unrolled_chunks!(c, a, b) + return c end else return quote - $(Expr(:meta, :inline)) - A_mul_B_blas!(C, A, B) - return C + @_inline_meta + A_mul_B_blas!(c, a, b) + return c end end else - if size(A,1) * size(A,2) * size(A,3) < 4*4*4 + if sa[1] * sa[2] * sb[2] < 4*4*4 return quote - $(Expr(:meta, :inline)) - A_mul_B_unrolled!(C, A, B) - return C + @_inline_meta + A_mul_B_unrolled!(c, a, b) + return c end else return quote - $(Expr(:meta, :inline)) - A_mul_B_unrolled_chunks!(C, A, B) - return C + @_inline_meta + A_mul_B_unrolled_chunks!(c, a, b) + return c end end end end -@generated function A_mul_B_blas!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - s = size(C) - T = eltype(C) - - if sB[1] != sA[2] || sA[1] != s[1] || sB[2] != s[2] - error("Dimension mismatch") +@generated function A_mul_B_blas!(::Size{s}, c::StaticMatrix{T}, ::Size{sa}, ::Size{sb}, a::StaticMatrix{T}, b::StaticMatrix{T}) where {s,sa,sb, T <: BlasFloat} + if sb[1] != sa[2] || sa[1] != s[1] || sb[2] != s[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $s")) end - if sA[1] > 0 && sA[2] > 0 && sB[2] > 0 && T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}} + if sa[1] > 0 && sa[2] > 0 && sb[2] > 0 # This code adapted from `gemm!()` in base/linalg/blas.jl if T == Float64 @@ -684,16 +282,16 @@ end end return quote - alpha = $(one(T)) - beta = $(zero(T)) + alpha = one(T) + beta = zero(T) transA = 'N' transB = 'N' - m = $(sA[1]) - ka = $(sA[2]) - kb = $(sB[1]) - n = $(sB[2]) - strideA = $(sA[1]) - strideB = $(sB[1]) + m = $(sa[1]) + ka = $(sa[2]) + kb = $(sb[1]) + n = $(sb[2]) + strideA = $(sa[1]) + strideB = $(sb[1]) strideC = $(s[1]) ccall((Base.BLAS.@blasfunc($gemm), Base.BLAS.libblas), Void, @@ -702,68 +300,49 @@ end Ptr{$T}, Ptr{Base.BLAS.BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{Base.BLAS.BlasInt}), &transA, &transB, &m, &n, - &ka, &alpha, A, &strideA, - B, &strideB, &beta, C, + &ka, &alpha, a, &strideA, + b, &strideB, &beta, c, &strideC) - return C + return c end else - error("Cannot call BLAS gemm with $C = $A * $B") + throw(DimensionMismatch("Cannot call BLAS gemm with zero-dimension arrays, attempted $sa * $sb -> $s.")) end end -@generated function A_mul_B_unrolled!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - T = eltype(C) - - s = (sA[1], sB[2]) - - if sB[1] != sA[2] - error("Dimension mismatch") - end - if s != size(C) - error("Dimension mismatch") +@generated function A_mul_B_unrolled!(::Size{sc}, c::StaticMatrix, ::Size{sa}, ::Size{sb}, a::StaticMatrix, b::StaticMatrix) where {sa, sb, sc} + if sb[1] != sa[2] || sa[1] != s[1] || sb[2] != sc[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) end - if sA[2] != 0 - exprs = [:(C[$(sub2ind(s, k1, k2))] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k1, j))]*B[$(sub2ind(sB, j, k2))]) for j = 1:sA[2]]))) for k1 = 1:sA[1], k2 = 1:sB[2]] + if sa[2] != 0 + exprs = [:(c[$(sub2ind(sc, k1, k2))] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(a[$(sub2ind(sa, k1, j))]*b[$(sub2ind(sb, j, k2))]) for j = 1:sa[2]]))) for k1 = 1:sa[1], k2 = 1:sb[2]] else - exprs = [:(C[$(sub2ind(s, k1, k2))] = $(zero(T))) for k1 = 1:sA[1], k2 = 1:sB[2]] + exprs = [:(c[$(sub2ind(sc, k1, k2))] = zero(eltype(c))) for k1 = 1:sa[1], k2 = 1:sb[2]] end return quote - $(Expr(:meta,:inline)) + @_inline_meta @inbounds $(Expr(:block, exprs...)) end end -@generated function A_mul_B_unrolled_chunks!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix) - sA = size(A) - sB = size(B) - TA = eltype(A) - TB = eltype(B) - - s = size(C) - T = eltype(C) - - if sB[1] != sA[2] || sA[1] != s[1] || sB[2] != s[2] - error("Dimension mismatch") +@generated function A_mul_B_unrolled_chunks!(::Size{sc}, c::StaticMatrix, ::Size{sa}, ::Size{sb}, a::StaticMatrix, b::StaticMatrix) where {sa, sb, sc} + if sb[1] != sa[2] || sa[1] != s[1] || sb[2] != sc[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) end #vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, B[:, $k2])) for k2 = 1:sB[2]] - # Do a custom B[:, k2] to return a SVector (an isbits type) rather than a mutable type. Avoids allocation == faster - tmp_type = SVector{sB[1], T} - vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, $(Expr(:call, tmp_type, [Expr(:ref, :B, sub2ind(s, i, k2)) for i = 1:sB[1]]...)))) for k2 = 1:sB[2]] + # Do a custom b[:, k2] to return a SVector (an isbits type) rather than a mutable type. Avoids allocation == faster + tmp_type = SVector{sb[1], eltype(c)} + vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(a, $(Expr(:call, tmp_type, [Expr(:ref, :b, sub2ind(s, i, k2)) for i = 1:sb[1]]...)))) for k2 = 1:sb[2]] - exprs = [:(C[$(sub2ind(s, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sA[1], k2 = 1:sB[2]] + exprs = [:(c[$(sub2ind(s, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sa[1], k2 = 1:sb[2]] return quote - Expr(:meta,:inline) + @_inline_meta @inbounds $(Expr(:block, vect_exprs...)) @inbounds $(Expr(:block, exprs...)) end diff --git a/src/util.jl b/src/util.jl index cabc174f..d8bf6a93 100644 --- a/src/util.jl +++ b/src/util.jl @@ -1,5 +1,5 @@ # For convenience -@compat TupleN{T,N} = NTuple{N,T} +TupleN{T,N} = NTuple{N,T} # Cast any Tuple to an TupleN{T} @inline convert_ntuple{T}(::Type{T},d::T) = T # For zero-dimensional arrays @@ -28,6 +28,9 @@ end end end + +# TODO: the below seems to be type piracy... +#= # some convenience functions for non-static arrays, generators, etc... @inline convert{T}(::Type{Tuple}, a::AbstractArray{T}) = (a...)::Tuple{Vararg{T}} @inline function convert{N,T}(::Type{NTuple{N,Any}}, a::AbstractArray{T}) @@ -46,6 +49,7 @@ end @inbounds return ntuple(i -> convert(T1,a[i]), Val{N}) end + if VERSION < v"0.5+" # TODO try and make this generate fast code @inline convert(::Type{Tuple}, g::Base.Generator) = (g...) @@ -57,7 +61,7 @@ if VERSION < v"0.5+" @inbounds return ntuple(i -> g.f(g.iter[i]), Val{N}) end end - +=# #= @generated function convert{N}(::Type{NTuple{N,Any}}, g::Base.Generator) exprs = [:(g.f(g.iter[$j])) for j=1:N] diff --git a/test/abstractarray.jl b/test/abstractarray.jl index c73653e7..1e788b3d 100644 --- a/test/abstractarray.jl +++ b/test/abstractarray.jl @@ -25,18 +25,18 @@ @test @inferred(similar_type(SVector{2,Int}, Float64, Size(3,3,3))) == SArray{(3,3,3), Float64, 3, 27} # Some specializations for the mutable case - @test @inferred(similar_type(MVector{3,Int}, Float64)) == MVector{3,Float64} - @test @inferred(similar_type(MMatrix{3,3,Int,9}, Size(2))) == MVector{2, Int} - @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64, Size(2))) == MVector{2, Float64} - @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64, Size(2))) == MVector{2, Float64} + @test @inferred(similar_type(MVector{3,Int}, Float64)) == SVector{3,Float64} + @test @inferred(similar_type(MMatrix{3,3,Int,9}, Size(2))) == SVector{2, Int} + @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64, Size(2))) == SVector{2, Float64} + @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64, Size(2))) == SVector{2, Float64} - @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64)) == MMatrix{3, 3, Float64, 9} - @test @inferred(similar_type(MVector{2,Int}, Size(3,3))) == MMatrix{3, 3, Int, 9} - @test @inferred(similar_type(MVector{2,Int}, Float64, Size(3,3))) == MMatrix{3, 3, Float64, 9} + @test @inferred(similar_type(MMatrix{3,3,Int,9}, Float64)) == SMatrix{3, 3, Float64, 9} + @test @inferred(similar_type(MVector{2,Int}, Size(3,3))) == SMatrix{3, 3, Int, 9} + @test @inferred(similar_type(MVector{2,Int}, Float64, Size(3,3))) == SMatrix{3, 3, Float64, 9} - @test @inferred(similar_type(MArray{(4,4,4),Int,3,64}, Float64)) == MArray{(4,4,4), Float64, 3, 64} - @test @inferred(similar_type(MVector{2,Int}, Size(3,3,3))) == MArray{(3,3,3), Int, 3, 27} - @test @inferred(similar_type(MVector{2,Int}, Float64, Size(3,3,3))) == MArray{(3,3,3), Float64, 3, 27} + @test @inferred(similar_type(MArray{(4,4,4),Int,3,64}, Float64)) == SArray{(4,4,4), Float64, 3, 64} + @test @inferred(similar_type(MVector{2,Int}, Size(3,3,3))) == SArray{(3,3,3), Int, 3, 27} + @test @inferred(similar_type(MVector{2,Int}, Float64, Size(3,3,3))) == SArray{(3,3,3), Float64, 3, 27} end @testset "similar" begin diff --git a/test/arraymath.jl b/test/arraymath.jl index a3d054fe..d28b1cff 100644 --- a/test/arraymath.jl +++ b/test/arraymath.jl @@ -25,8 +25,6 @@ m1 = @SMatrix [1 2; 3 4] m2 = @SMatrix [4 3; 2 1] - @test @inferred(.-(m1)) === @SMatrix [-1 -2; -3 -4] - @test @inferred(m1 .+ m2) === @SMatrix [5 5; 5 5] @test @inferred(m1 .* m2) === @SMatrix [4 6; 6 4] @test @inferred(m1 .- m2) === @SMatrix [-3 -1; 1 3] diff --git a/test/indexing.jl b/test/indexing.jl index d2b70a30..41003983 100644 --- a/test/indexing.jl +++ b/test/indexing.jl @@ -2,7 +2,7 @@ @testset "Linear getindex() on SVector" begin sv = SVector{4}(4,5,6,7) - # Tuple + # SVector @test (@inferred getindex(sv, SVector(4,3,2,1))) === SVector((7,6,5,4)) # Colon @@ -12,13 +12,13 @@ @testset "Linear getindex()/setindex!() on MVector" begin vec = @SVector [4,5,6,7] - # Tuple + # SVector mv = MVector{4,Int}() - @test (mv[(1,2,3,4)] = vec; (@inferred getindex(mv, (4,3,2,1)))::MVector{4,Int} == MVector((7,6,5,4))) + @test (mv[SVector(1,2,3,4)] = vec; (@inferred getindex(mv, SVector(4,3,2,1)))::SVector{4,Int} == SVector((7,6,5,4))) # Colon mv = MVector{4,Int}() - @test (mv[:] = vec; (@inferred getindex(mv, :))::MVector{4,Int} == MVector((4,5,6,7))) + @test (mv[:] = vec; (@inferred getindex(mv, :))::SVector{4,Int} == SVector((4,5,6,7))) end @testset "Linear getindex()/setindex!() with a SVector on an Array" begin @@ -36,15 +36,15 @@ @test sm[1,2] === 3 @test sm[2,2] === 4 - # Tuple, scalar - @test (@inferred getindex(sm, (2,1), (2,1))) === @SMatrix [4 2; 3 1] - @test (@inferred getindex(sm, 1, (1,2))) === @SVector [1,3] - @test (@inferred getindex(sm, (1,2), 1)) === @SVector [1,2] + # SVector, scalar + @test (@inferred getindex(sm, SVector(2,1), SVector(2,1))) === @SMatrix [4 2; 3 1] + @test (@inferred getindex(sm, 1, SVector(1,2))) === @SVector [1,3] + @test (@inferred getindex(sm, SVector(1,2), 1)) === @SVector [1,2] # Colon @test (@inferred getindex(sm, :, :)) === @SMatrix [1 3; 2 4] - @test (@inferred getindex(sm, (2,1), :)) === @SMatrix [2 4; 1 3] - @test (@inferred getindex(sm, :, (2,1))) === @SMatrix [3 1; 4 2] + @test (@inferred getindex(sm, SVector(2,1), :)) === @SMatrix [2 4; 1 3] + @test (@inferred getindex(sm, :, SVector(2,1))) === @SMatrix [3 1; 4 2] @test (@inferred getindex(sm, 1, :)) === @SVector [1,3] @test (@inferred getindex(sm, :, 1)) === @SVector [1,2] end @@ -53,16 +53,16 @@ sm = @MMatrix [1 3; 2 4] # Tuple, scalar - @test (mm = MMatrix{2,2,Int}(); mm[(2,1),(2,1)] = sm[(2,1),(2,1)]; (@inferred getindex(mm, (2,1), (2,1)))::MMatrix == @MMatrix [4 2; 3 1]) - @test (mm = MMatrix{2,2,Int}(); mm[1,(1,2)] = sm[1,(1,2)]; (@inferred getindex(mm, 1, (1,2)))::MVector == @MVector [1,3]) - @test (mm = MMatrix{2,2,Int}(); mm[(1,2),1] = sm[(1,2),1]; (@inferred getindex(mm, (1,2), 1))::MVector == @MVector [1,2]) + @test (mm = MMatrix{2,2,Int}(); mm[SVector(2,1),SVector(2,1)] = sm[SVector(2,1),SVector(2,1)]; (@inferred getindex(mm, SVector(2,1), SVector(2,1)))::SMatrix == @SMatrix [4 2; 3 1]) + @test (mm = MMatrix{2,2,Int}(); mm[1,SVector(1,2)] = sm[1,SVector(1,2)]; (@inferred getindex(mm, 1, SVector(1,2)))::SVector == @SVector [1,3]) + @test (mm = MMatrix{2,2,Int}(); mm[SVector(1,2),1] = sm[SVector(1,2),1]; (@inferred getindex(mm, SVector(1,2), 1))::SVector == @SVector [1,2]) # Colon - @test (mm = MMatrix{2,2,Int}(); mm[:,:] = sm[:,:]; (@inferred getindex(mm, :, :))::MMatrix == @MMatrix [1 3; 2 4]) - @test (mm = MMatrix{2,2,Int}(); mm[(2,1),:] = sm[(2,1),:]; (@inferred getindex(mm, (2,1), :))::MMatrix == @MMatrix [2 4; 1 3]) - @test (mm = MMatrix{2,2,Int}(); mm[:,(2,1)] = sm[:,(2,1)]; (@inferred getindex(mm, :, (2,1)))::MMatrix == @MMatrix [3 1; 4 2]) - @test (mm = MMatrix{2,2,Int}(); mm[1,:] = sm[1,:]; (@inferred getindex(mm, 1, :))::MVector == @MVector [1,3]) - @test (mm = MMatrix{2,2,Int}(); mm[:,1] = sm[:,1]; (@inferred getindex(mm, :, 1))::MVector == @MVector [1,2]) + @test (mm = MMatrix{2,2,Int}(); mm[:,:] = sm[:,:]; (@inferred getindex(mm, :, :))::SMatrix == @MMatrix [1 3; 2 4]) + @test (mm = MMatrix{2,2,Int}(); mm[SVector(2,1),:] = sm[SVector(2,1),:]; (@inferred getindex(mm, SVector(2,1), :))::SMatrix == @SMatrix [2 4; 1 3]) + @test (mm = MMatrix{2,2,Int}(); mm[:,SVector(2,1)] = sm[:,SVector(2,1)]; (@inferred getindex(mm, :, SVector(2,1)))::SMatrix == @SMatrix [3 1; 4 2]) + @test (mm = MMatrix{2,2,Int}(); mm[1,:] = sm[1,:]; (@inferred getindex(mm, 1, :))::SVector == @SVector [1,3]) + @test (mm = MMatrix{2,2,Int}(); mm[:,1] = sm[:,1]; (@inferred getindex(mm, :, 1))::SVector == @SVector [1,2]) end @testset "3D scalar indexing" begin diff --git a/test/linalg.jl b/test/linalg.jl index da271fc1..4f568fad 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -17,10 +17,10 @@ v3 = [2,4,6,8] v4 = [4,3,2,1] - @test @inferred(v1 - v4) === @SVector [-2, 1, 4, 7] - @test @inferred(v3 - v2) === @SVector [-2, 1, 4, 7] - @test @inferred(v1 - v4) === @SVector [-2, 1, 4, 7] - @test @inferred(v3 - v2) === @SVector [-2, 1, 4, 7] + @test_broken @inferred(v1 + v4) === @SVector [6, 7, 8, 9] + @test_broken @inferred(v3 + v2) === @SVector [6, 7, 8, 9] + @test_broken @inferred(v1 - v4) === @SVector [-2, 1, 4, 7] + @test_broken @inferred(v3 - v2) === @SVector [-2, 1, 4, 7] end @testset "Interaction with `UniformScaling`" begin @@ -73,11 +73,11 @@ @testset "transpose() and conj()" begin @test @inferred(conj(SVector(1+im, 2+im))) === SVector(1-im, 2-im) - @test @inferred(transpose(@SVector([1, 2, 3]))) === @SMatrix([1 2 3]) + @test @inferred(transpose(@SVector([1, 2, 3]))) === RowVector(@SVector([1, 2, 3])) @test @inferred(transpose(@SMatrix([1 2; 0 3]))) === @SMatrix([1 0; 2 3]) @test @inferred(transpose(@SMatrix([1 2 3; 4 5 6]))) === @SMatrix([1 4; 2 5; 3 6]) - @test @inferred(ctranspose(@SVector([1, 2, 3]))) === @SMatrix([1 2 3]) + @test @inferred(ctranspose(@SVector([1, 2, 3]))) === RowVector(@SVector([1, 2, 3])) @test @inferred(ctranspose(@SMatrix([1 2; 0 3]))) === @SMatrix([1 0; 2 3]) @test @inferred(ctranspose(@SMatrix([1 2 3; 4 5 6]))) === @SMatrix([1 4; 2 5; 3 6]) end diff --git a/test/mapreduce.jl b/test/mapreduce.jl index df1ba64b..efa5c984 100644 --- a/test/mapreduce.jl +++ b/test/mapreduce.jl @@ -9,8 +9,8 @@ @test @inferred(map(-, v1)) === @SVector [-2, -4, -6, -8] @test @inferred(map(+, v1, v2)) === @SVector [6, 7, 8, 9] - @test @inferred(map(+, normal_v1, v2)) === @SVector [6, 7, 8, 9] - @test @inferred(map(+, v1, normal_v2)) === @SVector [6, 7, 8, 9] + @test_broken @inferred(map(+, normal_v1, v2)) === @SVector [6, 7, 8, 9] + @test_broken @inferred(map(+, v1, normal_v2)) === @SVector [6, 7, 8, 9] map!(+, mv, v1, v2) @test mv == @MVector [6, 7, 8, 9] diff --git a/test/matrix_multiply.jl b/test/matrix_multiply.jl index 7a3aeb35..ecae737a 100644 --- a/test/matrix_multiply.jl +++ b/test/matrix_multiply.jl @@ -10,27 +10,27 @@ @test x == @SVector [CartesianIndex((7,5)), CartesianIndex((15,13))] v3 = [1, 2] - @test m*v3 === @SVector [5, 11] + @test_broken m*v3 === @SVector [5, 11] m2 = @MMatrix [1 2; 3 4] v4 = @MVector [1, 2] - @test (m2*v4)::MVector == @MVector [5, 11] + @test (m2*v4)::SVector == @SVector [5, 11] m3 = @SArray [1 2; 3 4] v5 = @SArray [1, 2] - @test m3*v5 === @SArray [5, 11] + @test m3*v5 === @SVector [5, 11] m4 = @MArray [1 2; 3 4] v6 = @MArray [1, 2] - @test (m4*v6)::MArray == @MArray [5, 11] + @test (m4*v6) === @SVector [5, 11] m5 = @SMatrix [1.0 2.0; 3.0 4.0] v7 = [1.0, 2.0] - @test (m5*v7)::SVector ≈ @SVector [5.0, 11.0] + @test_broken (m5*v7)::SVector ≈ @SVector [5.0, 11.0] m6 = @SMatrix Float32[1.0 2.0; 3.0 4.0] v8 = Float64[1.0, 2.0] - @test (m6*v8)::SVector{2,Float64} ≈ @SVector [5.0, 11.0] + @test_broken (m6*v8)::SVector{2,Float64} ≈ @SVector [5.0, 11.0] end @@ -58,15 +58,15 @@ m = @MMatrix [1 2; 3 4] n = @MMatrix [2 3; 4 5] - @test (m*n)::MMatrix == @MMatrix [10 13; 22 29] + @test (m*n) === @SMatrix [10 13; 22 29] m = @SArray [1 2; 3 4] n = @SArray [2 3; 4 5] - @test m*n === @SArray [10 13; 22 29] + @test m*n === @SMatrix [10 13; 22 29] m = @MArray [1 2; 3 4] n = @MArray [2 3; 4 5] - @test (m*n)::MArray == @MArray [10 13; 22 29] + @test (m*n) == @SMatrix [10 13; 22 29] # Alternative methods used between 8 < n <= 14 and n > 14 m_array = rand(1:10, 10, 10) @@ -92,7 +92,7 @@ m = MMatrix{10,10}(m_array) n = MMatrix{10,10}(n_array) - @test (m*n)::MMatrix == a_array + @test (m*n)::SMatrix == a_array m_array = rand(1:10, 16, 16) n_array = rand(1:10, 16, 16) @@ -100,7 +100,7 @@ m = MMatrix{16,16}(m_array) n = MMatrix{16,16}(n_array) - @test (m*n)::MMatrix == a_array + @test (m*n)::SMatrix == a_array # Mutating BLAS types follow yet different behaviour m_array = randn(4, 4) @@ -128,6 +128,7 @@ @test m*n::MMatrix ≈ a_array end +#= @testset "A_mul_B!" begin v = @SVector [2, 4] v2 = [2, 4] @@ -215,5 +216,5 @@ outvec2f = Vector{Float64}(2) A_mul_B!(outvec2f, mf, vf2) @test outvec2f ≈ [10.0, 22.0] - end + end =# end diff --git a/test/runtests.jl b/test/runtests.jl index 1407cf77..98b64fe8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,5 @@ using StaticArrays using Base.Test -using Compat @testset "StaticArrays" begin include("SVector.jl") @@ -18,13 +17,11 @@ using Compat include("mapreduce.jl") include("arraymath.jl") include("linalg.jl") - include("matrix_multiply.jl") + include("matrix_multiply.jl") #= include("det.jl") include("inv.jl") include("solve.jl") include("eigen.jl") include("deque.jl") - if VERSION < v"0.6.0-dev.1671" - include("fixed_size_arrays.jl") - end + #include("fixed_size_arrays.jl") =# end From 4773d4c204e837eed44552cc166f2376a1399daa Mon Sep 17 00:00:00 2001 From: Chris Foster Date: Mon, 20 Mar 2017 23:08:42 +1000 Subject: [PATCH 3/9] Fix minor bug in matrix multiply --- src/matrix_multiply.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/matrix_multiply.jl b/src/matrix_multiply.jl index da1fef23..fca8ae9b 100644 --- a/src/matrix_multiply.jl +++ b/src/matrix_multiply.jl @@ -138,7 +138,7 @@ end end -@generated function A_mul_B_loop(::Size{sa}, ::Size{sb}, b::StaticMatrix{Ta}, a::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} +@generated function A_mul_B_loop(::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} if sb[1] != sa[2] throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb")) end From 28f405a08d3b70dc906341b37b3774abc4fd9fee Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Tue, 21 Mar 2017 14:27:17 +1000 Subject: [PATCH 4/9] Finished v0.6 fixes (hopefully!) --- src/SVector.jl | 2 +- src/StaticArrays.jl | 16 ++++---- src/abstractarray.jl | 2 + src/cholesky.jl | 6 +-- src/deque.jl | 84 ++++++++++++++++++++++------------------- src/det.jl | 18 ++++----- src/eigen.jl | 19 ++++++---- src/linalg.jl | 1 + src/matrix_multiply.jl | 49 +++++++++++++++++++----- src/solve.jl | 66 +++++++++++++------------------- test/SVector.jl | 1 + test/abstractarray.jl | 2 + test/matrix_multiply.jl | 9 ++++- test/runtests.jl | 4 +- 14 files changed, 155 insertions(+), 124 deletions(-) diff --git a/src/SVector.jl b/src/SVector.jl index a83e12e9..c69152fe 100644 --- a/src/SVector.jl +++ b/src/SVector.jl @@ -25,7 +25,7 @@ immutable SVector{S, T} <: StaticVector{T} end end -@inline (::Type{SVector}){S}(x::NTuple{S}) = SVector{S}(x) +@inline (::Type{SVector}){S}(x::NTuple{S,Any}) = SVector{S}(x) @inline (::Type{SVector{S}}){S, T}(x::NTuple{S,T}) = SVector{S,T}(x) @inline (::Type{SVector{S}}){S, T <: Tuple}(x::T) = SVector{S,promote_tuple_eltype(T)}(x) diff --git a/src/StaticArrays.jl b/src/StaticArrays.jl index 5cb7acf3..39ad7e16 100644 --- a/src/StaticArrays.jl +++ b/src/StaticArrays.jl @@ -4,7 +4,7 @@ module StaticArrays import Base: @_inline_meta, @_propagate_inbounds_meta, @_pure_meta, @propagate_inbounds, @pure -import Base: getindex, setindex!, size, similar, +import Base: getindex, setindex!, size, similar, vec, length, convert, promote_op, map, map!, reduce, reducedim, mapreducedim, mapreduce, broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill, @@ -86,16 +86,14 @@ include("mapreduce.jl") include("arraymath.jl") include("linalg.jl") include("matrix_multiply.jl") -#include("solve.jl") -#include("deque.jl") -#include("det.jl") -#include("inv.jl") -#include("eigen.jl") -#include("cholesky.jl") - +include("det.jl") +include("inv.jl") +include("solve.jl") +include("eigen.jl") +include("cholesky.jl") +include("deque.jl") #include("FixedSizeArrays.jl") # Currently defunct include("ImmutableArrays.jl") - end # module diff --git a/src/abstractarray.jl b/src/abstractarray.jl index 2e618549..ec7c6f9e 100644 --- a/src/abstractarray.jl +++ b/src/abstractarray.jl @@ -123,6 +123,8 @@ end reshape(a::Array, s::Size{S}) where {S} = s(a) +@inline vec(a::StaticArray) = reshape(a, Size(prod(Size(typeof(a))))) + @inline copy(a::StaticArray) = similar_type(a)(Tuple(a)) # TODO permutedims? diff --git a/src/cholesky.jl b/src/cholesky.jl index 50482c6d..cc1d1645 100644 --- a/src/cholesky.jl +++ b/src/cholesky.jl @@ -4,13 +4,9 @@ _chol(Size(A), A) end -@inline function Base.chol{T<:Real, SM <: StaticMatrix}(A::Base.LinAlg.RealHermSymComplexHerm{T,SM}) +@inline function Base.chol(A::Base.LinAlg.RealHermSymComplexHerm{<:Real, <:StaticMatrix}) _chol(Size(A), A.data) end -#= -@inline function Base.chol{T<:Real,SM<:StaticMatrix}(A::Symmetric{T,SM}) - _chol(Size(A), A.data) -end=# @generated function _chol(::Size{(1,1)}, A::StaticMatrix) @assert size(A) == (1,1) diff --git a/src/deque.jl b/src/deque.jl index 73ecaa59..fefce5f5 100644 --- a/src/deque.jl +++ b/src/deque.jl @@ -1,62 +1,68 @@ -@generated function push(vec::StaticVector, x) - newtype = similar_type(vec, Size(length(vec) + 1)) - exprs = vcat([:(vec[$i]) for i = 1:length(vec)], :x) +@inline push(vec::StaticVector, x) = _push(Size(vec), vec, x) +@generated function _push(::Size{s}, vec::StaticVector, x) where {s} + newlen = s[1] + 1 + exprs = vcat([:(vec[$i]) for i = 1:s[1]], :x) return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end -@generated function unshift(vec::StaticVector, x) - newtype = similar_type(vec, Size(length(vec) + 1)) - exprs = vcat(:x, [:(vec[$i]) for i = 1:length(vec)]) +@inline unshift(vec::StaticVector, x) = _unshift(Size(vec), vec, x) +@generated function _unshift(::Size{s}, vec::StaticVector, x) where {s} + newlen = s[1] + 1 + exprs = vcat(:x, [:(vec[$i]) for i = 1:s[1]]) return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end -@generated function insert(vec::StaticVector, index, x) - newtype = similar_type(vec, Size(length(vec) + 1)) +@propagate_inbounds insert(vec::StaticVector, index, x) = _insert(Size(vec), vec, index, x) +@generated function _insert(::Size{s}, vec::StaticVector, index, x) where {s} + newlen = s[1] + 1 exprs = [(i == 1 ? :(ifelse($i < index, vec[$i], x)) : - i == length(vec)+1 ? :(ifelse($i == index, x, vec[$i-1])) : - :(ifelse($i < index, vec[$i], ifelse($i == index, x, vec[$i-1])))) for i = 1:length(vec) + 1] + i == newlen ? :(ifelse($i == index, x, vec[$i-1])) : + :(ifelse($i < index, vec[$i], ifelse($i == index, x, vec[$i-1])))) for i = 1:newlen] return quote - $(Expr(:meta, :inline)) - @boundscheck if (index < 1 || index > $(length(vec)+1)) + @_inline_meta + @boundscheck if (index < 1 || index > $newlen) throw(BoundsError(vec, index)) end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end -@generated function pop(vec::StaticVector) - newtype = similar_type(vec, Size(length(vec) - 1)) - exprs = [:(vec[$i]) for i = 1:length(vec)-1] +@inline pop(vec::StaticVector) = _pop(Size(vec), vec) +@generated function _pop(::Size{s}, vec::StaticVector) where {s} + newlen = s[1] - 1 + exprs = [:(vec[$i]) for i = 1:s[1]-1] return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end -@generated function shift(vec::StaticVector) - newtype = similar_type(vec, Size(length(vec) - 1)) - exprs = [:(vec[$i]) for i = 2:length(vec)] +@inline shift(vec::StaticVector) = _shift(Size(vec), vec) +@generated function _shift(::Size{s}, vec::StaticVector) where {s} + newlen = s[1] - 1 + exprs = [:(vec[$i]) for i = 2:s[1]] return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @_inline_meta + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end -@generated function deleteat(vec::StaticVector, index) - newtype = similar_type(vec, Size(length(vec) - 1)) - exprs = [:(ifelse($i < index, vec[$i], vec[$i+1])) for i = 1:length(vec) - 1] +@propagate_inbounds deleteat(vec::StaticVector, index) = _deleteat(Size(vec), vec, index) +@generated function _deleteat(::Size{s}, vec::StaticVector, index) where {s} + newlen = s[1] - 1 + exprs = [:(ifelse($i < index, vec[$i], vec[$i+1])) for i = 1:newlen] return quote - $(Expr(:meta, :inline)) - @boundscheck if (index < 1 || index > $(length(vec)+1)) + @_inline_meta + @boundscheck if (index < 1 || index > $(s[1])) throw(BoundsError(vec, index)) end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return similar_type(vec, Size($newlen))(tuple($(exprs...))) end end @@ -66,15 +72,15 @@ end # Immutable version of setindex!(). Seems similar in nature to the above, but # could also be justified to live in src/indexing.jl -@generated function setindex{T}(a::StaticArray{T}, x::T, index::Int) - newtype = a - exprs = [:(ifelse($i == index, x, a[$i])) for i = 1:length(a)] +@inline setindex(a::StaticArray{T}, x::T, index::Int) where {T} = _setindex(Size(a), a, x, index) +@generated function _setindex(::Size{s}, a::StaticArray{T}, x::T, index::Int) where {s, T} + exprs = [:(ifelse($i == index, x, a[$i])) for i = 1:s[1]] return quote - $(Expr(:meta, :inline)) - @boundscheck if (index < 1 || index > $(length(a))) + @_inline_meta + @boundscheck if (index < 1 || index > $(s[1])) throw(BoundsError(a, index)) end - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) + @inbounds return typeof(a)(tuple($(exprs...))) end end diff --git a/src/det.jl b/src/det.jl index 3e4a5288..c3164781 100644 --- a/src/det.jl +++ b/src/det.jl @@ -1,39 +1,39 @@ @inline det(A::StaticMatrix) = _det(Size(A), A) -@inline _det(::Size{(1,1)}, A::AbstractMatrix) = @inbounds return A[1] +@inline _det(::Size{(1,1)}, A::StaticMatrix) = @inbounds return A[1] -@inline function _det(::Size{(2,2)}, A::AbstractMatrix) +@inline function _det(::Size{(2,2)}, A::StaticMatrix) @inbounds return A[1]*A[4] - A[3]*A[2] end -@inline function _det{T<:Unsigned}(::Size{(2,2)}, A::AbstractMatrix{T}) +@inline function _det(::Size{(2,2)}, A::StaticMatrix{T}) where {T <: Unsigned} @inbounds return Signed(A[1]*A[4]) - Signed(A[3]*A[2]) end -@inline function _det(::Size{(3,3)}, A::AbstractMatrix) +@inline function _det(::Size{(3,3)}, A::StaticMatrix) @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)) end -@inline function _det{T<:Unsigned}(::Size{(3,3)}, A::AbstractMatrix{T}) +@inline function _det{T<:Unsigned}(::Size{(3,3)}, A::StaticMatrix{T}) @inbounds x0 = SVector(Signed(A[1]), Signed(A[2]), Signed(A[3])) @inbounds x1 = SVector(Signed(A[4]), Signed(A[5]), Signed(A[6])) @inbounds x2 = SVector(Signed(A[7]), Signed(A[8]), Signed(A[9])) return vecdot(x0, cross(x1, x2)) end -@generated function _det{S,T}(::Size{S}, A::AbstractMatrix{T}) +@generated function _det{S,T}(::Size{S}, A::StaticMatrix{T}) if S[1] != S[2] throw(DimensionMismatch("matrix is not square")) end - T2 = typeof((one(T)*zero(T) + zero(T))/one(T)) return quote # Implementation from Base + T2 = typeof((one(T)*zero(T) + zero(T))/one(T)) if istriu(A) || istril(A) - return convert($T2, det(UpperTriangular(A)))::$T2 # Is this a Julia bug that a convert is not type stable?? + return convert(T2, det(UpperTriangular(A))) # Is this a Julia bug that a convert is not type stable?? end - AA = convert(Array{$T2}, A) + AA = convert(Array{T2}, A) return det(lufact(AA)) end end diff --git a/src/eigen.jl b/src/eigen.jl index 41685890..258b98fc 100644 --- a/src/eigen.jl +++ b/src/eigen.jl @@ -1,5 +1,8 @@ -@inline eigvals{T<:Real,SA<:StaticArray}(a::Base.LinAlg.RealHermSymComplexHerm{T,SA},; permute::Bool=true, scale::Bool=true) = _eigvals(Size(SA), a, permute, scale) +@inline function eigvals(a::Base.LinAlg.RealHermSymComplexHerm{T,SA},; permute::Bool=true, scale::Bool=true) where {T <: Real, SA <: StaticArray} + _eigvals(Size(SA), a, permute, scale) +end + @inline function eigvals(a::StaticArray; permute::Bool=true, scale::Bool=true) if ishermitian(a) _eigvals(Size(a), Hermitian(a), permute, scale) @@ -10,7 +13,7 @@ end @inline _eigvals(::Size{(1,1)}, a, permute, scale) = @inbounds return SVector(real(a.data[1])) -@inline function _eigvals{T<:Real}(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eigvals(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} a = A.data if A.uplo == 'U' @@ -30,7 +33,7 @@ end end end -@inline function _eigvals{T<:Real}(::Size{(3,3)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eigvals(::Size{(3,3)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} S = typeof((one(T)*zero(T) + zero(T))/one(T)) Sreal = real(S) @@ -106,7 +109,7 @@ end _eig(Size(A), A, permute, scale) end -@inline function eig{T, SM <: StaticMatrix}(A::Base.LinAlg.HermOrSym{T,SM}; permute::Bool=true, scale::Bool=true) +@inline function eig(A::Base.LinAlg.HermOrSym{<:Any, SM}; permute::Bool=true, scale::Bool=true) where {SM <: StaticMatrix} _eig(Size(SM), A, permute, scale) end @@ -120,18 +123,18 @@ end end end -@inline function _eig{T<:Real}(s::Size, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eig(s::Size, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} eigen = eigfact(Hermitian(Array(parent(A))); permute=permute, scale=scale) return (s(eigen.values), s(eigen.vectors)) # Return a SizedArray end -@inline function _eig{T<:Real}(::Size{(1,1)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eig(::Size{(1,1)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} @inbounds return (SVector{1,T}((A[1],)), eye(SMatrix{1,1,T})) end # TODO adapt the below to be complex-safe? -@inline function _eig{T<:Real}(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eig(::Size{(2,2)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} a = A.data if A.uplo == 'U' @@ -190,7 +193,7 @@ end # A small part of the code in the following method was inspired by works of David # Eberly, Geometric Tools LLC, in code released under the Boost Software # License (included at the end of this file). -@inline function _eig{T<:Real}(::Size{(3,3)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) +@inline function _eig(::Size{(3,3)}, A::Base.LinAlg.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real} S = typeof((one(T)*zero(T) + zero(T))/one(T)) Sreal = real(S) diff --git a/src/linalg.jl b/src/linalg.jl index f78dbcde..bae91705 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -341,6 +341,7 @@ end end # TODO same for `RowVector`? +@inline Size(::Union{RowVector{T, SA}, Type{RowVector{T, SA}}}) where {T, SA <: StaticArray} = Size(1, Size(SA)[1]) @inline Size(::Union{Symmetric{T,SA}, Type{Symmetric{T,SA}}}) where {T,SA<:StaticArray} = Size(SA) @inline Size(::Union{Hermitian{T,SA}, Type{Hermitian{T,SA}}}) where {T,SA<:StaticArray} = Size(SA) diff --git a/src/matrix_multiply.jl b/src/matrix_multiply.jl index fca8ae9b..547c940e 100644 --- a/src/matrix_multiply.jl +++ b/src/matrix_multiply.jl @@ -34,10 +34,12 @@ promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero @inline *(A::StaticMatrix, B::StaticVector) = _A_mul_B(Size(A), Size(B), A, B) @inline *(A::StaticMatrix, B::StaticMatrix) = _A_mul_B(Size(A), Size(B), A, B) @inline *(A::StaticVector, B::StaticMatrix) = *(reshape(A, Size(Size(A)[1], 1)), B) +@inline *(A::StaticVector, B::RowVector{<:Any, <:StaticVector}) = _A_mul_B(Size(A), Size(B), A, B) @inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticVector) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B) @inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticMatrix) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B) @inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::StaticMatrix) = A_mul_B!(dest, reshape(A, Size(Size(A)[1], 1)), B) +@inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::RowVector{<:Any, <:StaticVector}) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B) #@inline *{TA<:Base.LinAlg.BlasFloat,Tb}(A::StaticMatrix{TA}, b::StaticVector{Tb}) @@ -61,6 +63,18 @@ promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero end end +# outer product +@generated function _A_mul_B(::Size{sa}, ::Size{sb}, a::StaticVector{Ta}, b::RowVector{Tb, <:StaticVector}) where {sa, sb, Ta, Tb} + newsize = (sa[1], sb[2]) + exprs = [:(a[$i]*b[$j]) for i = 1:sa[1], j = 1:sb[2]] + + return quote + @_inline_meta + T = promote_op(*, Ta, Tb) + @inbounds return similar_type(b, T, Size($newsize))(tuple($(exprs...))) + end +end + # TODO: I removed StaticMatrix * AbstractVector. Reinstate? @generated function _A_mul_B(Sa::Size{sa}, Sb::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} @@ -219,29 +233,44 @@ end return quote @_inline_meta @inbounds $(Expr(:block, exprs...)) + return c + end +end + +@generated function _A_mul_B!(::Size{sc}, c::StaticMatrix, ::Size{sa}, ::Size{sb}, a::StaticVector, b::RowVector{<:Any, <:StaticVector}) where {sa, sb, sc} + if sa[1] != sc[1] || sb[2] != sc[2] + throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) + end + + exprs = [:(c[$(sub2ind(sc, i, j))] = a[$i] * b[$j]) for i = 1:sa[1], j = 1:sb[2]] + + return quote + @_inline_meta + @inbounds $(Expr(:block, exprs...)) + return c end end -@generated function _A_mul_B!(::Size{sc}, c::StaticMatrix{Tc}, ::Size{sa}, ::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, sc, Ta, Tb, Tc} +@generated function _A_mul_B!(Sc::Size{sc}, c::StaticMatrix{Tc}, Sa::Size{sa}, Sb::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, sc, Ta, Tb, Tc} can_blas = Tc == Ta && Tc == Tb && Tc <: BlasFloat if can_blas if sa[1] * sa[2] * sb[2] < 4*4*4 return quote @_inline_meta - A_mul_B_unrolled!(c, a, b) + A_mul_B_unrolled!(Sc, c, Sa, Sb, a, b) return c end elseif sa[1] * sa[2] * sb[2] < 14*14*14 # Something seems broken for this one with large matrices (becomes allocating) return quote @_inline_meta - A_mul_B_unrolled_chunks!(c, a, b) + A_mul_B_unrolled_chunks!(Sc, c, Sa, Sb, a, b) return c end else return quote @_inline_meta - A_mul_B_blas!(c, a, b) + A_mul_B_blas!(Sc, c, Sa, Sb, a, b) return c end end @@ -249,13 +278,13 @@ end if sa[1] * sa[2] * sb[2] < 4*4*4 return quote @_inline_meta - A_mul_B_unrolled!(c, a, b) + A_mul_B_unrolled!(Sc, c, Sa, Sb, a, b) return c end else return quote @_inline_meta - A_mul_B_unrolled_chunks!(c, a, b) + A_mul_B_unrolled_chunks!(Sc, c, Sa, Sb, a, b) return c end end @@ -312,7 +341,7 @@ end @generated function A_mul_B_unrolled!(::Size{sc}, c::StaticMatrix, ::Size{sa}, ::Size{sb}, a::StaticMatrix, b::StaticMatrix) where {sa, sb, sc} - if sb[1] != sa[2] || sa[1] != s[1] || sb[2] != sc[2] + if sb[1] != sa[2] || sa[1] != sc[1] || sb[2] != sc[2] throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) end @@ -329,7 +358,7 @@ end end @generated function A_mul_B_unrolled_chunks!(::Size{sc}, c::StaticMatrix, ::Size{sa}, ::Size{sb}, a::StaticMatrix, b::StaticMatrix) where {sa, sb, sc} - if sb[1] != sa[2] || sa[1] != s[1] || sb[2] != sc[2] + if sb[1] != sa[2] || sa[1] != sc[1] || sb[2] != sc[2] throw(DimensionMismatch("Tried to multiply arrays of size $sa and $sb and assign to array of size $sc")) end @@ -337,9 +366,9 @@ end # Do a custom b[:, k2] to return a SVector (an isbits type) rather than a mutable type. Avoids allocation == faster tmp_type = SVector{sb[1], eltype(c)} - vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(a, $(Expr(:call, tmp_type, [Expr(:ref, :b, sub2ind(s, i, k2)) for i = 1:sb[1]]...)))) for k2 = 1:sb[2]] + vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply($(Size(sa)), $(Size(sb[1])), a, $(Expr(:call, tmp_type, [Expr(:ref, :b, sub2ind(sc, i, k2)) for i = 1:sb[1]]...)))) for k2 = 1:sb[2]] - exprs = [:(c[$(sub2ind(s, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sa[1], k2 = 1:sb[2]] + exprs = [:(c[$(sub2ind(sc, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sa[1], k2 = 1:sb[2]] return quote @_inline_meta diff --git a/src/solve.jl b/src/solve.jl index ddc63d1e..155a1569 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -1,40 +1,28 @@ -@generated function (\){T}(A::StaticMatrix{T}, b::StaticVector{T}) - S = typeof((one(T)*zero(T) + zero(T))/one(T)) - newtype = similar_type(b, S) +@inline (\)(a::StaticMatrix{T}, b::StaticVector{T}) where {T} = solve(Size(a), Size(b), a, b) - if size(A) == (1,1) && length(b) == 1 - return quote - $(Expr(:meta, :inline)) - @inbounds return $newtype( b[1] / A[1,1] ) - end - elseif size(A) == (2,2) && length(b) == 2 - return quote - $(Expr(:meta, :inline)) - d = det(A) - @inbounds return $newtype( - A[2,2]*b[1] - A[1,2]*b[2], - A[1,1]*b[2] - A[2,1]*b[1] ) / d - end - elseif size(A) == (3,3) && length(b) == 3 - return quote - $(Expr(:meta, :inline)) - d = det(A) - @inbounds return $newtype( - (A[2,2]*A[3,3] - A[2,3]*A[3,2])*b[1] + - (A[1,3]*A[3,2] - A[1,2]*A[3,3])*b[2] + - (A[1,2]*A[2,3] - A[1,3]*A[2,2])*b[3], - (A[2,3]*A[3,1] - A[2,1]*A[3,3])*b[1] + - (A[1,1]*A[3,3] - A[1,3]*A[3,1])*b[2] + - (A[1,3]*A[2,1] - A[1,1]*A[2,3])*b[3], - (A[2,1]*A[3,2] - A[2,2]*A[3,1])*b[1] + - (A[1,2]*A[3,1] - A[1,1]*A[3,2])*b[2] + - (A[1,1]*A[2,2] - A[1,2]*A[2,1])*b[3] ) / d - end - else - # FixMe! Unsatisfactory ineffective but requires some infrastructure - # to make efficient so we fall back on inv for now - quote - inv(A)*b - end - end -end \ No newline at end of file +# TODO: Ineffective but requires some infrastructure (e.g. LU or QR) to make efficient so we fall back on inv for now +@inline solve(::Size, ::Size, a, b) = inv(a) * b + +@inline solve(::Size{(1,1)}, ::Size{(1,)}, a, b) = similar_type(b, typeof(b[1] \ a[1]))(b[1] \ a[1]) + +@inline function solve(::Size{(2,2)}, ::Size{(2,)}, a::StaticMatrix{Ta}, b::StaticVector{Tb}) where {Ta, Tb} + d = det(a) + T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/d) + @inbounds return similar_type(b, T)((a[2,2]*b[1] - a[1,2]*b[2])/d, + (a[1,1]*b[2] - a[2,1]*b[1])/d) +end + +@inline function solve(::Size{(3,3)}, ::Size{(3,)}, a::StaticMatrix{Ta}, b::StaticVector{Tb}) where {Ta, Tb} + d = det(a) + T = typeof((one(Ta)*zero(Tb) + one(Ta)*zero(Tb))/d) + @inbounds return similar_type(b, T)( + ((a[2,2]*a[3,3] - a[2,3]*a[3,2])*b[1] + + (a[1,3]*a[3,2] - a[1,2]*a[3,3])*b[2] + + (a[1,2]*a[2,3] - a[1,3]*a[2,2])*b[3]) / d, + ((a[2,3]*a[3,1] - a[2,1]*a[3,3])*b[1] + + (a[1,1]*a[3,3] - a[1,3]*a[3,1])*b[2] + + (a[1,3]*a[2,1] - a[1,1]*a[2,3])*b[3]) / d, + ((a[2,1]*a[3,2] - a[2,2]*a[3,1])*b[1] + + (a[1,2]*a[3,1] - a[1,1]*a[3,2])*b[2] + + (a[1,1]*a[2,2] - a[1,2]*a[2,1])*b[3]) / d ) +end diff --git a/test/SVector.jl b/test/SVector.jl index 6b33380c..0258a7af 100644 --- a/test/SVector.jl +++ b/test/SVector.jl @@ -17,6 +17,7 @@ @test SVector((1.0,)).data === (1.0,) @test SVector(1).data === (1,) + @test SVector(1,1.0).data === (1.0,1.0) @test ((@SVector [1.0])::SVector{1}).data === (1.0,) @test ((@SVector [1, 2, 3])::SVector{3}).data === (1, 2, 3) diff --git a/test/abstractarray.jl b/test/abstractarray.jl index 1e788b3d..c77447fd 100644 --- a/test/abstractarray.jl +++ b/test/abstractarray.jl @@ -63,5 +63,7 @@ @testset "reshape" begin @test @inferred(reshape(SVector(1,2,3,4), Size(2,2))) === SMatrix{2,2}(1,2,3,4) @test @inferred(reshape([1,2,3,4], Size(2,2)))::SizedArray{(2,2),Int,2,1} == [1 3; 2 4] + + @test @inferred(vec(SMatrix{2, 2}([1 2; 3 4])))::SVector{4,Int} == [1, 3, 2, 4] end end diff --git a/test/matrix_multiply.jl b/test/matrix_multiply.jl index ecae737a..cda5c033 100644 --- a/test/matrix_multiply.jl +++ b/test/matrix_multiply.jl @@ -9,6 +9,12 @@ @test isa(x, SVector{2,CartesianIndex{2}}) @test x == @SVector [CartesianIndex((7,5)), CartesianIndex((15,13))] + # inner product + @test @inferred(v'*v) === 5 + + # outer product + @test @inferred(v*v') === @SMatrix [1 2; 2 4] + v3 = [1, 2] @test_broken m*v3 === @SVector [5, 11] @@ -128,7 +134,6 @@ @test m*n::MMatrix ≈ a_array end -#= @testset "A_mul_B!" begin v = @SVector [2, 4] v2 = [2, 4] @@ -216,5 +221,5 @@ outvec2f = Vector{Float64}(2) A_mul_B!(outvec2f, mf, vf2) @test outvec2f ≈ [10.0, 22.0] - end =# + end end diff --git a/test/runtests.jl b/test/runtests.jl index 98b64fe8..34c28cdf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,11 +17,11 @@ using Base.Test include("mapreduce.jl") include("arraymath.jl") include("linalg.jl") - include("matrix_multiply.jl") #= + include("matrix_multiply.jl") include("det.jl") include("inv.jl") include("solve.jl") include("eigen.jl") include("deque.jl") - #include("fixed_size_arrays.jl") =# + #include("fixed_size_arrays.jl") end From ed6fb3b6596d9ebf6129a15202655fe7a4fcc677 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Tue, 21 Mar 2017 14:50:03 +1000 Subject: [PATCH 5/9] Updated travis, appveyer --- .travis.yml | 8 ++++---- appveyor.yml | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.travis.yml b/.travis.yml index 6fd84510..4247baed 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,14 +4,14 @@ os: - linux - osx julia: - - 0.5 + - 0.6 - nightly notifications: email: false -matrix: - allow_failures: +# matrix: + # allow_failures: # - julia: nightly - - os: osx + # - os: osx # uncomment the following lines to override the default test script #script: # - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi diff --git a/appveyor.yml b/appveyor.yml index b1a56438..4f74dab0 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,7 +1,7 @@ environment: matrix: - - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe" - - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe" +# - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe" +# - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe" - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe" - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe" From 0f3494b08fdcf15dd9c61dbf75be2605585f2da2 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Tue, 21 Mar 2017 14:58:26 +1000 Subject: [PATCH 6/9] Enable code coverage of v0.6 --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 4247baed..843bb4e7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -17,6 +17,6 @@ notifications: # - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi # - julia -e 'Pkg.clone(pwd()); Pkg.build("StaticArrays"); Pkg.test("StaticArrays"; coverage=true)' after_success: - - if [ $TRAVIS_JULIA_VERSION = "0.5" ] && [ $TRAVIS_OS_NAME = "linux" ]; then + - if [ $TRAVIS_JULIA_VERSION = "0.6" ] && [ $TRAVIS_OS_NAME = "linux" ]; then julia -e 'cd(Pkg.dir("StaticArrays")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'; fi From 003c40fc23cce43e4e206a1717ac355e2bc9eb4c Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Wed, 22 Mar 2017 17:25:53 +1000 Subject: [PATCH 7/9] Code review fixes --- src/MMatrix.jl | 19 ++++++-------- src/SMatrix.jl | 38 ++++++---------------------- src/abstractarray.jl | 5 ++-- src/arraymath.jl | 6 +---- src/broadcast.jl | 4 +-- src/convert.jl | 60 -------------------------------------------- src/det.jl | 5 ++-- test/MMatrix.jl | 4 +-- test/SMatrix.jl | 6 ++--- 9 files changed, 28 insertions(+), 119 deletions(-) diff --git a/src/MMatrix.jl b/src/MMatrix.jl index 6e34f0ae..0125bc53 100644 --- a/src/MMatrix.jl +++ b/src/MMatrix.jl @@ -34,27 +34,26 @@ type MMatrix{S1, S2, T, L} <: StaticMatrix{T} end end -@generated function check_MMatrix_params{S1,S2,L}(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}}) - if !(T <: DataType) # I think the way types are handled in generated fnctions might have changed in 0.5? - return :(error("MMatrix: Parameter T must be a DataType. Got $T")) - end +function check_MMatrix_params(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}}) where {S1,S2,L} + throw(ArgumentError("MMatrix: Parameter T must be a Type. Got $T")) +end +@generated function check_MMatrix_params(::Type{Val{S1}}, ::Type{Val{S2}}, ::Type{T}, ::Type{Val{L}}) where {S1,S2,L,T} if !isa(S1, Int) || !isa(S2, Int) || !isa(L, Int) || S1 < 0 || S2 < 0 || L < 0 - return :(error("MMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)")) + throw(ArgumentError("MMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)")) end if S1*S2 == L return nothing else - str = "Size mismatch in MMatrix. S1 = $S1, S2 = $S2, but recieved $L elements" - return :(error($str)) + throw(ArgumentError("Size mismatch in MMatrix. S1 = $S1, S2 = $S2, but recieved $L elements")) end end @generated function (::Type{MMatrix{S1}}){S1,L}(x::NTuple{L}) S2 = div(L, S1) if S1*S2 != L - error("Incorrect matrix sizes. $S1 does not divide $L elements") + throw(DimensionMismatch("Incorrect matrix sizes. $S1 does not divide $L elements")) end T = promote_tuple_eltype(x) @@ -90,13 +89,9 @@ end @inline convert{S1,S2,T}(::Type{MMatrix{S1,S2}}, a::StaticArray{T}) = MMatrix{S1,S2,T}(Tuple(a)) @inline MMatrix(a::StaticMatrix) = MMatrix{size(typeof(a),1),size(typeof(a),2)}(Tuple(a)) - # Some more advanced constructor-like functions @inline one{N}(::Type{MMatrix{N}}) = one(MMatrix{N,N}) @inline eye{N}(::Type{MMatrix{N}}) = eye(MMatrix{N,N}) -#@inline eye{N,M}(::Type{MMatrix{N,M}}) = eye(MMatrix{N,M,Float64}) -#@inline zeros{N,M}(::Type{MMatrix{N,M}}) = zeros(MMatrix{N,M,Float64}) -#@inline ones{N,M}(::Type{MMatrix{N,M}}) = ones(MMatrix{N,M,Float64}) ##################### ## MMatrix methods ## diff --git a/src/SMatrix.jl b/src/SMatrix.jl index 3c64d74f..af78162d 100644 --- a/src/SMatrix.jl +++ b/src/SMatrix.jl @@ -28,27 +28,26 @@ immutable SMatrix{S1, S2, T, L} <: StaticMatrix{T} end end -@generated function check_smatrix_params{S1,S2,L}(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}}) - if !(T <: DataType) # I think the way types are handled in generated fnctions might have changed in 0.5? - return :(error("SMatrix: Parameter T must be a DataType. Got $T")) - end +function check_smatrix_params(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}}) where {S1,S2,L} + throw(ArgumentError("SMatrix: Parameter T must be a Type. Got $T")) +end +@generated function check_smatrix_params(::Type{Val{S1}}, ::Type{Val{S2}}, ::Type{T}, ::Type{Val{L}}) where {S1,S2,L,T} if !isa(S1, Int) || !isa(S2, Int) || !isa(L, Int) || S1 < 0 || S2 < 0 || L < 0 - return :(error("SMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)")) + throw(ArgumentError("SMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)")) end if S1*S2 == L return nothing else - str = "Size mismatch in SMatrix. S1 = $S1, S2 = $S2, but recieved $L elements" - return :(error($str)) + throw(ArgumentError("Size mismatch in SMatrix. S1 = $S1, S2 = $S2, but recieved $L elements")) end end @generated function (::Type{SMatrix{S1}}){S1,L}(x::NTuple{L,Any}) S2 = div(L, S1) if S1*S2 != L - error("Incorrect matrix sizes. $S1 does not divide $L elements") + throw(DimensionMismatch("Incorrect matrix sizes. $S1 does not divide $L elements")) end T = promote_tuple_eltype(x) @@ -85,32 +84,9 @@ end @inline convert{S1,S2,T}(::Type{SMatrix{S1,S2}}, a::StaticArray{T}) = SMatrix{S1,S2,T}(Tuple(a)) @inline SMatrix(a::StaticMatrix) = SMatrix{size(typeof(a),1),size(typeof(a),2)}(Tuple(a)) -#= -@inline (::Type{SMatrix{S1}}){S1}(x1) = SMatrix{S1}((x1,)) -@inline (::Type{SMatrix{S1}}){S1}(x1,x2) = SMatrix{S1}((x1,x2)) -@inline (::Type{SMatrix{S1}}){S1}(x1,x2,x3) = SMatrix{S1}((x1,x2,x3)) -@inline (::Type{SMatrix{S1}}){S1}(x1,x2,x3,x4) = SMatrix{S1}((x1,x2,x3,x4)) -@inline (::Type{SMatrix{S1}}){S1}(x...) = SMatrix{S1}(x) - -@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1) = SMatrix{S1,S2}((x1,)) -@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2) = SMatrix{S1,S2}((x1,x2)) -@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2,x3) = SMatrix{S1,S2}((x1,x2,x3)) -@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2,x3,x4) = SMatrix{S1,S2}((x1,x2,x3,x4)) -@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x...) = SMatrix{S1,S2}(x) - -@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1) = SMatrix{S1,S2,T}((x1,)) -@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2) = SMatrix{S1,S2,T}((x1,x2)) -@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2,x3) = SMatrix{S1,S2,T}((x1,x2,x3)) -@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2,x3,x4) = SMatrix{S1,S2,T}((x1,x2,x3,x4)) -@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x...) = SMatrix{S1,S2,T}(x) -=# - # Some more advanced constructor-like functions @inline one{N}(::Type{SMatrix{N}}) = one(SMatrix{N,N}) @inline eye{N}(::Type{SMatrix{N}}) = eye(SMatrix{N,N}) -#@inline eye{N,M}(::Type{SMatrix{N,M}}) = eye(SMatrix{N,M,Float64}) -#@inline zeros{N,M}(::Type{SMatrix{N,M}}) = zeros(SMatrix{N,M,Float64}) -#@inline ones{N,M}(::Type{SMatrix{N,M}}) = ones(SMatrix{N,M,Float64}) ##################### ## SMatrix methods ## diff --git a/src/abstractarray.jl b/src/abstractarray.jl index ec7c6f9e..129ac3ff 100644 --- a/src/abstractarray.jl +++ b/src/abstractarray.jl @@ -105,8 +105,9 @@ similar{A<:Array,T,S}(::Type{A},::Type{T},s::Size{S}) = sizedarray_similar_type( @inline reshape(a::StaticArray, s::Size) = similar_type(a, s)(Tuple(a)) -@generated function reshape(a::AbstractArray, s::Size{S}) where {S} - if IndexStyle(a) == IndexLinear() # TODO this isn't "hyperpure" +@inline reshape(a::AbstractArray, s::Size) = _reshape(s, IndexStyle(a), a) +@generated function _reshape(a::AbstractArray, indexstyle, s::Size{S}) where {S} + if indexstyle == IndexLinear exprs = [:(a[$i]) for i = 1:prod(S)] else exprs = [:(a[$(inds...)]) for inds ∈ CartesianRange(S)] diff --git a/src/arraymath.jl b/src/arraymath.jl index 93d28ec9..1d308fcb 100644 --- a/src/arraymath.jl +++ b/src/arraymath.jl @@ -1,7 +1,3 @@ -# Support for elementwise ops on AbstractArray{S<:StaticArray} with Number -#Base.promote_op{Op,A<:StaticArray,T<:Number}(op::Op, ::Type{A}, ::Type{T}) = similar_type(A, promote_op(op, eltype(A), T)) -#Base.promote_op{Op,T<:Number,A<:StaticArray}(op::Op, ::Type{T}, ::Type{A}) = similar_type(A, promote_op(op, T, eltype(A))) - @inline zeros(::SA) where {SA <: StaticArray} = zeros(SA) @generated function zeros(::Type{SA}) where {SA <: StaticArray} T = eltype(SA) @@ -34,7 +30,7 @@ end end end -@inline fill(val, ::SA) where {SA <: StaticArray} = ones(val, SA) +@inline fill(val, ::SA) where {SA <: StaticArray} = _fill(val, Size(SA), SA) @inline fill(val, ::Type{SA}) where {SA <: StaticArray} = _fill(val, Size(SA), SA) @generated function _fill(val, ::Size{s}, ::Type{SA}) where {s, SA <: StaticArray} v = [:val for i = 1:prod(s)] diff --git a/src/broadcast.jl b/src/broadcast.jl index 22696818..d73dde7e 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -60,7 +60,7 @@ end exprs_vals = [(a[i] <: Number ? :(a[$i]) : :(a[$i][$(broadcasted_index(sizes[i], current_ind))])) for i = 1:length(sizes)] exprs[current_ind...] = :(f($(exprs_vals...))) - # increment current_ind + # increment current_ind (maybe use CartesianRange?) current_ind[1] += 1 for i ∈ 1:length(newsize) if current_ind[i] > newsize[i] @@ -122,7 +122,7 @@ end exprs_vals = [(a[i] <: Number ? :(a[$i]) : :(a[$i][$(broadcasted_index(sizes[i], current_ind))])) for i = 1:length(sizes)] exprs[current_ind...] = :(dest[$j] = f($(exprs_vals...))) - # increment current_ind + # increment current_ind (maybe use CartesianRange?) current_ind[1] += 1 for i ∈ 1:length(newsize) if current_ind[i] > newsize[i] diff --git a/src/convert.jl b/src/convert.jl index d799aa6c..57778f92 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -30,63 +30,3 @@ length_val(a::Type{T}) where {T<:StaticArray} = length_val(Size(T)) @inbounds return $(Expr(:tuple, exprs...)) end end - -#= -@generated function convert(::Type{Tuple}, a::StaticArray) - n = length(a) - exprs = [:(a[$j]) for j = 1:n] - quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:tuple, exprs...)) - end -end -=# - - - -# People might not want to use Tuple for everything (TODO: check this with FieldVector...) -# Generic case, with least 2 inputs -#@inline (::Type{SA}){SA<:StaticArray}(x1,x2,xs...) = SA((x1,x2,xs...)) - - -#= -function convert{T,N}(::Type{Array}, sa::StaticArray{T,N}) - out = Array{T,N}(size(sa)) - @inbounds for i = 1:length(sa) - out[i] = sa[i] - end - return out -end - -function convert{T,N}(::Type{Array{T}}, sa::StaticArray{T,N}) - out = Array{T,N}(size(sa)) - @inbounds for i = 1:length(sa) - out[i] = sa[i] - end - return out -end - -function convert{T,N}(::Type{Array{T,N}}, sa::StaticArray{T,N}) - out = Array{T,N}(size(sa)) - @inbounds for i = 1:length(sa) - out[i] = sa[i] - end - return out -end - -function convert{T}(::Type{Matrix}, sa::StaticMatrix{T}) - out = Matrix{T}(size(sa)) - @inbounds for i = 1:length(sa) - out[i] = sa[i] - end - return out -end - -function convert{T}(::Type{Vector}, sa::StaticVector{T}) - out = Vector{T}(size(sa)) - @inbounds for i = 1:length(sa) - out[i] = sa[i] - end - return out -end -=# diff --git a/src/det.jl b/src/det.jl index c3164781..18e698cc 100644 --- a/src/det.jl +++ b/src/det.jl @@ -6,7 +6,7 @@ @inbounds return A[1]*A[4] - A[3]*A[2] end -@inline function _det(::Size{(2,2)}, A::StaticMatrix{T}) where {T <: Unsigned} +@inline function _det(::Size{(2,2)}, A::StaticMatrix{<:Unsigned}) @inbounds return Signed(A[1]*A[4]) - Signed(A[3]*A[2]) end @@ -17,7 +17,7 @@ end return vecdot(x0, cross(x1, x2)) end -@inline function _det{T<:Unsigned}(::Size{(3,3)}, A::StaticMatrix{T}) +@inline function _det(::Size{(3,3)}, A::StaticMatrix{<:Unsigned}) @inbounds x0 = SVector(Signed(A[1]), Signed(A[2]), Signed(A[3])) @inbounds x1 = SVector(Signed(A[4]), Signed(A[5]), Signed(A[6])) @inbounds x2 = SVector(Signed(A[7]), Signed(A[8]), Signed(A[9])) @@ -29,6 +29,7 @@ end throw(DimensionMismatch("matrix is not square")) end return quote # Implementation from Base + @_inline_meta T2 = typeof((one(T)*zero(T) + zero(T))/one(T)) if istriu(A) || istril(A) return convert(T2, det(UpperTriangular(A))) # Is this a Julia bug that a convert is not type stable?? diff --git a/test/MMatrix.jl b/test/MMatrix.jl index 88af0af9..f7178de8 100644 --- a/test/MMatrix.jl +++ b/test/MMatrix.jl @@ -13,8 +13,8 @@ # Bad parameters @test_throws Exception MMatrix{1,1,Int,2}((1,)) @test_throws Exception MMatrix{1,1,1,1}((1,)) - @test_throws Exception MMatrix{1,2,Int,1}((1,)) - @test_throws Exception MMatrix{2,1,Int,1}((1,)) + @test_throws ArgumentError MMatrix{1,2,Int,1}((1,)) + @test_throws ArgumentError MMatrix{2,1,Int,1}((1,)) end @testset "Outer constructors and macro" begin diff --git a/test/SMatrix.jl b/test/SMatrix.jl index 8e41a171..50ad8b1d 100644 --- a/test/SMatrix.jl +++ b/test/SMatrix.jl @@ -11,9 +11,9 @@ # Bad parameters @test_throws Exception SMatrix{1,1,Int,2}((1,)) - @test_throws Exception SMatrix{1,1,1,1}((1,)) - @test_throws Exception SMatrix{1,2,Int,1}((1,)) - @test_throws Exception SMatrix{2,1,Int,1}((1,)) + @test_throws ArgumentError SMatrix{1,1,1,1}((1,)) + @test_throws ArgumentError SMatrix{1,2,Int,1}((1,)) + @test_throws ArgumentError SMatrix{2,1,Int,1}((1,)) end @testset "Outer constructors and macro" begin From a0d6658d89e1dd72927bb632d3c90aa436f1ebfe Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Fri, 24 Mar 2017 09:32:14 +1000 Subject: [PATCH 8/9] Fixed #123, inner constructor ambiguous with error message constructor --- src/convert.jl | 2 +- test/custom_types.jl | 7 +++++++ test/runtests.jl | 1 + 3 files changed, 9 insertions(+), 1 deletion(-) create mode 100644 test/custom_types.jl diff --git a/src/convert.jl b/src/convert.jl index 57778f92..2eab77c1 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -1,4 +1,4 @@ -(::Type{SA})(x::Tuple) where {SA <: StaticArray} = error("No precise constructor for $SA found. Length of input was $(length(x)).") +(::Type{SA})(x::Tuple{Tuple{Tuple{<:Tuple}}}) where {SA <: StaticArray} = error("No precise constructor for $SA found. Length of input was $(length(x[1][1][1])).") @inline (::Type{SA})(x...) where {SA <: StaticArray} = SA(x) @inline (::Type{SA})(a::AbstractArray) where {SA <: StaticArray} = convert(SA, a) # Is this a good idea? diff --git a/test/custom_types.jl b/test/custom_types.jl new file mode 100644 index 00000000..251779a5 --- /dev/null +++ b/test/custom_types.jl @@ -0,0 +1,7 @@ +@testset "Custom types" begin + # Issue 123 + @eval (struct MyType{N, T} <: StaticVector{T} + data::NTuple{N, T} + end) + @test (MyType(3, 4) isa MyType{2, Int}) +end diff --git a/test/runtests.jl b/test/runtests.jl index 34c28cdf..addd19c3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,7 @@ using Base.Test include("MArray.jl") include("FieldVector.jl") include("Scalar.jl") + include("custom_types.jl") include("core.jl") include("abstractarray.jl") From 5ebaf5cf6c06f4e394c4350193836babe056893f Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Fri, 24 Mar 2017 17:01:34 +1000 Subject: [PATCH 9/9] More code review fixes --- src/convert.jl | 1 + src/linalg.jl | 51 ++--------------------- src/matrix_multiply.jl | 91 ++++++++++++++++++++---------------------- src/util.jl | 46 --------------------- test/custom_types.jl | 6 +++ test/linalg.jl | 1 + 6 files changed, 55 insertions(+), 141 deletions(-) diff --git a/src/convert.jl b/src/convert.jl index 2eab77c1..d5f50a4c 100644 --- a/src/convert.jl +++ b/src/convert.jl @@ -6,6 +6,7 @@ # this covers most conversions and "statically-sized reshapes" @inline convert(::Type{SA}, sa::StaticArray) where {SA<:StaticArray} = SA(Tuple(sa)) @inline convert(::Type{SA}, sa::SA) where {SA<:StaticArray} = sa +@inline convert(::Type{SA}, x::Tuple) where {SA<:StaticArray} = SA(x) # convert -> constructor. Hopefully no loops... # A general way of going back to a tuple @inline function convert(::Type{Tuple}, a::StaticArray) diff --git a/src/linalg.jl b/src/linalg.jl index bae91705..3ce0ff73 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -54,34 +54,9 @@ end # Transpose, conjugate, etc -# TODO different methods for v0.5, v0.6 (due to `RowVector`) @inline conj(a::StaticArray) = map(conj, a) - -#= -@generated function transpose(v::StaticVector) - n = length(v) - newtype = similar_type(v, Size(1,n)) - exprs = [:(v[$j]) for j = 1:n] - - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end - -@generated function ctranspose(v::StaticVector) - n = length(v) - newtype = similar_type(v, Size(1,n)) - exprs = [:(conj(v[$j])) for j = 1:n] - - return quote - $(Expr(:meta, :inline)) - @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))) - end -end -=# - @inline transpose(m::StaticMatrix) = _transpose(Size(m), m) +# note: transpose of StaticVector is a RowVector, handled by Base @generated function _transpose(::Size{S}, m::StaticMatrix) where {S} Snew = (S[2], S[1]) @@ -214,34 +189,16 @@ _cross(::Size{S}, a::StaticVector, b::StaticVector) where {S} = error("Cross pro @inbounds return a[1]*b[2] - a[2]*b[1] end @inline function _cross(::Size{(3,)}, a::StaticVector, b::StaticVector) - T = typeof(a[2]*b[3]-a[3]*b[2]) @inbounds return similar_type(a, typeof(a[2]*b[3]-a[3]*b[2]))((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 @inline function _cross(::Size{(2,)}, a::StaticVector{<:Unsigned}, b::StaticVector{<:Unsigned}) @inbounds return Signed(a[1]*b[2]) - Signed(a[2]*b[1]) end @inline function _cross(::Size{(3,)}, a::StaticVector{<:Unsigned}, b::StaticVector{<:Unsigned}) - T = typeof(a[2]*b[3]-a[3]*b[2]) @inbounds return similar_type(a, typeof(Signed(a[2]*b[3])-Signed(a[3]*b[2])))(((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 -@inline dot(a::StaticVector, b::StaticVector) = _dot(same_size(a, b), a, b) -@generated function _dot(::Size{S}, a::StaticVector, b::StaticVector) where {S} - if S[1] == 0 - return :(zero(promote_op(*, eltype(a), eltype(b)))) - end - - expr = :(conj(a[1]) * b[1]) - for j = 2:S[1] - expr = :($expr + conj(a[$j]) * b[$j]) - end - - return quote - @_inline_meta - @inbounds return $expr - end -end - +@inline dot(a::StaticVector, b::StaticVector) = _vecdot(same_size(a, b), a, b) @inline vecdot(a::StaticArray, b::StaticArray) = _vecdot(same_size(a, b), a, b) @generated function _vecdot(::Size{S}, a::StaticArray, b::StaticArray) where {S} if prod(S) == 0 @@ -270,9 +227,9 @@ end return zero(real(eltype(a))) end - expr = :(real(conj(a[1]) * a[1])) + expr = :(abs2(a[1])) for j = 2:prod(S) - expr = :($expr + real(conj(a[$j]) * a[$j])) + expr = :($expr + abs2(a[$j])) end return quote diff --git a/src/matrix_multiply.jl b/src/matrix_multiply.jl index 547c940e..bc9c3fb3 100644 --- a/src/matrix_multiply.jl +++ b/src/matrix_multiply.jl @@ -63,6 +63,8 @@ promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero end end +# TODO: I removed StaticMatrix * AbstractVector. Reinstate? + # outer product @generated function _A_mul_B(::Size{sa}, ::Size{sb}, a::StaticVector{Ta}, b::RowVector{Tb, <:StaticVector}) where {sa, sb, Ta, Tb} newsize = (sa[1], sb[2]) @@ -75,58 +77,51 @@ end end end -# TODO: I removed StaticMatrix * AbstractVector. Reinstate? - @generated function _A_mul_B(Sa::Size{sa}, Sb::Size{sb}, a::StaticMatrix{Ta}, b::StaticMatrix{Tb}) where {sa, sb, Ta, Tb} - can_mutate = a.mutable && b.mutable # TODO this probably isn't safe. Maybe a trait?? - can_blas = Ta == Tb && Ta <: BlasFloat + # Heuristic choice for amount of codegen + if sa[1]*sa[2]*sb[2] <= 8*8*8 + return quote + @_inline_meta + return A_mul_B_unrolled(Sa, Sb, a, b) + end + elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 + return quote + @_inline_meta + return A_mul_B_unrolled_chunks(Sa, Sb, a, b) + end + else + return quote + @_inline_meta + return A_mul_B_loop(Sa, Sb, a, b) + end + end +end - if can_mutate - S = Size(sa[1], sb[2]) +@generated function _A_mul_B(Sa::Size{sa}, Sb::Size{sb}, a::Union{SizedMatrix{T}, MMatrix{T}, MArray{T}}, b::Union{SizedMatrix{T}, MMatrix{T}, MArray{T}}) where {sa, sb, T <: BlasFloat} + S = Size(sa[1], sb[2]) - # Heuristic choice between BLAS and explicit unrolling (or chunk-based unrolling) - if can_blas && sa[1]*sa[2]*sb[2] >= 14*14*14 - return quote - @_inline_meta - T = promote_matprod(Ta, Tb) - C = similar(a, T, $S) - A_mul_B_blas!($S, C, Sa, Sb, a, b) - return C - end - elseif sa[1]*sa[2]*sb[2] < 8*8*8 - return quote - @_inline_meta - return A_mul_B_unrolled(Sa, Sb, a, b) - end - elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 - return quote - @_inline_meta - T = promote_matprod(Ta, Tb) - return similar_type(a, T, $S)(A_mul_B_unrolled_chunks(Sa, Sb, a, b)) - end - else - return quote - @_inline_meta - return A_mul_B_loop(Sa, Sb, a, b) - end + # Heuristic choice between BLAS and explicit unrolling (or chunk-based unrolling) + if sa[1]*sa[2]*sb[2] >= 14*14*14 + return quote + @_inline_meta + C = similar(a, T, $S) + A_mul_B_blas!($S, C, Sa, Sb, a, b) + return C end - else # both are isbits type... - # Heuristic choice for amount of codegen - if sa[1]*sa[2]*sb[2] <= 8*8*8 - return quote - @_inline_meta - return A_mul_B_unrolled(Sa, Sb, a, b) - end - elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 - return quote - @_inline_meta - return A_mul_B_unrolled_chunks(Sa, Sb, a, b) - end - else - return quote - @_inline_meta - return A_mul_B_loop(Sa, Sb, a, b) - end + elseif sa[1]*sa[2]*sb[2] < 8*8*8 + return quote + @_inline_meta + return A_mul_B_unrolled(Sa, Sb, a, b) + end + elseif sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14 + return quote + @_inline_meta + return similar_type(a, T, $S)(A_mul_B_unrolled_chunks(Sa, Sb, a, b)) + end + else + return quote + @_inline_meta + return A_mul_B_loop(Sa, Sb, a, b) end end end diff --git a/src/util.jl b/src/util.jl index d8bf6a93..c2cbccda 100644 --- a/src/util.jl +++ b/src/util.jl @@ -27,49 +27,3 @@ end $t end end - - -# TODO: the below seems to be type piracy... -#= -# some convenience functions for non-static arrays, generators, etc... -@inline convert{T}(::Type{Tuple}, a::AbstractArray{T}) = (a...)::Tuple{Vararg{T}} -@inline function convert{N,T}(::Type{NTuple{N,Any}}, a::AbstractArray{T}) - @boundscheck if length(a) != N - error("Array of length $(length(a)) cannot be converted to a $N-tuple") - end - - @inbounds return ntuple(i -> a[i], Val{N}) -end - -@inline function convert{N,T1,T2}(::Type{NTuple{N,T1}}, a::AbstractArray{T2}) - @boundscheck if length(a) != N - error("Array of length $(length(a)) cannot be converted to a $N-tuple") - end - - @inbounds return ntuple(i -> convert(T1,a[i]), Val{N}) -end - - -if VERSION < v"0.5+" - # TODO try and make this generate fast code - @inline convert(::Type{Tuple}, g::Base.Generator) = (g...) - @inline function convert{N}(::Type{NTuple{N,Any}}, g::Base.Generator) - @boundscheck if length(g.iter) != N - error("Array of length $(length(a)) cannot be converted to a $N-tuple") - end - - @inbounds return ntuple(i -> g.f(g.iter[i]), Val{N}) - end -end -=# -#= -@generated function convert{N}(::Type{NTuple{N,Any}}, g::Base.Generator) - exprs = [:(g.f(g.iter[$j])) for j=1:N] - return quote - @boundscheck if length(g.iter) != N - error("Array of length $(length(a)) cannot be converted to a $N-tuple") - end - - @inbounds return $(Expr(:tuple, exprs...)) - end -end=# diff --git a/test/custom_types.jl b/test/custom_types.jl index 251779a5..dda6ffde 100644 --- a/test/custom_types.jl +++ b/test/custom_types.jl @@ -4,4 +4,10 @@ data::NTuple{N, T} end) @test (MyType(3, 4) isa MyType{2, Int}) + + # Issue 110 + @eval (struct Polly{N,T} + data::SVector{N,T} + end) + @test (Polly{2,Float64}((1.0, 0.0)) isa Polly) end diff --git a/test/linalg.jl b/test/linalg.jl index 4f568fad..1029b2a0 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -17,6 +17,7 @@ v3 = [2,4,6,8] v4 = [4,3,2,1] + # We broke "inferrable" sizes of AbstractVectors for vector+vector, matrix*vector, etc... @test_broken @inferred(v1 + v4) === @SVector [6, 7, 8, 9] @test_broken @inferred(v3 + v2) === @SVector [6, 7, 8, 9] @test_broken @inferred(v1 - v4) === @SVector [-2, 1, 4, 7]