Update for latest VectorizationBase and SIMDPirates using 1-based indexing, so that LoopVectorization now uses the same indexing as the original arrays, hopefully simplifying code generation.

Author: chriselrod

Manifest.toml

@@ -32,9 +32,9 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[OffsetArrays]]
-git-tree-sha1 = "707e34562700b81e8aa13548eb6b23b18112e49b"
+git-tree-sha1 = "6a35d9446b40ae5004cd7bd0f1ae3505528c7fd6"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.0.2"
+version = "1.0.3"
 
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
@@ -54,15 +54,15 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[SIMDPirates]]
 deps = ["VectorizationBase"]
-git-tree-sha1 = "f91198b7ef74b04028f98e0eed7c556b93538a2e"
+git-tree-sha1 = "8705f93f66f789605baa670e8df0244ad4c80d40"
 uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a"
-version = "0.6.6"
+version = "0.7.0"
 
 [[SLEEFPirates]]
 deps = ["Libdl", "SIMDPirates", "VectorizationBase"]
-git-tree-sha1 = "66014d5d3f77753a70b706aeaede1ad70a942915"
+git-tree-sha1 = "d53dea0b025c01e5a40c827236505507d4572a05"
 uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
-version = "0.4.1"
+version = "0.4.2"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -76,6 +76,6 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[VectorizationBase]]
 deps = ["CpuId", "LinearAlgebra"]
-git-tree-sha1 = "8abb5697fb64cadccd1bba444c955942d3181e5c"
+git-tree-sha1 = "34004c00c561c6d4de6d30970bbd9c0a5c2be0a3"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.7.1"
+version = "0.8.1"

Project.toml

@@ -14,9 +14,9 @@ VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 [compat]
 OffsetArrays = "1"
 Parameters = "0"
-SIMDPirates = "~0.6.6"
-SLEEFPirates = "~0.4"
-VectorizationBase = "~0.7.1"
+SIMDPirates = "0.7"
+SLEEFPirates = "0.4"
+VectorizationBase = "0.8"
 julia = "1.1"
 
 [extras]

src/LoopVectorization.jl

@@ -5,7 +5,8 @@ using VectorizationBase: REGISTER_SIZE, REGISTER_COUNT, extract_data, num_vector
     mask, masktable, pick_vector_width_val, valmul, valrem, valmuladd, valadd, valsub, _MM,
     maybestaticlength, maybestaticsize, staticm1, subsetview, vzero, stridedpointer_for_broadcast,
     Static, StaticUnitRange, StaticLowerUnitRange, StaticUpperUnitRange, unwrap, maybestaticrange,
-    PackedStridedPointer, SparseStridedPointer, RowMajorStridedPointer, StaticStridedPointer, StaticStridedStruct
+    PackedStridedPointer, SparseStridedPointer, RowMajorStridedPointer, StaticStridedPointer, StaticStridedStruct,
+    maybestaticfirst, maybestaticlast
 using SIMDPirates: VECTOR_SYMBOLS, evadd, evmul, vrange, reduced_add, reduced_prod, reduce_to_add, reduce_to_prod,
     sizeequivalentfloat, sizeequivalentint, vadd!, vsub!, vmul!, vfdiv!, vfmadd!, vfnmadd!, vfmsub!, vfnmsub!,
     vfmadd231, vfmsub231, vfnmadd231, vfnmsub231, #prefetch,

src/add_loads.jl

@@ -56,17 +56,13 @@ end
 
 function add_loopvalue!(ls::LoopSet, arg::Symbol, elementbytes::Int)
     # check for CSE opportunity
-    loopsym = Symbol("##LOOPSYMBOL##", arg)
-    ar = ArrayReference(loopsym, [arg])
-    for op ∈ operations(ls)
-        if isload(op) && op.ref.ref == ar
-            return op
-        end
+    instr = Instruction(arg, arg)
+    for op ∈ operations(ls)#check to CSE
+        (op.variable === arg && instr == instruction(op)) && return op
     end
-    pushpreamble!(ls, Expr(:(=), loopsym, LoopValue()))
-    loopsymop = add_simple_load!(ls, gensym(loopsym), ar, elementbytes, false)
-    push!(ls.syms_aliasing_refs, name(loopsymop))
-    push!(ls.refs_aliasing_syms, loopsymop.ref)
-    loopsymop
+    op = Operation(
+        length(operations(ls)), arg, elementbytes, instr, loopvalue, [arg]
+    )
+    pushop!(ls, op, arg)
 end
 

src/broadcast.jl

@@ -69,7 +69,7 @@ function add_broadcast!(
     pushpreamble!(ls, Expr(:(=), mB, Expr(:(.), bcname, QuoteNode(:b))))
     pushpreamble!(ls, Expr(:(=), K, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), Expr(:ref, Expr(:call, :size, mB), 1))))
     k = gensym(:k)
-    add_loop!(ls, Loop(k, 0, K), k)
+    add_loop!(ls, Loop(k, 1, K), k)
     m = loopsyms[1];
     if numdims(B) == 1
         bloopsyms = Symbol[k]
@@ -226,7 +226,7 @@ end
     sizes = Expr(:tuple)
     for (n,itersym) ∈ enumerate(loopsyms)
         Nsym = gensym(:N)
-        add_loop!(ls, Loop(itersym, 0, Nsym), itersym)
+        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
         push!(sizes.args, Nsym)
     end
     pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest)))
@@ -251,7 +251,7 @@ end
     sizes = Expr(:tuple)
     for (n,itersym) ∈ enumerate(loopsyms)
         Nsym = gensym(:N)
-        add_loop!(ls, Loop(itersym, 0, Nsym), itersym)
+        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
         push!(sizes.args, Nsym)
     end
     pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest′)))

src/condense_loopset.jl

@@ -46,7 +46,7 @@ function ArrayRefStruct(ls::LoopSet, mref::ArrayReferenceMeta, arraysymbolinds::
     ArrayRefStruct( index_types, indices, offsets )
 end
 
-struct OperationStruct
+struct OperationStruct <: AbstractLoopOperation
     # instruction::Instruction
     loopdeps::UInt64
     reduceddeps::UInt64
@@ -56,10 +56,8 @@ struct OperationStruct
     array::UInt8
     symid::UInt8
 end
-isload(os::OperationStruct) = os.node_type == memload
-isstore(os::OperationStruct) = os.node_type == memstore
-iscompute(os::OperationStruct) = os.node_type == compute
-isconstant(os::OperationStruct) = os.node_type == constant
+optype(os) = os.node_type
+
 function findmatchingarray(ls::LoopSet, mref::ArrayReferenceMeta)
     id = 0x01
     for r ∈ ls.refs_aliasing_syms
@@ -213,12 +211,8 @@ end
         assigned_names[i] = LHS
         d = (D[i])::Union{Nothing,Int}
         if d === nothing # stridedpointer
-            if ari == -1
-                RHS = Expr(:call, :LoopValue)
-            else
-                num_arrays += 1
-                RHS = Expr(:call, lv(:stridedpointer), Expr(:ref, :vargs, ari), Expr(:ref, :arraydescript, ari))
-            end
+            num_arrays += 1
+            RHS = Expr(:call, lv(:stridedpointer), Expr(:ref, :vargs, ari), Expr(:ref, :arraydescript, ari))
         else #subsetview
             j += 1
             RHS = Expr(:call, :subsetview, assigned_names[ari], Expr(:call, Expr(:curly, :Val, d)), Expr(:ref, :subsetvals, j))
@@ -317,8 +311,6 @@ function setup_call_noinline(ls::LoopSet, U = zero(Int8), T = zero(Int8))
                     push!(stridedpointerLHS, vp)
                     push!(vptrindices.args, findfirst(a -> vptr(a) == vp, ls.refs_aliasing_syms))
                 end
-            elseif ex.args[2] == LoopValue()
-                push!(loopvalueLHS, first(ex.args))
             end
         end
         push!(q.args, ex)

src/determinestrategy.jl

@@ -1,7 +1,46 @@
 
-# TODO: FIXME for general case
-# wrong for transposed matrices, and certain views/SubArrays.
-unitstride(op::Operation, s) = first(getindices(op)) === s
+function indexappearences(op::Operation, s::Symbol)
+    s ∉ loopdependencies(op) && return 0
+    appearences = 0
+    if isloopvalue(op)
+        return s === first(loopdependencies(op)) ? 1 : 0
+    elseif isload(op)
+        return 100
+    end
+    newapp = 0
+    for opp ∈ parents(op)
+        newapp += indexappearences(opp, s)
+    end
+    factor = instruction(op).instr ∈ (:+, :vadd, :add_fast, :evadd) ? 1 : 10
+    newapp * factor
+end
+function findparent(ls::LoopSet, s::Symbol)#opdict isn't filled when reconstructing
+    id = findfirst(op -> name(op) === s, operations(ls))
+    id === nothing && throw("$s not found")
+    operations(ls)[id]
+end
+function unitstride(ls::LoopSet, op::Operation, s::Symbol)
+    inds = getindices(op)
+    li = op.ref.loopedindex
+    # The first index is allowed to be indexed by `s`
+    fi = first(inds)
+    if fi === Symbol("##DISCONTIGUOUSSUBARRAY##")
+        return false
+    elseif !first(li)
+        # We must check if this
+        parent = findparent(ls, fi)
+        indexappearences(parent, s) > 1 && return false
+    end
+    for i ∈ 2:length(inds)
+        if li[i]
+            s === inds[i] && return false
+        else
+            parent = findparent(ls, inds[i])
+            s ∈ loopdependencies(parent) && return false
+        end
+    end
+    true
+end
 
 function register_pressure(op::Operation)
     if isconstant(op)
@@ -10,7 +49,7 @@ function register_pressure(op::Operation)
         instruction_cost(instruction(op)).register_pressure
     end
 end
-function cost(op::Operation, unrolled::Symbol, Wshift::Int, size_T::Int = op.elementbytes)
+function cost(ls::LoopSet, op::Operation, unrolled::Symbol, Wshift::Int, size_T::Int = op.elementbytes)
     isconstant(op) && return 0.0, 0, 1
     # Wshift == dependson(op, unrolled) ? Wshift : 0
     # c = first(cost(instruction(op), Wshift, size_T))::Int
@@ -27,7 +66,7 @@ function cost(op::Operation, unrolled::Symbol, Wshift::Int, size_T::Int = op.ele
         # either vbroadcast/reductionstore, vmov(a/u)pd, or gather/scatter
         # @show instr, unrolled, loopdependencies(op), unitstride(op, unrolled)
         if opisunrolled
-            if !unitstride(op, unrolled)# || !isdense(op) # need gather/scatter
+            if !unitstride(ls, op, unrolled)# || !isdense(op) # need gather/scatter
                 r = (1 << Wshift)
                 srt *= r
                 sl *= r
@@ -93,7 +132,7 @@ function evaluate_cost_unroll(
             hasintersection(rd, nested_loop_syms[1:end-length(rd)]) && return Inf
             included_vars[id] = true
             # @show op first(cost(op, vectorized, Wshift, size_T)), iter
-            total_cost += iter * first(cost(op, vectorized, Wshift, size_T))
+            total_cost += iter * first(cost(ls, op, vectorized, Wshift, size_T))
             total_cost > max_cost && return total_cost # abort if more expensive; we only want to know the cheapest
         end
     end
@@ -102,18 +141,18 @@ end
 
 # only covers vectorized ops; everything else considered lifted?
 function depchain_cost!(
-    skip::Vector{Bool}, op::Operation, vectorized::Symbol, Wshift::Int, size_T::Int, rt::Float64 = 0.0, sl::Int = 0
+    ls::LoopSet, skip::Vector{Bool}, op::Operation, vectorized::Symbol, Wshift::Int, size_T::Int, rt::Float64 = 0.0, sl::Int = 0
 )
     skip[identifier(op)] = true
     # depth first search
     for opp ∈ parents(op)
         skip[identifier(opp)] && continue
-        rt, sl = depchain_cost!(skip, opp, vectorized, Wshift, size_T, rt, sl)
+        rt, sl = depchain_cost!(ls, skip, opp, vectorized, Wshift, size_T, rt, sl)
     end
     # Basically assuming memory and compute don't conflict, but everything else does
     # Ie, ignoring the fact that integer and floating point operations likely don't either
     if iscompute(op)
-        rtᵢ, slᵢ = cost(op, vectorized, Wshift, size_T)
+        rtᵢ, slᵢ = cost(ls, op, vectorized, Wshift, size_T)
         rt += rtᵢ; sl += slᵢ
     end
     rt, sl
@@ -139,9 +178,9 @@ function unroll_no_reductions(ls, order, vectorized, Wshift, size_T)
     for op ∈ operations(ls)
         dependson(op, innermost) || continue
         if iscompute(op)
-            compute_rt += first(cost(op, vectorized, Wshift, size_T))
+            compute_rt += first(cost(ls, op, vectorized, Wshift, size_T))
         elseif isload(op)
-            load_rt += first(cost(op, vectorized, Wshift, size_T))
+            load_rt += first(cost(ls, op, vectorized, Wshift, size_T))
         end
     end
     # heuristic guess
@@ -181,13 +220,13 @@ function determine_unroll_factor(
     for op ∈ operations(ls)
         dependson(op, unrolled) || continue
         if isreduction(op)
-            rt, sl = depchain_cost!(visited_nodes, op, vectorized, Wshift, size_T)
+            rt, sl = depchain_cost!(ls, visited_nodes, op, vectorized, Wshift, size_T)
             latency = max(sl, latency)
             compute_recip_throughput += rt
         elseif isload(op)
-            load_recip_throughput += first(cost(op, vectorized, Wshift, size_T))
+            load_recip_throughput += first(cost(ls, op, vectorized, Wshift, size_T))
         elseif isstore(op)
-            store_recip_throughput += first(cost(op, vectorized, Wshift, size_T))
+            store_recip_throughput += first(cost(ls, op, vectorized, Wshift, size_T))
         end
     end
     recip_throughput = max(
@@ -424,7 +463,7 @@ function evaluate_cost_tile(
         opisininnerloop = descendentsininnerloop[id]
         isunrolled = unrolledtiled[1,id]
         istiled = unrolledtiled[2,id]
-        rt, lat, rp = cost(op, vectorized, Wshift, size_T)
+        rt, lat, rp = cost(ls, op, vectorized, Wshift, size_T)
         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