Fix but in LowDimArray{ntuple(_ -> false, Val(N))}(data), and add vmapreduce to resolve #110.

Author: chriselrod

src/broadcast.jl

@@ -115,6 +115,7 @@ struct LowDimArray{D,T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
     data::A
 end
 @inline Base.pointer(A::LowDimArray) = pointer(A.data)
+Base.@propagate_inbounds Base.getindex(A::LowDimArray, i...) = getindex(A.data, i...)
 Base.size(A::LowDimArray) = Base.size(A.data)
 @generated function VectorizationBase.stridedpointer(A::LowDimArray{D,T,N}) where {D,T,N}
     s = Expr(:tuple, [Expr(:ref, :strideA, n) for n ∈ 1+D[1]:N if D[n]]...)
@@ -125,11 +126,19 @@ end
 function LowDimArray{D}(data::A) where {D,T,N,A <: AbstractArray{T,N}}
     LowDimArray{D,T,N,A}(data)
 end
+function extract_all_1_array!(ls::LoopSet, bcname::Symbol, N::Int, elementbytes::Int)
+    refextract = gensym(bcname)
+    pushpreamble!(ls, Expr(:(=), refextract, Expr(:ref, bcname, [1 for n ∈ 1:N]...)))
+    return add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
+end
 function add_broadcast!(
     ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
     @nospecialize(LDA::Type{<:LowDimArray}), elementbytes::Int
 )
     D,T,N::Int,_ = LDA.parameters
+    if !any(D)
+        return extract_all_1_array!(ls, bcname, N, elementbytes)
+    end
     fulldims = Symbol[loopsyms[n] for n ∈ 1:N if D[n]::Bool]
     ref = ArrayReference(bcname, fulldims)
     add_simple_load!(ls, destname, ref, elementbytes, true, false )::Operation
@@ -139,12 +148,14 @@ function add_broadcast_adjoint_array!(
 ) where {T,N,A<:AbstractArray{T,N}}
     parent = gensym(:parent)
     pushpreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))
+    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
     ref = ArrayReference(parent, Symbol[loopsyms[N + 1 - n] for n ∈ 1:N])
     add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
 end
 function add_broadcast_adjoint_array!(
     ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{<:AbstractVector}, elementbytes::Int
 )
+    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
     ref = ArrayReference(bcname, Symbol[loopsyms[2]])
     add_simple_load!( ls, destname, ref, elementbytes, true, true )
 end

src/mapreduce.jl

@@ -0,0 +1,76 @@
+
+@inline vreduce(::typeof(+), v::VectorizationBase.AbstractSIMDVector) = vsum(v)
+@inline vreduce(::typeof(*), v::VectorizationBase.AbstractSIMDVector) = vprod(v)
+@inline vreduce(::typeof(max), v::VectorizationBase.AbstractSIMDVector) = vmaximum(v)
+@inline vreduce(::typeof(min), v::VectorizationBase.AbstractSIMDVector) = vminimum(v)
+@inline vreduce(op, v::VectorizationBase.AbstractSIMDVector) = _vreduce(op, v)
+@inline _vreduce(op, v::VectorizationBase.AbstractSIMDVector) = _reduce(op, SVec(v))
+@inline function _vreduce(op, v::SVec)
+    isone(length(v)) && return v[1]
+    a = op(v[1], v[2])
+    for i ∈ 3:length(v)
+        a = op(a, v[i])
+    end
+    a
+end
+
+function mapreduce_simple(f::F, op::OP, args::Vararg{DenseArray{T},A}) where {F,OP,T<:NativeTypes,A}
+    ptrargs = ntuple(a -> pointer(args[a]), Val(A))
+    N = length(first(args))
+    iszero(N) && throw("Length of vector is 0!")
+    a_0 = f(vload.(ptrargs)...); i = 1
+    while i < N
+        a_0 = op(a_0, f(vload.(ptrargs, i)...)); i += 1
+    end
+    a_0
+end
+
+
+"""
+    vmapreduce(f, op, A::DenseArray...)
+
+Vectorized version of `mapreduce`. Applies `f` to each element of the arrays `A`, and reduces the result with `op`.
+"""
+function vmapreduce(f::F, op::OP, args::Vararg{DenseArray{T},A}) where {F,OP,T<:NativeTypes,A}
+    N = length(first(args))
+    A > 1 && @assert all(isequal(length.(args)...))
+    W = VectorizationBase.pick_vector_width(T)
+    V = VectorizationBase.pick_vector_width_val(T)
+    N < W && return mapreduce_simple(f, op, args...)
+    ptrargs = pointer.(args)
+        
+    a_0 = f(vload.(V, ptrargs)...); i = W
+    if N ≥ 4W
+        a_1 = f(vload.(V, ptrargs, i)...); i += W
+        a_2 = f(vload.(V, ptrargs, i)...); i += W
+        a_3 = f(vload.(V, ptrargs, i)...); i += W
+        while i < N - ((W << 2) - 1)
+            a_0 = op(a_0, f(vload.(V, ptrargs, i)...)); i += W
+            a_1 = op(a_1, f(vload.(V, ptrargs, i)...)); i += W
+            a_2 = op(a_2, f(vload.(V, ptrargs, i)...)); i += W
+            a_3 = op(a_3, f(vload.(V, ptrargs, i)...)); i += W
+        end
+        a_0 = op(a_0, a_1)
+        a_2 = op(a_2, a_3)
+        a_0 = op(a_0, a_2)
+    end
+    while i < N - (W - 1)
+        a_0 = op(a_0, f(vload.(V, ptrargs, i)...)); i += W
+    end
+    if i < N
+        m = mask(T, N & (W - 1))
+        a_0 = vifelse(m, op(a_0, f(vload.(V, ptrargs, i)...)), a_0)
+    end
+    vreduce(op, a_0)
+end
+
+@inline vmapreduce(f, op, args...) = mapreduce(f, op, args...)
+
+
+"""
+    vreduce(op, destination, A::DenseArray...)
+
+Vectorized version of `reduce`. Reduces the array `A` using the operator `op`.
+"""
+@inline vreduce(op, arg) = vmapreduce(identity, op, arg)
+

test/broadcast.jl

@@ -28,6 +28,11 @@
         fill!(c2, 99999);
         @avx @. c2 = a + bl;
         @test c1 ≈ c2
+
+        xs = rand(T, M);
+        max_ = maximum(xs, dims=1)
+        @test (@avx exp.(xs .- LowDimArray{(false,)}(max_))) ≈ exp.(xs .- LowDimArray{(false,)}(max_))
+
 
         a = rand(R, M); B = rand(R, M, N); c = rand(R, N); c′ = c';
         d1 =      @. a + B * c′;

test/gemm.jl

@@ -634,7 +634,7 @@
             @test C ≈ 2C2
             AmuladdBavx!(C, A, B, -1)
             @test C ≈ C2
-            AmuladdBavx!(C, At', B, -2)
+            AmuladdBavx!(C, At', Bt', -2)
             @test C ≈ -C2
             AmuladdBavx!(C, At', B, 3, 2)
             @test C ≈ C2

test/mapreduce.jl

@@ -0,0 +1,25 @@
+
+@testset "mapreduce" begin
+
+    for T ∈ (Int32, Int64, Float32, Float64)
+        if T <: Integer
+            R = T(1):T(100)
+            x7 = rand(R, 7); y7 = rand(R, 7);
+            x = rand(R, 127); y = rand(R, 127);
+        else
+            x7 = rand(T, 7); y7 = rand(T, 7);
+            x = rand(T, 127); y = rand(T, 127);
+            @test vmapreduce(hypot, +, x, y) ≈ mapreduce(hypot, +, x, y)
+            @test vmapreduce(^, (a,b) -> a + b, x7, y7) ≈ mapreduce(^, (a,b) -> a + b, x7, y7)
+        end
+        @test vreduce(+, x7) ≈ sum(x7)
+        @test vreduce(+, x) ≈ sum(x)
+        @test_throws AssertionError vmapreduce(hypot, +, x7, x)
+        @test vmapreduce(a -> 2a, *, x) ≈ mapreduce(a -> 2a, *, x)
+        @test vmapreduce(sin, +, x7) ≈ mapreduce(sin, +, x7)
+        @test vmapreduce(log, +, x) ≈ mapreduce(log, +, x)
+        @test vmapreduce(abs2, +, x) ≈ mapreduce(abs2, +, x)
+    end
+
+end
+