Improve prefetching behavior slightly, update README to note loops will generally be faster than broadcasts.

Author: chriselrod

README.md

@@ -213,17 +213,56 @@ julia> buf1 = Matrix{Float64}(undef, size(C,1), size(C,2));
 julia> buf2 = similar(X1);
 
 julia> @btime $X1 .= view($A,1,:) .+ mul!($buf2, $B, ($buf1 .= $C .+ $D'));
-  7.896 μs (0 allocations: 0 bytes)
+  9.188 μs (0 allocations: 0 bytes)
 
 julia> @btime @avx $X2 .= view($A,1,:) .+ $B .*ˡ ($C .+ $D');
-  7.647 μs (0 allocations: 0 bytes)
+  6.751 μs (0 allocations: 0 bytes)
+
+julia> @test X1 ≈ X2
+Test Passed
+
+julia> AmulBtest!(X1, B, C, D, view(A,1,:))
+
+julia> AmulBtest2!(X2, B, C, D, view(A,1,:))
 
 julia> @test X1 ≈ X2
 Test Passed
 ```
 The lazy matrix multiplication operator `*ˡ` escapes broadcasts and fuses, making it easy to write code that avoids intermediates. However, I would recomend always checking if splitting the operation into pieces, or at least isolating the matrix multiplication, increases performance. That will often be the case, especially if the matrices are large, where a separate multiplication can leverage BLAS (and perhaps take advantage of threads).
 This may improve as the optimizations within LoopVectorization improve.
 
+Note that loops will be faster than broadcasting in general. This is because the behavior of broadcasts is determined by runtime information (i.e., dimensions other than the leading dimension of size `1` will be broadcasted; it is not known which these will be at compile time).
+```julia
+julia> function AmulBtest!(C,A,Bk,Bn,d)
+          @avx for m ∈ axes(A,1), n ∈ axes(Bk,2)
+             ΔCₘₙ = zero(eltype(C))
+             for k ∈ axes(A,2)
+                ΔCₘₙ += A[m,k] * (Bk[k,n] + Bn[n,k])
+             end
+             C[m,n] = ΔCₘₙ + d[m]
+          end
+       end
+AmulBtest! (generic function with 1 method)
+
+julia> AmulBtest!(X2, B, C, D, view(A,1,:))
+
+julia> @test X1 ≈ X2
+Test Passed
+
+julia> @benchmark AmulBtest!($X2, $B, $C, $D, view($A,1,:))
+BenchmarkTools.Trial:
+  memory estimate:  0 bytes
+  allocs estimate:  0
+  --------------
+  minimum time:     5.793 μs (0.00% GC)
+  median time:      5.816 μs (0.00% GC)
+  mean time:        5.824 μs (0.00% GC)
+  maximum time:     14.234 μs (0.00% GC)
+  --------------
+  samples:          10000
+  evals/sample:     6
+```
+
 </p>
 </details>
 

src/determinestrategy.jl

@@ -315,6 +315,7 @@ function solve_unroll(X, R, u₁L, u₂L, u₁step, u₂step)
     if !(isfinite(u₂float) & isfinite(u₁float)) # brute force
         u₁low = u₂low = 1
         u₁high = u₂high = REGISTER_COUNT == 32 ? 10 : 6#8
+        println("Fail")
         return solve_unroll_iter(X, R, u₁L, u₂L, u₁low:u₁step:u₁high, u₂low:u₂step:u₂high)
     end
     u₁low = floor(Int, u₁float)

src/lower_compute.jl

@@ -1,9 +1,22 @@
 
-function load_constrained(op, u₁loop, u₂loop)
-    unrolleddeps = Symbol[]
+function load_constrained(op, u₁loop, u₂loop, forprefetch = false)
     loopdeps = loopdependencies(op)
-    u₁loop ∈ loopdeps && push!(unrolleddeps, u₁loop)
-    u₂loop ∈ loopdeps && push!(unrolleddeps, u₂loop)
+    dependsonu₁ = u₁loop ∈ loopdeps
+    if u₂loop === Symbol("##undefined##")
+        if forprefetch
+            dependsonu₁ || return false
+        end
+        # unrolleddeps = [ u₁loop ]
+    else        
+        dependsonu₂ = u₂loop ∈ loopdeps
+        if forprefetch
+            (dependsonu₁ & dependsonu₂) || return false
+        end
+        # unrolleddeps = [ u₁loop, u₂loop ]
+    end
+    unrolleddeps = Symbol[]
+    dependsonu₁ && push!(unrolleddeps, u₁loop)
+    dependsonu₂ && push!(unrolleddeps, u₂loop)
     any(opp -> isload(opp) && all(in(loopdependencies(opp)), unrolleddeps), parents(op))
 end
 

src/lower_load.jl

@@ -89,7 +89,7 @@ function prefetchisagoodidea(ls::LoopSet, op::Operation, td::UnrollArgs)
             if prod(s -> length(getloop(ls, s)), @view(indices[1:innermostloopind-1])) ≥ 120 && length(getloop(ls, innermostloopsym)) ≥ 120
                 if op.ref.ref.offsets[innermostloopind] < 120
                     for opp ∈ operations(ls)
-                        iscompute(opp) && (innermostloopsym ∈ loopdependencies(opp)) && load_constrained(opp, u₁loopsym, u₂loopsym) && return 0
+                        iscompute(opp) && (innermostloopsym ∈ loopdependencies(opp)) && load_constrained(opp, u₁loopsym, u₂loopsym, true) && return 0
                     end
                     return innermostloopind
                 end