Efficient operations on vectors of static things

General Usage Performance
,
1

Hello.

I am (still) trying to make sense on how to efficiently traverse a vector containing static vectors and matrices and doing operations on them. This is in the context of molecular dynamics, where I have N particles, with vectors on them (positions, velocities), and matrices (e.g. inertia tensors).

So, to compute the norm of all vectors I have the naïve implementation (function r2) , which is greatly sped-up with the simple use of an inner variable r_i = r[i] (function r2_v2). Using a regular vector for this (function r2_v3) also speeds things up, but not as greatly. All this, due to fewer allocations, I guess.

  174.579 μs (4490 allocations: 140.58 KiB)
  589.067 ns (1 allocation: 7.94 KiB)
  8.175 μs (2 allocations: 8.02 KiB)

But then, when I compute the operation r·J·r (which I need, long story…), these tricks do not work. The number of allocations does not go down. Why?

  32.308 μs (1490 allocations: 31.20 KiB)
  46.945 μs (1979 allocations: 38.84 KiB)
  704.404 μs (6493 allocations: 328.33 KiB)

Here is the code. Thanks for any hints whatsoever.

using StaticArrays

using LinearAlgebra

using BenchmarkTools

N = 1000 ;

rr  =  Vector{SVector{2}}(undef,N);

for i in 1:N rr[i] = rand(2) end

function r2(  r  ) 
    N = length(r)

    rr = Vector{Float64}(undef,N)

    for i in 1:N
        rr[i] = dot( r[i] , r[i] )
    end

    return rr

end

@btime r2($rr) ;

function r2_v2(  r  ) 
    N = length(r)

    rr = Vector{Float64}(undef,N)

    for i in 1:N
        r_i = rr[i]
        rr[i] = dot( r_i , r_i )
    end

    return rr

end

@btime r2_v2($rr) ;

function r2_v3(  r  ) 
    N = length(r)

    rr = Vector{Float64}(undef,N)
    r_i = zeros(2)
    
    for i in 1:N
        r_i .= rr[i]
        rr[i] = dot( r_i , r_i )
    end

    return rr

end

@btime r2_v3($rr) ;

JJ  =  Vector{SMatrix{2, 2}}(undef,N);

for i in 1:N JJ[i] = rand(2,2) end

function rJr(  r , J ) 
    N = length(r)

    rjr = Vector{Float64}(undef,N)

    for i in 1:N
        rjr[i] = dot( r[i] , J[i], r[i] )
    end

    return rjr

end

@btime rJr($rr , $JJ) ;

function rJr_v2(  r , J ) 
    N = length(r)

    rjr = Vector{Float64}(undef,N)

    for i in 1:N
        r_i = rr[i]
        J_i = J[i]

        rjr[i] = dot( r_i , J_i, r_i )
    end

    return rjr

end

@btime rJr_v2($rr , $JJ) ;

function rJr_v3(  r , J ) 
    N = length(r)

    rjr = Vector{Float64}(undef,N)

    r_i = zeros(2)
    J_i = zeros(2,2)
    Jr_i = zeros(2)

    for i in 1:N
        r_i .= rr[i]
        J_i .= J[i]
        Jr_i.= J_i * r_i
        rjr[i] = dot( r_i , Jr_i )
    end

    return rjr

end

@btime rJr_v3($rr , $JJ) ;
3

This is an abstract type. You need SMatrix{2,2,Float64,4}, the 4 being redundant but unfortunately necessary.

3 Likes
4

Then this is

JJ = rand(SMatrix{2,2,Float64,4}, 1000)
3 Likes
5

Thank you. Yes, I had been fighting with an elegant initialization for JJ, and had to settle for a loop :smiley:

Anyway, your suggestion only seems to help with version 3 of the r·J·r code:

  36.568 μs (2490 allocations: 78.08 KiB)
  50.403 μs (2979 allocations: 85.72 KiB)
  141.808 μs (2493 allocations: 125.20 KiB)

It seems also strange that SVector should work, but not SMatrix, but these things do happen.

6

This is also abstract and should be SVector{2,Float64}

2 Likes
7

I’m on the phone, but the only allocations you should see on those loops, types fixed to be concrete, are in this line of the last example, because the product on the right will allocate a new array for these non-static vectors and matrices (to eliminate that one you would need do use mul!).

1 Like
8

Hi, and thanks again for the two suggestions. The suggestion about turning SVector{2} into SVector{2,Float64} has a great effect:

  971.286 ns (1 allocation: 7.94 KiB)
  558.457 ns (1 allocation: 7.94 KiB)
  8.301 μs (2 allocations: 8.02 KiB)
  4.077 μs (1 allocation: 7.94 KiB)
  54.752 μs (3979 allocations: 116.97 KiB)
  142.932 μs (3493 allocations: 156.45 KiB)

Even tough the numbers are exciting, it is still not clear to me that the function with the inner r_i, r2_v2(r) is slightly better, and does not allocate — while the
one with inner r_i and J_i, rJr_v2(r) is allocating like crazy.

You can see rJr_v3(r) is also not much better.

9

Be careful, all your v2 and v3 versions of the code dont do what you want to. Make sure when optimizing that you dont change the end result.

In particular, in your inner loops, you do r_i .= rr[i], where you should do r_i .= r[i]. I think that should solve the weird performance results (if you also fix the type of the static vector as pointed out above).

2 Likes
10

Oh, man, I can’t believe those globals have been messing with my functions again! The “r2” functions where ok, because the “rr” inside them overrode the global one. Not so for the rJr, where “rr” is a total bug.

Thank you very much, now results make complete sense

  668.019 ns (1 allocation: 7.94 KiB)
  594.266 ns (1 allocation: 7.94 KiB)
  8.133 μs (2 allocations: 8.02 KiB)
  3.779 μs (1 allocation: 7.94 KiB)
  3.703 μs (1 allocation: 7.94 KiB)
  76.301 μs (1004 allocations: 86.31 KiB)

And, seriously, is there a way to stop global variables from entering functions? I have been appending “gl_” to all my global variables in my code, but sometimes (as in this example) I forget.

11

Yes there are ways. Either you

  1. create a main() function and run that so you dont have any globals
  2. avoid using the same variable names in your functions as existing globals (avoid name shadowing). Perhaps you can use some naming convention like putting globals in all capitals or appending _global to the name (as you are doing)
  3. Sometimes cleaning up your globals by putting them in structs or named tuples can also help to limit the surface for these types of issues.

But (as far as i know) there is no macro or syntax to prohibit a function from using any (untyped) global variables.

1 Like
12

Thank you for your suggestions, and thanks to all for helping. Marking this one as solved. Sorry, I can only mark one answer as “solution”, but the global variable stuff is really my fault, and a bug.

13

These implementations are wrong.

14

Has been repeatedly bitten by globals forgotten in the body of a function. One of the ways of preventing it is (temporarily) putting your function into a let block:

foo = let 
    function foo(x)
        return x
    end
end

But actually, after inflicting myself more pain that was necessary, I came to conclusion that putting all my functions into a package from the very beginning, while keeping globals outside of it, solves this problem, too.

2 Likes