Unexpected allocations in loops

General Usage Performance
1

Consider the following reduced example:

using BenchmarkTools

const F = Float64

struct Species
  q::F
  ρ::Array{F}
end

addCharge!(dst::Array{F}, q::F, ρ::Array{F}) =
  for I ∈ eachindex(dst) dst[I] += q * ρ[I] end

addCharges_bad!(dst::Array{F}, species::Vector{Species}) = begin
  fill!(dst, 0)
  for s ∈ species
    for I ∈ eachindex(dst) dst[I] += s.q * s.ρ[I] end
  end
end

addCharges_good!(dst::Array{F}, species::Vector{Species}) = begin
  fill!(dst, 0)
  for s ∈ species
    addCharge!(dst, s.q, s.ρ)
  end
end

main() = begin
  a = rand(100, 100, 100)
  species = [Species(1., rand(size(a)...)), Species(1., rand(size(a)...))]

  @btime addCharges_bad!($a, $species)
  @btime addCharges_good!($a, $species)

  return
end

main()

Results:

  84.034 ms (3998978 allocations: 61.02 MiB)
  2.804 ms (2 allocations: 32 bytes)

I expected addCharges_bad! not to allocate, but it does badly.
There is no reason to allocate here since no aliasing occurs.
Instead, I have to hoist the inner for loop to a function addCharge! to perform as expected.

Is there a mechanism in Julia allowing to mark dst with something like @noalias, thus ensuring no intermediate copies ? Or is there something I have missed ?

1 Like
2 
Array{F}

Is an abstract type. Probably you want just Vector{F}.

Edit: actually since you want to define the dimension at runtime, you may need to define:

Such that you get a concrete species.

2 Likes
3

Thanks @lmiq.

This works:

struct Species{N}
  q::F
  ρ::Array{F, N}
end

since isconcretype(Array{Float64, 3}) == true.

Since this part comes from a much more complex code (and took some time to figure out), are you aware of any tool in Julia itself or an external package that could warn against these kind of mistakes ? Some kind of annotation maybe ?

4

Reduced case:

const F = Float64

struct Species{N}
  q::F
  ρ::Array{F, N}
end

struct BadSpecies
  q::F
  ρ::Array{F}
end

main() = begin
  a = rand(10, 10, 10)
  @code_llvm Species(1., rand(size(a)...))
  @code_llvm BadSpecies(1., rand(size(a)...))
end

main()

Results in:

;  @ [...].jl:7 within `Species`
define void @julia_Species_374({ double, {}* }* noalias nocapture sret({ double, {}* }) %0, [1 x {}*]* noalias nocapture %1, {}* nonnull readonly %2, double %3, {}* nonnull align 16 dereferenceable(40) %4) #0 {
top:
  %5 = getelementptr inbounds [1 x {}*], [1 x {}*]* %1, i64 0, i64 0
  store {}* %4, {}** %5, align 8
  %.repack = getelementptr inbounds { double, {}* }, { double, {}* }* %0, i64 0, i32 0
  store double %3, double* %.repack, align 8
  %.repack1 = getelementptr inbounds { double, {}* }, { double, {}* }* %0, i64 0, i32 1
  store {}* %4, {}** %.repack1, align 8
  ret void
}
;  @ [...].jl:12 within `BadSpecies`
define void @julia_BadSpecies_388({ double, {}* }* noalias nocapture sret({ double, {}* }) %0, [1 x {}*]* noalias nocapture %1, double %2, {}* nonnull align 16 dereferenceable(40) %3) #0 {
top:
  %4 = getelementptr inbounds [1 x {}*], [1 x {}*]* %1, i64 0, i64 0
  store {}* %3, {}** %4, align 8
  %.repack = getelementptr inbounds { double, {}* }, { double, {}* }* %0, i64 0, i32 0
  store double %2, double* %.repack, align 8
  %.repack1 = getelementptr inbounds { double, {}* }, { double, {}* }* %0, i64 0, i32 1
  store {}* %3, {}** %.repack1, align 8
  ret void
}

There is an extra {}* nonnull readonly %2 argument in the Species emitted code.

5

I’m not at the computer, but @code_warnype probably warns you about something. There is also JET.jl (never used it, but people talk positively about it here) and Chutlhu.jl.

At the same time, with some practice one starts to see these things, as they are central for the Julia programing model.

1 Like
6

Thanks for your help.
Indeed @code_warntype hints (in red) with:

│   %20 = Base.getproperty(s, :ρ)::Array{Float64}
versus
│   %20 = Base.getproperty(s, :ρ)::Array{Float64, 3}

JET does not report errors.

1 Like
7

FWIW, there is also the good old https://github.com/JunoLab/Traceur.jl which should / might warn about this.

1 Like
8

Thanks @carstenbauer, this really helps: this is the first tool that reports something by analysing the whole program, and not just a specific call (I wish @code_warntype could be recursive ? : @code_warntype main2() does not report anything).

→ Cthulhu does it: @descend_code_warntype main2()

using Traceur

const F = Float64

struct Species{N}
  q::F
  ρ::Array{F, N}
end

struct BadSpecies
  q::F
  ρ::Array{F}
end

addCharges!(dst::Array, species::Vector) = begin
  fill!(dst, 0)
  for s ∈ species
    for I ∈ eachindex(dst) dst[I] += s.q * s.ρ[I] end
  end
end

main2() = begin
  a = ones(10, 10, 10)
  addCharges!(a, [BadSpecies(1., zeros(size(a)...))])
end

main3() = begin
  a = ones(10, 10, 10)
  addCharges!(a, [Species(1., zeros(size(a)...))])
end

julia> @trace main2()
[...]
┌ Warning:  is assigned as Union{Nothing, Tuple{BadSpecies, Int64}}
┌ Warning: dynamic dispatch to Base.getindex(Base.getfield(φ (%48 => %44, %103 => %98), ρ), φ (%63 => %60, %84 => %79))
└ @ [...].jl:24

julia> @trace main3()
[...]
┌ Warning:  is assigned as Union{Nothing, Tuple{Species{3}, Int64}}

I had to change rand() calls since I hit github.com/JunoLab/Traceur.jl/issues/35 though.

The trap I’ve fallen into.

1 Like
9

I’m completing my answer with the most automated way (julia >= 1.7):

julia> using JET
julia> @report_opt main2()
═════ 1 possible error found ═════
┌ @ [...]/incomplete_allocs.jl:70 Main.addCharges!(a, Base.vect(Main.BadSpecies(1.0, Main.zeros(Main.size(a)...))))
│┌ @ [...]/incomplete_allocs.jl:64 Base.getindex(%69, %65)
││ runtime dispatch detected: Base.getindex(%69::Array{Float64}, %65::Int64)
│└──────────────────────────────────────────────────────────────