Right way to "static" allocate vectors in a function


#1

Hi guys!

I have some functions that were called to fill a 2000 x 2000 matrix. Some functions need to allocate 2000x1 vectors to do the computations. However, those functions are called like 2000x2000 times :sweat_smile: Something like:

function func1(...)
    den_ar = Vector{Float64}(undef,M)
    int_fe     = Vector{Float64}(undef,M)
    ablacao_fe = Vector{Float64}(undef,M)
    ...

The problem is that the allocations are very, very high. My measurements are showing 145 GB in allocations. The execution time is 135s vs. 46s in FORTRAN.

In C, this is solved by making those vectors static. I read some posts here and even a discussion in GitHub about a static variable in a Julia function. However, I was not be able to find what is the best approach to it.

Can anyone help me please?


#2

First of all, I’d suggest not doing this at all. Instead, create a Func1Workspace type (or whatever you want to call it) which holds those vectors. Then change your method to something like:

function func1!(work::Func1Workspace, ...)
  ...
end

You can even make a convenient definition that constructs a workspace (for situations where convenience is more important than performance):

function func1(...)
  work = Func1Workspace()
  func1!(work, ...)
end

In this way, you gain a number of benefits:

  1. Your code is safe for parallel execution (while the static approach definitely is not)
  2. It’s easy to handle workspaces with different sizes, which the static approach can’t easily do.
  3. It’s very obvious what is going on and when allocation is happening, which makes things easy for users of your code to understand

That said, if you really want static variables, you can use a let block and an anonymous function:

julia> const foo = let x = zeros(100)
         function (y)
           x .= y
           x
         end
       end
#3 (generic function with 1 method)

julia> foo(1)
100-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 ⋮  
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

julia> @allocated foo(2)
0

#3

Thanks @rdeits ! I will go with the workspace approach. Seems way cleaner.