Why does this matrix increment allocate?

New to Julia
1 
module somemodule

using BenchmarkTools

struct MyType
    A::AbstractMatrix
    B::AbstractMatrix
end

myinc(A::AbstractMatrix, B::AbstractMatrix) = A .+= B
mytypeinc(t::MyType, B::AbstractMatrix) = t.A .+= B
mytypeinc2(t::MyType, B::AbstractMatrix) = myinc(t.A, B)

function perf()
    H = [1.0 0 ; 1.0 -1.0];
    Q = [0.01 0 ; 0 0.01];
    t = MyType(H, Q);

    @btime myinc($H, $Q)
    @btime mytypeinc($t, $Q)
    @btime mytypeinc2($t, $Q)
end

perf()

end

This outputs

  20.435 ns (0 allocations: 0 bytes)
  445.258 ns (2 allocations: 64 bytes)
  28.728 ns (0 allocations: 0 bytes)

Initially I had implemented in my code something like mytypeinc, I was surprised to see that it allocates, eventually tried myinc and that doesn’t allocate. So I tried to write an “improved” version of mytypeinc called mytypeinc2 which just calls myinc, and that works to avoid the allocation. But why is it that such a workaround is necessary? Why doesn’t mytypeinc do the obvious right thing? Clearly my mental model is missing something crucial.

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4570 CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, haswell)
Environment:
  JULIA_EDITOR = "/home/me/.vscode-server/bin/7db1a2b88f7557e0a43fec75b6ba7e50b3e9f77e/node"
  JULIA_NUM_THREADS =
2

Probably because these are abstract containers. Use Matrix{T}, or a full parametric type.

1 Like
3

That is a function barrier, which allows the inner method to specialize

4 
module somemodule

using BenchmarkTools

struct MyType{T}
    A::Matrix{T}
    B::Matrix{T}
end

myinc(A, B) = A .+= B
mytypeinc(t::MyType, B) = t.A .+= B
mytypeinc2(t::MyType, B) = myinc(t.A, B)

function perf()
    H = [1.0 0 ; 1.0 -1.0];
    Q = [0.01 0 ; 0 0.01];
    t = MyType{Float64}(H, Q);

    @btime myinc($H, $Q)
    @btime mytypeinc($t, $Q)
    @btime mytypeinc2($t, $Q)
end

perf()

end

  20.553 ns (0 allocations: 0 bytes)
  13.706 ns (0 allocations: 0 bytes)
  15.100 ns (0 allocations: 0 bytes)
1 Like
5

Thanks for the info, but what lesson should I take from here? (other than that I don’t understand abstract types and need to study them). That I should not make fields of structs abstract if I can avoid it?

For more context, my use case is very simple: I want to have a type which is a container for vectors and matrices, and I want to define on that type a specific pattern of linear algebra operations on those vectors and matrices, Given is, should I just make the change above to my type?

1 Like
6

You could also write:

struct MyType{T<:AbstractMatrix}
    A::T
    B::T
end
2 Likes
7

Yes, that is one important performance tip.

8

See docs

https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-fields-with-abstract-type

https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-fields-with-abstract-containers

3 Likes