Strategies to reuse memory

An interesting pattern for this is to define the function with keyword parameters defining the auxiliary arrays (I will use a tuple here as suggested above):

julia> function f!(x;aux=(zeros(3),zeros(3)))
           a, b = aux
           a .= 2 .* x
           b .= x .^ 2
           x .= x .+ a .* b
           return x
       end
f (generic function with 1 method)

julia> x = rand(3);

julia> @btime f!(x) setup=(x=copy(xin)) evals=1
  135.000 ns (2 allocations: 160 bytes)
3-element Vector{Float64}:
 0.9603701960196425
 1.0196162985304351
 2.727934952974922

julia> auxin = (zeros(3),zeros(3));

julia> @btime f!(x,aux=aux) setup=(x=copy(xin);aux=deepcopy(auxin)) evals=1
  66.000 ns (0 allocations: 0 bytes)
3-element Vector{Float64}:
 0.9603701960196425
 1.0196162985304351
 2.727934952974922

The nice thing of this pattern is that you can develop the code without worrying about these allocations and then add them as a performance optimization where needed. Thus the complexity appears only where really needed.

It is possible, of course, to create constant global variables to act as a buffer, but that will break modularity and at the end increase code maintenance complexity.

Concerning this point, I think one useful pattern is to define two functions, taking the advantage of multiple dispatch, one that receives preallocated vectors and return them, and the other that does not receive them and does not return them. Then you can again only take of care preallocations when needed, as an optimization:

julia> function f!(x,aux)
           a, b = aux
           x .= x .+ a .* b
           return x, aux
       end
f! (generic function with 2 methods)

julia> function f!(x)
           aux = (zeros(3),zeros(3))
           f!(x,aux)
           return x
       end
f! (generic function with 2 methods)

julia> x = rand(3); aux = (zeros(3),zeros(3))
([0.0, 0.0, 0.0], [0.0, 0.0, 0.0])

julia> x, aux = f!(x,aux)
([0.01274472197573262, 0.9823930168953983, 0.48667397158434544], ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0]))

julia> x = f!(x)
3-element Vector{Float64}:
 0.01274472197573262
 0.9823930168953983
 0.48667397158434544

(ps: since I wrote everything mutating in place, returning the variables is optional here, in these examples the functions could well return nothing instead)

7 Likes