`sizehint!` style approach for matrices

My program has the following structure:

for i in 1:1000
    u = foo(i) # A 3xN Matrix{Float64}
    bar(u)
end

where N is usually between 100 and 150.

Through profiling, I have learned that allocating memory in line 2 is a bottleneck for this program, and I would like to reduce these allocations. The size of foo(i) is not known until runtime, so I can’t use my usual approach of setting u = zeros(3, 150) and then u[:] = foo().

In cases where u is a vector, you can improve performance in this kind of situation with sizehint! and empty!: something like

u = Float64[]
sizehint!(u, 150)
for i in 1:1000
    append!(u, foo(i)) # An N Vector{Float64}
    bar(u)
    empty!(u)
    # sizehint!(u)
end

(I think the commented sizehint! is unnecessary per Documenting performance model of `empty!`, `sizehint!`, `push!` & friends .)

Can I do something similar for a matrix? Or is this wishful thinking because foo() will have to allocate its own memory anyway? (foo() is actually a call to solve() from DifferentialEquations.jl and I am doing some post analysis on the trajectory.)

How about preallocating a buffer for u that is the biggest possible size and then just using views or something? If N is really bounded by 150 or whatever, you could just allocate u = Array{Float64}(undef, 3, 150) or something and then for each loop do something like

for i in 1:1000
  N = get_N(i)
  uv = view(u, :, 1:N)
  [...]
end

EDIT: I just remembered that there is a package for resizing multi-index arrays like you’re asking. So I guess the most direct answer to your question is to use ElasticArrays.jl. But I do think preallocation a matrix and just using views in the loop is also something worth considering, depending on the details of your problem.

2 Likes

Rather than a 3xN Matrix{Float64}, consider using a Vector{SVector{3,Float64}} with the StaticArrays.jl package. I’ll speculate that (like me when I want this pattern) what you’re really doing is making a collection of 3-vectors anyway, rather than an actual “matrix”. You can sizehint!, push!, and empty! this Vector – which you can’t do with a Matrix.

If you need the matrix form at some point, you can reinterpret(reshape,Float64,vector_of_svectors) to get essentially the same thing you would have built using your original mechanism.

1 Like