Change handling of reductions to better accomdate their use within the loop body, and don't unroll inner reductions where the unrolled loop is the reduction loop if that reduction is already tiled.

Author: chriselrod

Manifest.toml

@@ -42,9 +42,9 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[SIMDPirates]]
 deps = ["VectorizationBase"]
-git-tree-sha1 = "1a9cbe1be1f5d43ac49eeb38ca64dd78e50a0cc6"
+git-tree-sha1 = "8f89aa38f5e4e89f2a474ffdc850fc21d6ab9ed4"
 uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a"
-version = "0.7.3"
+version = "0.7.4"
 
 [[SLEEFPirates]]
 deps = ["Libdl", "SIMDPirates", "VectorizationBase"]
@@ -69,6 +69,6 @@ version = "0.1.0"
 
 [[VectorizationBase]]
 deps = ["CpuId", "LinearAlgebra"]
-git-tree-sha1 = "83e32d4835fc4f4ecfd43eb59fa7fc00854b3d41"
+git-tree-sha1 = "2a83ab02d3fb6ad5de0c6d04103b0ca403d9a7d8"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.9.4"
+version = "0.9.5"

Project.toml

@@ -1,7 +1,7 @@
 name = "LoopVectorization"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 authors = ["Chris Elrod <[email protected]>"]
-version = "0.6.23"
+version = "0.6.24"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -13,10 +13,10 @@ VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 
 [compat]
 OffsetArrays = "1"
-SIMDPirates = "0.7.3"
+SIMDPirates = "0.7.4"
 SLEEFPirates = "0.4"
 UnPack = "0"
-VectorizationBase = "0.9.3"
+VectorizationBase = "0.9.5"
 julia = "1.1"
 
 [extras]

docs/src/assets/bench_AmulB_v1.svg

@@ -33,7 +33,7 @@ function setdiffv!(s4::AbstractVector{T}, s3::AbstractVector{T}, s1::AbstractVec
     end
 end
 function update_deps!(deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, parent::Operation)
-    mergesetv!(deps, loopdependencies(parent))#, reduceddependencies(parent))        
+    mergesetv!(deps, loopdependencies(parent))#, reduceddependencies(parent))
     if !(isload(parent) || isconstant(parent)) && !isreductcombineinstr(parent)
         mergesetv!(reduceddeps, reduceddependencies(parent))
     end
@@ -101,6 +101,18 @@ function isreductzero(op::Operation, ls::LoopSet, reduct_zero::Symbol)
     false
 end
 
+# function substitute_op_in_parents!(vparents::Vector{Operation}, replacer::Operation, replacee::Operation)
+#     for i ∈ eachindex(vparents)
+#         opp = vparents[i]
+#         if opp === replacee
+#             vparents[i] = replacer
+#         else
+#             substitute_op_in_parents!(parents(opp), replacer, replacee)
+#         end
+#     end
+# end
+
+
 function add_reduction_update_parent!(
     vparents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet,
     parent::Operation, instr::Symbol, directdependency::Bool, elementbytes::Int
@@ -109,11 +121,14 @@ function add_reduction_update_parent!(
     isouterreduction = parent.instruction === LOOPCONSTANT
     Instr = instruction(ls, instr)
     instrclass = reduction_instruction_class(Instr) # key allows for faster lookups
+    reduct_zero = reduction_zero(instrclass)
     # if parent is not an outer reduction...
-    if !isouterreduction
-        # and parent is not a reduction_zero
-        reduct_zero = reduction_zero(instrclass)
+    # if !isouterreduction && !isreductzero(parent, ls, reduct_zero)
+    add_reduct_instruct = !isouterreduction && !isconstant(parent)
+    if add_reduct_instruct
+        # We add 
         reductcombine = reduction_scalar_combine(instrclass)
+        # reductcombine = :identity
         reductsym = gensym(:reduction)
         reductinit = add_constant!(ls, gensym(:reductzero), loopdependencies(parent), reductsym, elementbytes, :numericconstant)
         if reduct_zero === :zero
@@ -128,13 +143,13 @@ function add_reduction_update_parent!(
             end
             pushpreamble!(ls, op, name, reductinit)
         end
-        if isreductzero(parent, ls, reduct_zero)
-            reductcombine = reduction_combine_to(instrclass)
-        end
+        # if 
+            # reductcombine = reduction_combine_to(instrclass)
+        # end
     else
         reductinit = parent
         reductsym = var
-        reductcombine = Symbol("")
+        reductcombine = :identity#Symbol("")
     end
     combineddeps = copy(deps); mergesetv!(combineddeps, reduceddeps)
     # directdependency && pushparent!(vparents, deps, reduceddeps, reductinit)#parent) # deps and reduced deps will not be disjoint
@@ -145,16 +160,20 @@ function add_reduction_update_parent!(
         else
             push!(vparents, reductinit)
         end
+    # elseif !isouterreduction
+        # substitute_op_in_parents!(vparents, reductinit, parent)
     end
     update_reduction_status!(vparents, reduceddeps, name(reductinit))
     # this is the op added by add_compute
     op = Operation(length(operations(ls)), reductsym, elementbytes, instr, compute, deps, reduceddeps, vparents)
     parent.instruction === LOOPCONSTANT && push!(ls.outer_reductions, identifier(op))
     opout = pushop!(ls, op, var) # note this overwrites the entry in the operations dict, but not the vector
     # isouterreduction || iszero(length(reduceddeps)) && return opout
+    # return opout
     isouterreduction && return opout
     # create child op, which is the reduction combination
-    childrdeps = Symbol[]; childparents = Operation[ op, parent ]
+    childrdeps = Symbol[]; childparents = Operation[ op ]#, parent ]
+    add_reduct_instruct && push!(childparents, parent)
     childdeps = loopdependencies(reductinit)
     setdiffv!(childrdeps, loopdependencies(op), childdeps)
     child = Operation(

docs/src/assets/bench_AmulBt_v1.svg

@@ -198,9 +198,9 @@ const COST = Dict{Instruction,InstructionCost}(
     Instruction(:adjoint) => InstructionCost(0,0.0,0.0,0),
     Instruction(:transpose) => InstructionCost(0,0.0,0.0,0),
     Instruction(:prefetch) => InstructionCost(0,0.0,0.0,0),
+    Instruction(:prefetch0) => InstructionCost(0,0.0,0.0,0),
     Instruction(:prefetch1) => InstructionCost(0,0.0,0.0,0),
-    Instruction(:prefetch2) => InstructionCost(0,0.0,0.0,0),
-    Instruction(:prefetch3) => InstructionCost(0,0.0,0.0,0)
+    Instruction(:prefetch2) => InstructionCost(0,0.0,0.0,0)
 )
 @inline prefetch0(x, i) = SIMDPirates.prefetch(gep(stridedpointer(x), (extract_data(i) - 1,)), Val{3}(), Val{0}())
 @inline prefetch0(x, i, j) = SIMDPirates.prefetch(gep(stridedpointer(x), (extract_data(i) - 1, extract_data(j) - 1)), Val{3}(), Val{0}())
@@ -274,7 +274,7 @@ reduction_instruction_class(instr::Symbol) = get(REDUCTION_CLASS, instr, NaN)
 reduction_instruction_class(instr::Instruction) = reduction_instruction_class(instr.instr)
 function reduction_to_single_vector(x::Float64)
     # x == 1.0 ? :evadd : x == 2.0 ? :evmul : x == 3.0 ? :vor : x == 4.0 ? :vand : x == 5.0 ? :max : x == 6.0 ? :min : throw("Reduction not found.")
-    x == 1.0 ? :evadd : x == 2.0 ? :evmul : x == 5.0 ? :max : x == 6.0 ? :min : throw("Reduction not found.")
+    x == ADDITIVE_IN_REDUCTIONS ? :evadd : x == MULTIPLICATIVE_IN_REDUCTIONS ? :evmul : x == MAX ? :max : x == MIN ? :min : throw("Reduction not found.")
 end
 reduction_to_single_vector(x) = reduction_to_single_vector(reduction_instruction_class(x))
 # function reduction_to_scalar(x::Float64)
@@ -284,17 +284,17 @@ reduction_to_single_vector(x) = reduction_to_single_vector(reduction_instruction
 # reduction_to_scalar(x) = reduction_to_scalar(reduction_instruction_class(x))
 function reduction_scalar_combine(x::Float64)
     # x == 1.0 ? :reduced_add : x == 2.0 ? :reduced_prod : x == 3.0 ? :reduced_any : x == 4.0 ? :reduced_all : x == 5.0 ? :reduced_max : x == 6.0 ? :reduced_min : throw("Reduction not found.")
-    x == 1.0 ? :reduced_add : x == 2.0 ? :reduced_prod : x == 5.0 ? :reduced_max : x == 6.0 ? :reduced_min : throw("Reduction not found.")
+    x == ADDITIVE_IN_REDUCTIONS ? :reduced_add : x == MULTIPLICATIVE_IN_REDUCTIONS ? :reduced_prod : x == MAX ? :reduced_max : x == MIN ? :reduced_min : throw("Reduction not found.")
 end
 reduction_scalar_combine(x) = reduction_scalar_combine(reduction_instruction_class(x))
 function reduction_combine_to(x::Float64)
     # x == 1.0 ? :reduce_to_add : x == 2.0 ? :reduce_to_prod : x == 3.0 ? :reduce_to_any : x == 4.0 ? :reduce_to_all : x == 5.0 ? :reduce_to_max : x == 6.0 ? :reduce_to_min : throw("Reduction not found.")
-    x == 1.0 ? :reduce_to_add : x == 2.0 ? :reduce_to_prod : x == 5.0 ? :reduce_to_max : x == 6.0 ? :reduce_to_min : throw("Reduction not found.")
+    x == ADDITIVE_IN_REDUCTIONS ? :reduce_to_add : x == MULTIPLICATIVE_IN_REDUCTIONS ? :reduce_to_prod : x == MAX ? :reduce_to_max : x == MIN ? :reduce_to_min : throw("Reduction not found.")
 end
 reduction_combine_to(x) = reduction_combine_to(reduction_instruction_class(x))
 function reduction_zero(x::Float64)
     # x == 1.0 ? :zero : x == 2.0 ? :one : x == 3.0 ? :false : x == 4.0 ? :true : x == 5.0 ? :typemin : x == 6.0 ? :typemax : throw("Reduction not found.")
-    x == 1.0 ? :zero : x == 2.0 ? :one : x == 5.0 ? :typemin : x == 6.0 ? :typemax : throw("Reduction not found.")
+    x == ADDITIVE_IN_REDUCTIONS ? :zero : x == MULTIPLICATIVE_IN_REDUCTIONS ? :one : x == MAX ? :typemin : x == MIN ? :typemax : throw("Reduction not found.")
 end
 reduction_zero(x) = reduction_zero(reduction_instruction_class(x))
 

docs/src/assets/bench_Amulvb_v1.svg

@@ -50,18 +50,18 @@ function register_pressure(op::Operation)
     end
 end
 function cost(ls::LoopSet, op::Operation, vectorized::Symbol, Wshift::Int, size_T::Int = op.elementbytes)
-    isconstant(op) && return 0.0, 0, Int(length(loopdependencies(op)) > 0)
-    isloopvalue(op) && return 0.0, 0, 0
+    isconstant(op) && return 0.0, 0, Float64(length(loopdependencies(op)) > 0)
+    isloopvalue(op) && return 0.0, 0, 0.0
     # Wshift == dependson(op, vectorized) ? Wshift : 0
     # c = first(cost(instruction(op), Wshift, size_T))::Int
     instr = Instruction(:LoopVectorization, instruction(op).instr)
     # instr = instruction(op)
     if length(parents(op)) == 1
         if instr == Instruction(:-) || instr === Instruction(:vsub) || instr == Instruction(:+) || instr == Instruction(:vadd)
-            return 0.0, 0, 0
+            return 0.0, 0, 0.0
         end
     elseif iscompute(op) && all(opp -> (isloopvalue(opp) | isconstant(opp)), parents(op))
-        return 0.0, 0, 0
+        return 0.0, 0, 0.0
     end
     opisvectorized = dependson(op, vectorized)
     srt, sl, srp = opisvectorized ? vector_cost(instr, Wshift, size_T) : scalar_cost(instr)
@@ -80,7 +80,7 @@ function cost(ls::LoopSet, op::Operation, vectorized::Symbol, Wshift::Int, size_
             sl *= 3
         end
     end
-    srt, sl, srp
+    srt, sl, Float64(srp)
 end
 
     # Base._return_type()
@@ -336,7 +336,7 @@ function maybedemotesize(U::Int, N::Int)
     um1rep > urep ? U : Um1
 end
 function maybedemotesize(T::Int, N::Int, U::Int, Uloop::Loop, maxTbase::Int)
-    T > 1 || return T
+    T > 1 || return 1
     T == N && return T
     T = maybedemotesize(T, N)
     if !(isstaticloop(Uloop) && length(Uloop) == U)
@@ -451,11 +451,11 @@ function convolution_cost_factor(ls::LoopSet, op::Operation, u1::Symbol, u2::Sym
             if isstaticloop(loop) && length(loop) ≤ 4
                 itersym = loop.itersymbol
                 if itersym !== u1 && itersym !== u2
-                    return (0.25, 1.0)
+                    return (0.25, VectorizationBase.REGISTER_COUNT == 32 ? 2.0 : 1.0)
                 end
             end
         end
-        (0.25, 0.5)
+        (0.25, VectorizationBase.REGISTER_COUNT == 32 ? 1.2 : 1.0)
     else
         (1.0, 1.0)
     end
@@ -529,7 +529,7 @@ function evaluate_cost_tile(
             factor1, factor2 = convolution_cost_factor(ls, op, unrolled, tiled, vectorized)
             rt *= factor1; rp *= factor2;
         end
-        # @show op rt, lat, rp
+        # @show isunrolled, istiled, op rt, lat, rp
         rp = opisininnerloop ? rp : 0 # we only care about register pressure within the inner most loop
         rt *= iters[id]
         if isunrolled && istiled # no cost decrease; cost must be repeated
@@ -546,6 +546,7 @@ function evaluate_cost_tile(
             reg_pressure[4] += rp
         end
     end
+    # @show reg_pressure
     costpenalty = (sum(reg_pressure) > VectorizationBase.REGISTER_COUNT) ? 2 : 1
     # @show order, vectorized cost_vec reg_pressure
     # @show solve_tilesize(ls, unrolled, tiled, cost_vec, reg_pressure)
@@ -665,12 +666,15 @@ function choose_order(ls::LoopSet)
     end
 end
 
-function register_pressure(ls::LoopSet)
-    order, unroll, vec, U, T = choose_order(ls)
+function register_pressure(ls::LoopSet, U, T)
     if T == -1
         sum(register_pressure, operations(ls))
     else
-        rp = @view ls.reg_pressure[:,1]
-        tU * tT * rp[1] + tU * rp[2] + rp[3] + rp[4]
+        rp = @view ls.reg_pres[:,1]
+        U * T * rp[1] + U * rp[2] + rp[3] + rp[4]
     end
 end
+function register_pressure(ls::LoopSet)
+    order, unroll, tile, vec, U, T = choose_order(ls)
+    register_pressure(ls, U, T)
+end

docs/src/assets/bench_AplusAt_v1.svg

@@ -1,41 +1,4 @@
 
-# """
-# ShortVector{T} simply wraps a Vector{T}, but uses a different hash function that is faster for short vectors to support using it as the keys of a Dict.
-# This hash function scales O(N) with length of the vectors, so it is slow for long vectors.
-# """
-# struct ShortVector{T} <: DenseVector{T}
-#     data::Vector{T}
-# end
-# Base.@propagate_inbounds Base.getindex(x::ShortVector, I...) = x.data[I...]
-# Base.@propagate_inbounds Base.setindex!(x::ShortVector, v, I...) = x.data[I...] = v
-# ShortVector{T}(::UndefInitializer, N::Integer) where {T} = ShortVector{T}(Vector{T}(undef, N))
-# @inbounds Base.length(x::ShortVector) = length(x.data)
-# @inbounds Base.size(x::ShortVector) = size(x.data)
-# @inbounds Base.strides(x::ShortVector) = strides(x.data)
-# @inbounds Base.push!(x::ShortVector, v) = push!(x.data, v)
-# @inbounds Base.append!(x::ShortVector, v) = append!(x.data, v)
-# function Base.hash(x::ShortVector, h::UInt)
-#     @inbounds for n ∈ eachindex(x)
-#         h = hash(x[n], h)
-#     end
-#     h
-# end
-# function Base.isequal(a::ShortVector{T}, b::ShortVector{T}) where {T}
-#     length(a) == length(b) || return false
-#     @inbounds for i ∈ 1:length(a)
-#         a[i] === b[i] || return false
-#     end
-#     true
-# end
-# Base.convert(::Type{Vector}, sv::ShortVector) = sv.data
-# Base.convert(::Type{Vector{T}}, sv::ShortVector{T}) where {T} = sv.data
-
-
-# For passing options like array types and mask
-# struct LoopSetOptions
-
-# end
-
 struct UnrollSpecification
     unrolledloopnum::Int
     tiledloopnum::Int

docs/src/assets/bench_AtmulB_v1.svg

@@ -8,9 +8,10 @@ function lower_compute!(
     opunrolled = unrolled ∈ loopdependencies(op)
 )
 
-    var = op.variable
-    mvar = mangledvar(op)
+    var = name(op)
+    instr = instruction(op)
     parents_op = parents(op)
+    mvar = mangledvar(op)
     nparents = length(parents_op)
     parentstiled = if suffix === nothing
         optiled = false
@@ -25,13 +26,20 @@ function lower_compute!(
         optiled = true
         [tiled ∈ loopdependencies(opp) for opp ∈ parents_op]
     end
-    parentsunrolled = [unrolled ∈ loopdependencies(opp) || unrolled ∈ reducedchildren(opp) for opp ∈ parents_op]
+    parentsunrolled = isunrolled_sym.(parents_op, unrolled, tiled)
+    if instr.instr === :identity && name(first(parents_op)) === var && isone(length(parents_op))
+        if (opunrolled == first(parentsunrolled)) && ((!isnothing(suffix)) == first(parentstiled))
+            return
+        end
+    end
+    unrollsym = isunrolled_sym(op, unrolled, suffix)
     if !opunrolled && any(parentsunrolled)
         parents_op = copy(parents_op)
         for i ∈ eachindex(parentsunrolled)
             parentsunrolled[i] || continue
             parentsunrolled[i] = false
             parentop = parents_op[i]
+            # @show op, parentop
             newparentop = Operation(
                 parentop.identifier, gensym(parentop.variable), parentop.elementbytes, parentop.instruction, parentop.node_type,
                 parentop.dependencies, parentop.reduced_deps, parentop.parents, parentop.ref, parentop.reduced_children
@@ -43,14 +51,18 @@ function lower_compute!(
                 parentname = Symbol(parentname, suffix_)
                 newparentname = Symbol(newparentname, suffix_)
             end
-            for u ∈ 0:U-1
-                push!(q.args, Expr(:(=), Symbol(newparentname, u), Symbol(parentname, u)))
+            if isconstant(newparentop)
+                push!(q.args, Expr(:(=), newparentname, Symbol(parentname, 0)))
+                continue
+            else
+                for u ∈ 0:U-1
+                    push!(q.args, Expr(:(=), Symbol(newparentname, u), Symbol(parentname, u)))
+                end
+                reduce_expr!(q, newparentname, Instruction(reduction_to_single_vector(instruction(newparentop))), U)
+                push!(q.args, Expr(:(=), newparentname, Symbol(newparentname, 0)))
             end
-            reduce_expr!(q, newparentname, Instruction(reduction_to_single_vector(instruction(newparentop))), U)
-            push!(q.args, Expr(:(=), newparentname, Symbol(newparentname, 0)))
         end
     end
-    instr = op.instruction
     # cache unroll and tiling check of parents
     # not broadcasted, because we use frequent checks of individual bools
     # making BitArrays inefficient.
@@ -84,7 +96,7 @@ function lower_compute!(
             # modsuffix = ((u + suffix*U) & 3)
             modsuffix = (suffix & 3)
             Symbol(mvar, modsuffix)
-        elseif opunrolled
+        elseif unrollsym
             Symbol(mvar, u)
         else
             mvar
@@ -125,7 +137,11 @@ function lower_compute!(
                 continue
             end
         end
-        push!(q.args, Expr(:(=), varsym, instrcall))
+        if instr.instr === :identity && isone(length(parents_op))
+            push!(q.args, Expr(:(=), varsym, instrcall.args[2]))
+        else
+            push!(q.args, Expr(:(=), varsym, instrcall))
+        end
     end
 end