Broadcasting is much slower than a for loop

[YingboMa]

Here is a minimal working example.

julia> using BenchmarkTools

julia> function foo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}, e::Vector{T}) where T
           @. a = b + 0.1 * (0.2c + 0.3d + 0.4e)
           nothing
       end
foo (generic function with 1 method)

julia> function goo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}, e::Vector{T}) where T
           @assert length(a) == length(b) == length(c) == length(d) == length(e)
           @inbounds for i in eachindex(a)
               a[i] = b[i] + 0.1 * (0.2c[i] + 0.3d[i] + 0.4e[i])
           end
           nothing
       end
goo (generic function with 1 method)

julia> a,b,c,d,e=(rand(1000) for i in 1:5)
Base.Generator{UnitRange{Int64},getfield(Main, Symbol("##9#10"))}(getfield(Main, Symbol("##9#10"))(), 1:5)

julia> @btime foo($a,$b,$c,$d,$e)
  1.277 μs (0 allocations: 0 bytes)

julia> @btime goo($a,$b,$c,$d,$e)
  345.568 ns (0 allocations: 0 bytes)

julia> versioninfo()
Julia Version 0.7.0-beta2.12
Commit a878341 (2018-07-15 15:57 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
  JULIA_PKG3_PRECOMPILE = 1

[gasagna]

Maybe related to #27988 ?

[KristofferC]

There is no inference problem, the inbounds loop manages to simd while the broadcast doesn't. Removing the inbounds and they are equally fast.

[YingboMa]

No, I don't think so. After applying this diff and rebuilding Julia, I got the same result.

julia> @btime foo($a,$b,$c,$d,$e)
  1.274 μs (0 allocations: 0 bytes)

julia> @btime goo($a,$b,$c,$d,$e)
  348.897 ns (0 allocations: 0 bytes)

julia> versioninfo()
Julia Version 0.7.0-beta2.12
Commit a878341* (2018-07-15 15:57 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
Environment:
  JULIA_PKG3_PRECOMPILE = 1

[jebej]

Don't dot-ops do operations inbounds?

[YingboMa]

Maybe related to #22255?

julia> using BenchmarkTools

julia> a,b,c,d,e=(rand(1000) for i in 1:5)
Base.Generator{UnitRange{Int64},getfield(Main, Symbol("##5#6"))}(getfield(Main, Symbol("##5#6"))(), 1:5)

julia> function foo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}) where T
           @. a = b + 0.2c + 0.3d
           nothing
       end
foo (generic function with 1 method)

julia> function goo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}) where T
           @assert length(a) == length(b) == length(c) == length(d)
           @inbounds for i in eachindex(a)
               a[i] = b[i] + 0.2c[i] + 0.3d[i]
           end
           nothing
       end
goo (generic function with 1 method)

julia> function foo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}, e::Vector{T}) where T
           @. a = b + 0.2c + 0.3d + 0.4e
           nothing
       end
foo (generic function with 2 methods)

julia> function goo(a::Vector{T}, b::Vector{T}, c::Vector{T}, d::Vector{T}, e::Vector{T}) where T
           @assert length(a) == length(b) == length(c) == length(d) == length(e)
           @inbounds for i in eachindex(a)
               a[i] = b[i] + 0.2c[i] + 0.3d[i] + 0.4e[i]
           end
           nothing
       end
goo (generic function with 2 methods)

julia> @btime foo($a,$b,$c,$d)
  168.767 ns (0 allocations: 0 bytes)

julia> @btime goo($a,$b,$c,$d)
  166.587 ns (0 allocations: 0 bytes)

julia> @btime foo($a,$b,$c,$d,$e)
  1.135 μs (0 allocations: 0 bytes)

julia> @btime goo($a,$b,$c,$d,$e)
  342.142 ns (0 allocations: 0 bytes)

[Keno]

Looks like LLVM is having some trouble with the array bounds in the broadcast example:

remark: <unknown>:0:0: loop not vectorized: cannot identify array bounds

[Keno]

If I'm reading this correctly, the array pointer load is not getting LICM'ed. @vtjnash was looking at something related recently:

  %171 = load double*, double* addrspace(11)* %152, align 8
  %172 = getelementptr inbounds double, double* %171, i64 %value_phi1517.us.us
  store double %170, double* %172, align 8

[Keno]

Actually, my bad, I accidentally dropped TBAA.

[Keno]

Ok, the array bounds are still the issue though. Not sure why LLVM is unable to find them. We should probably fix that. However, there is a larger issue here. We go through all this trouble on the julia side to unalias the destination array from any of the arguments, but we don't actually tell LLVM that we did that, so it tries to insert a dynamic memory range check for which it then can't find the bounds. The better fix here is probably to figure out how to inform LLVM that the pointers don't alias.

[ChrisRackauckas]

the inbounds loop manages to simd while the broadcast doesn't

If this is the case, why? Why not assert that the sizes work and then just @inbounds @simd (or let it auto-vectorize)?

[ChrisRackauckas]

Maybe related to #22255?

It's not #22255 because if it was then the reduced size example would have no regression and the larger size example would have allocations.

[chriselrod]

If this is the case, why? Why not assert that the sizes work and then just @inbounds @simd (or let it auto-vectorize)?

I don't know the details of how broadcasting is implemented (and all the LLVM interfacing is way beyond me), but as Keno elaborated above, LLVM isn't aware that the array's aren't aliased while broadcasting, and therefore doesn't auto-vectorize (if b[i+1] refers to the same element as a[i] in @. a = b + 0.1 * (0.2c + 0.3d + 0.4e), you'd obviously have to calculate each a[i] in order, instead of in SIMD chunks).