About memory allocation of auxiliary array in successive calls of a function

General Usage Performance
1

If I call a function that internally requires the allocation of an array for some computation, and afterwards I call this same function with the exact same type of input variables, is the auxiliary array which was allocated the first time allocated again?

Example:

x = rand(3)
function h!(x)
   A = rand(3,3)
   x = A*x
end
# first call (function is compiled and certainly A is allocated)
h!(x)
# second call (function is not compiled, but is A reallocated?)
h!(x)

My question is because I am not sure if, to improve performance, it is a good idea to pass A as well, with:

x = rand(3)
A = Matrix{Float64}(undef,3,3)
function h!(x,A)
   @. A = rand()
   x = A*x
end

Or if the compiler/running-time does that already, allowing a cleaner code. From a few tests I performed here it seems that the second option is needed for maximum performance, but I would like to hear from someone that really knows what is under the hood.

(ps: I am coming from Fortran, in which everything is allocated on execution, such that the declaration of auxiliary arrays inside a function certainly does not lead to multiple allocations, and the explicit passing of the the auxiliary arrays which are only used internally and have fixed sizes is not needed).

Thank you.

2

A is allocated every time you call h!(x)

You can do this

function h!(x,A)
          Random.rand!(A)
          mul!(x,A,x);
       end
3 Likes
3

Great. But ‘x=A*x’ doesn’t allocate anything, wright? Using mult! there makes any difference?

4

x = A*x allocates a vector for x.

2 Likes
5

Uhm… thanks. That independence of names from actual values still confuses me from time to time. That is the same as ‘y=A*x’, of course.

Edit: added this comment just to remember how things work, and perhaps helps other newbies:

julia> function f(x)
           x = 2*x
           return x
       end
f (generic function with 1 method)

julia> x = [ 1, 2 ]
2-element Array{Int64,1}:
 1
 2

julia> y = f(x)
2-element Array{Int64,1}:
 2
 4

julia> x
2-element Array{Int64,1}:
 1
 2

Note that the x which enters f(x) is not modified.

1 Like
6

I don’t think mul! can be used as suggested without using a 2nd vector. If you look at the docs, it says that the 1st and 3rd arguments can not be aliases of each other, which appears to be correct – when I try it, I get a zero vector as a result.

Julia does have a solution for “static” arrays so that you can call a function multiple times without having to reallocate or explicitly pass an array. I’d prefer Julia had a cleaner static syntax, but this works:

using LinearAlgebra
using Random
let y = zeros(3), A = zeros(3,3)
  global function h!(x)
    rand!(A)
    mul!(y,A,x)
    x .= y
  end
end

julia> x = rand(3)
3-element Array{Float64,1}:
 0.528643672437896
 0.6445978634354903
 0.13390924294137552

julia> @time h!(x)   # notice zero allocations reported
  0.000005 seconds
3-element Array{Float64,1}:
 0.9482372634076288
 0.9135226035098386
 0.6381315118219281

julia> x
3-element Array{Float64,1}:
 0.9482372634076288
 0.9135226035098386
 0.6381315118219281
7

If A is of fixed and small dimension, you can use StaticArrays to move everything from the heap onto the stack:

julia> function h!(x)
          A = @SMatrix rand(3,3)
          x .= A*x
       end
       x = @MVector rand(3)
       @time h!(x)  # This compiles and hence allocates
       @time h!(x)  # Look Ma, no allocations!
  0.010066 seconds (23.40 k allocations: 977.899 KiB)
  0.000004 seconds
3-element MArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 0.3388153014844782
 0.9433372130533029
 0.8048580079870412
2 Likes