Static Arrays Memory & Allocations

Below two functions

function generateLocalMatrixWorks(element::Element)
h = element.area
Iloc = SVector(element.e1, element.e1, element.e2, element.e2)
Jloc = SVector(element.e1, element.e2, element.e1, element.e2)
Aloc = SVector(1/h, -1/h, -1/h, 1/h)
return Iloc, Jloc, Aloc
end

function generateLocalMatrix(element::Element)
h = element.area
aux1d = SVector(element.e1, element.e2)
Iloc = SVector{4}(repeat(aux1d,2,1))
Jloc = Iloc;
Aloc = SVector(1/h, -1/h, -1/h, 1/h)
return Iloc, Jloc, Aloc
end

The first does not allocate. The second one does. I fail to understand what I am doing wrong. Thanks.

Edit: this seems to be related to having to create static matrices in one go as referred to, e.g., here Replacing normal arrays with static arrays - #2 by foobar_lv2

Please format your code so we could read it easily. See here for more: Please read: make it easier to help you

1 Like

Probably the repeat call is returning a vector, and thus allocating.

Thanks. I agree that the repeat function is likely the bad guy. I verified that comprehension does allocate as well. What alternative do you see?

Indeed repeat returns a regular array. You could do

Iloc = [aux1d; aux1d]

To get an SVector. The bracket syntax has a specialization for SVectors.

1 Like

FYI I prefer fill over repeat in 99% of cases as it does not allocate. IIRC it doesn’t create an intermediate array, but may not be appropriate in this case.

This is an elegant answer!

Thanks! How do I find the other operations that have a specialization for SVectors?

I rely on experimentation. Just check the type of the output for various functions. Sometimes it matters how you call certain functions, for example map(x->2*x, v) returns an SVector if v is an SVector but map(i->2*v[i], 1:length(v)) will return a Vector (not a surprise if you think about it). If you are curious, you can also check what code is being called with @less, for example try @less [aux1d; aux1d].

It probably returns a matrix. To get a vector it should be repeat(aux1d, 2) instead.

@Domenico_Lahaye, here you don’t want an Array at all, but in general you should be aware that Julia (unlike e.g. Matlab) distinguishes vectors from Nx1 matrices. Therefore, you use zeros(N), ones(N), rand(N) and repeat(x, N) instead of zeros(N, 1), ones(N,1), rand(N, 1) and repeat(x, N, 1), etc.

@tverho Thx! How should I make sense of the output that @less generate?

@DNF Thx! Clear. Matlab habbits are sometimes hard to leave behind.

I would like to add local variables to my function. The first variant works fine. The second variant fails. I fail to understand this behavior.

function generateLocalMatrix(element::Element)
    h     = element.area 
    e1    = element.e1
    e2    = element.e2
    Iloc  = SVector(e1, e1, e2, e2) 
    Jloc  = SVector(e1, e2, e1, e2) 
    # Kloc: local variable to be used in the future: goes fine - *no* memory allocations
    Kloc  = SVector(1., 2., 3., 4.) 
    Aloc  = SVector(1/h, -1/h, -1/h, 1/h) 
    return Iloc, Jloc, Aloc
end
function generateLocalMatrix(element::Element)
    h     = element.area 
    e1    = element.e1
    e2    = element.e2
    Iloc  = SVector(e1, e1, e2, e2) 
    Jloc  = SVector(e1, e2, e1, e2) 
    # Kloc: local variable to be used in the future: goes fine - *no* memory allocations
    Kloc  = SVector(1., 2., 3., 4.)
    # Lloc: another local variable to be used in the future: fials - causes memory allocations
    Lloc  = SMatrix{2,2}([1. 2.;3. 4.])  
    Aloc  = SVector(1/h, -1/h, -1/h, 1/h) 
    return Iloc, Jloc, Aloc
end

This defined a matrix, explicitly, and causes the allocation.

Try

SMatrix{2,2}(1.0, 3.0, 2.0, 4.0)

or

@SMatrix [1.0 2.0; 3.0 4.0]

or

SA[1.0 2.0; 3.0 4.0] 

Thx again! Amazing stuff!

@less shows the source file where the function is defined, at the function definition. It’s maybe more clear if you use it with a regular function. [aux1d; aux1d] is translated by Julia directly into vcat(aux1d, aux1d) so that’s why you see the definition of (one method for) vcat.