Fix for reductions while unrolling 5x-7x with static loops.

Author: chriselrod

Project.toml

@@ -1,7 +1,7 @@
 name = "LoopVectorization"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 authors = ["Chris Elrod <[email protected]>"]
-version = "0.8.19"
+version = "0.8.20"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -15,10 +15,10 @@ VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 [compat]
 DocStringExtensions = "0.8"
 OffsetArrays = "1"
-SIMDPirates = "0.8.16"
+SIMDPirates = "0.8.19"
 SLEEFPirates = "0.5.4"
 UnPack = "0,1"
-VectorizationBase = "0.12.24"
+VectorizationBase = "0.12.28"
 julia = "1.1"
 
 [extras]

src/condense_loopset.jl

@@ -293,7 +293,6 @@ end
 
 make_fast(q) = Expr(:macrocall, Symbol("@fastmath"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), q)
 make_crashy(q) = Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), q)
-make_fast_and_crashy(q) = q |> make_fast |> make_crashy
 
 function setup_call_inline(ls::LoopSet, inline::Int8 = zero(Int8), U::Int8 = zero(Int8), T::Int8 = zero(Int8))
     call = generate_call(ls, (inline,U,T))
@@ -335,5 +334,5 @@ function setup_call(ls::LoopSet, q = nothing, inline::Int8 = zero(Int8), check_e
     call = setup_call_inline(ls, inline, u₁, u₂)
     call = check_empty ? check_if_empty(ls, call) : call
     isnothing(q) && return Expr(:block, ls.prepreamble, call)
-    Expr(:block, ls.prepreamble, Expr(:if, check_args_call(ls), call, make_fast_and_crashy(q)))
+    Expr(:block, ls.prepreamble, Expr(:if, check_args_call(ls), call, make_crashy(make_fast(q))))
 end

src/lower_store.jl

@@ -57,7 +57,7 @@ function reduce_range!(q::Expr, toreduct::Symbol, instr::Instruction, Uh::Int, U
             push!(q.args, Expr(:(=), Symbol(toreduct, (u>>>1)), instrexpr))
         end
     else
-        for u ∈ Uh:Uh2-1
+        for u ∈ Uh:Uh2-2
             tru = Symbol(toreduct, u - Uh)
             instrexpr = callexpr(instr)
             push!(instrexpr.args, tru)

src/lowering.jl

@@ -652,11 +652,13 @@ function lsexpr(ls::LoopSet, q)
 end
 
 function calc_Ureduct(ls::LoopSet, us::UnrollSpecification)
-    @unpack u₁loopnum, u₁, u₂ = us
+    @unpack u₁loopnum, u₁, u₂, vectorizedloopnum = us
     if iszero(length(ls.outer_reductions))
         -1
     elseif u₂ == -1
-        min(u₁, 4)
+        loopisstatic = isstaticloop(getloop(ls, names(ls)[u₁loopnum]))
+        loopisstatic &= ((vectorizedloopnum != u₁loopnum) | (!iszero(ls.vector_width[])))
+        loopisstatic ? u₁ : min(u₁, 4)
     else
         8#u₂#u₁
     # elseif num_loops(ls) == u₁loopnum

test/dot.jl

@@ -210,6 +210,20 @@ using Test
         end
         s
     end
+    function dot33(a,b)
+        s = zero(eltype(a))
+        @avx for i ∈ 1:33
+            s += a[i] * b[i]
+        end
+        s
+    end
+    function dot17(a,b)
+        s = zero(eltype(a))
+        @avx for i ∈ 1:17
+            s += a[i] * b[i]
+        end
+        s
+    end
     # @macroexpand @_avx for i = 1:length(a_re) - 1
     #     c_re[i] = b_re[i] * a_re[i + 1] - b_im[i] * a_im[i + 1]
     #     c_im[i] = b_re[i] * a_im[i + 1] + b_im[i] * a_re[i + 1]
@@ -239,6 +253,9 @@ using Test
         @test myselfdot_avx(a) ≈ s
         @test myselfdotavx(a) ≈ s
 
+        @test dot17(a,b) ≈ @view(a[1:17])' * @view(b[1:17])
+        @test dot33(a,b) ≈ @view(a[1:33])' * @view(b[1:33])
+
         if T <: Union{Float32,Float64}
             πest = pi(a, b)
             @test πest == piavx(a, b)

test/miscellaneous.jl

@@ -1059,3 +1059,38 @@ end
 end
 
 
+function mul1!(y::Vector{T}, A::Matrix{UInt8}, x::Vector{T}) where T 
+    packedstride = size(A, 1)
+    m, n = size(A)
+    @avx for j ∈ eachindex(x)
+        for i ∈ eachindex(y)
+            k = 2 * ((i-1) & 3)
+            block = A[(j-1) * packedstride + ((i-1) >> 2) + 1]
+            Aij = (block >> k) & 3
+            y[i] += (((Aij >= 2) + (Aij >= 3))) * x[j]
+        end
+    end
+    y
+end
+function mul2!(y::Vector{T}, A::Matrix{UInt8}, x::Vector{T}) where T 
+    packedstride = size(A, 1)
+    m, n = size(A)
+    for j ∈ eachindex(x)
+        for i ∈ eachindex(y)
+            k = 2 * ((i-1) & 3)
+            block = A[(j-1) * packedstride + ((i-1) >> 2) + 1]
+            Aij = (block >> k) & 3
+            y[i] += (((Aij >= 2) + (Aij >= 3))) * x[j]
+        end
+    end
+    y
+end
+@testset "UInt8 mul" begin
+    for n in 1:200
+        v1 = rand(n); v3 =copy(v1);
+        v2 = rand(n);
+        A = rand(UInt8, (length(v1)>>2) + (length(v1)%4 != 0), length(v2))
+        @test mul1!(v1, A, v2) ≈ mul2!(v3, A, v2)
+    end
+end
+