Prevent array re-allocation inside a loop (I am not allowed to allocate before the loop)

General Usage Performance
1

Hi,
I am executing a code similar to the simple one below under the constraint that I cannot allocate anything before the loop due to some constraints of the calling function.

for it=1:ntime
     my_fun()
end

where

my_fun()    
    F           = zeros(T, npoin, nelem)
    ...

As you can see, this reallocates memory at every loop iteration. I would like to avoid this and only allocate it the first time and reset the array to zero the following times.

The problem would not exist if I could allocate F outside of the loop, but I cannot do that.

Is the a way to achieve this?
thanks

2

Can you modify my_fun()?

If so, even if not very good for practice, you can define in the global scope:

const global_F = zeros(T, npoin, nelem)

and then inside my_fun() do:

function my_fun()
    F = fill!(global_F, zero(T))

but that comes with a lot of possible issues to care about (maintainability, races, etc.).

3

I can modify my_fun() and everything else before it’s called, but nothing out of the loop.
I can try this.
Feel free to add more options if anything else comes to mind. I don’t particularly like the const option of a global variable.

thanks @lmiq

4

Can you explain a little bit this restriction?

(it would prevent even the definition of a global buffer, as I suggested - and complicates things quite a bit, since the variables created in each loop iteration are local to that iteration, so they woudn’t survive for the next iteration to be reused).

One alternative, make my_fun a callable struct that carries its buffer:

struct MyFun{T}
    buffer::Matrix{T}
end

function (my_fun::MyFun)(args)
    F = fill!(my_fun.buffer, 0.0)
    ...
end

const my_fun = MyFun()

That is a little bit more clean that the global buffer, as at least the buffer is carried with the struct, but still has to be declared in the global scope (otherwise if you instantiate a new MyFun object in each iteration of the loop, the buffer will be recreated the same).

1 Like
5 
for it=1:ntime
    If it == 1
        X = ... # allocate
    end
    my_fun(X)
end

?

1 Like
6

Doesn’t work, the x is local to the loop iteration:

julia> function f()
           for i in 1:2
               if i == 1
                   x = zeros(2)
               end
               if i == 2
                   println(" iteration 2 ")
                   x .= 0.0
               end
           end
           return x
       end
f (generic function with 1 method)

julia> f()
 iteration 2 
ERROR: UndefVarError: `x` not defined
Stacktrace:
 [1] f()
   @ Main ./REPL[3]:8
 [2] top-level scope
   @ REPL[4]:1

One would need to add something like local x outside the loop, but that’s forbidden by the rules :slight_smile: .

7

This is somewhat standard:

let F = zeros(T, npoin, nelem)
    global function my_fun()
        F .= zero(T)
        ...
    end
end

but I would rather go for a callable struct.

8

@Imiq, the issue with the calling function is the following:
I am using DifferentialEquations.jl to solve du/dt = RHS. This packge is used as defined below, where u is the solution array and rhs! is the function that calculates the RHS vector at every time step.

function rhs!() is the only function that I can control when it comes to allocations and array definitions etc. So, if rhs!() allocates anything, it does so at every time iteration
of solve.


 prob = ODEProblem(rhs!,
                      u,
                      tspan,
                      params);
    
solution = solve(prob, solver, dt=Δt...);

I cannot pass other than function rhs!(), u, tspan, params to ODEProblem().

1 Like
9

You can preallocate everything in Params.
Note that you can provide any type you like there so an array or a tuple or your own struct are all fine

2 Likes
10

You can pass in pre-allocated arrays via params etcetera.

See also the PreallocationTools.jl package, which is specifically designed for use with ODEProblem and similar.

1 Like
11

Thank you for clarifying this. PreallocationTools.jl does exactly what I need!