# Efficient operations on vectors of static things
**URL:** https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387
**Category:** Performance
**Tags:** linearalgebra, simulations
**Created:** [June 9, 2024, 9:19am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387 "2024-06-09T09:19:00Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@ddcampayo](https://discourse.julialang.org/u/ddcampayo)
#### Post date: [June 9, 2024, 9:19am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/1 "2024-06-09T09:19:00Z")
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.
```julia
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?
```julia
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.
```julia
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) ;
```
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [June 9, 2024, 9:44am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/3 "2024-06-09T09:44:31Z")
> [@ddcampayo](#):
>
> `Vector{SMatrix{2, 2}}`
This is an abstract type. You need `SMatrix{2,2,Float64,4}`, the 4 being redundant but unfortunately necessary.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [June 9, 2024, 9:51am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/4 "2024-06-09T09:51:43Z")
> [@ddcampayo](#):
>
> `for i in 1:N JJ[i] = rand(2,2) end`
Then this is
```julia
JJ = rand(SMatrix{2,2,Float64,4}, 1000)
```
---
### Author:  [@ddcampayo](https://discourse.julialang.org/u/ddcampayo)
#### Post date: [June 9, 2024, 5:13pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/5 "2024-06-09T17:13:31Z")
Thank you. Yes, I had been fighting with an elegant initialization for JJ, and had to settle for a loop 😃
Anyway, your suggestion only seems to help with version 3 of the r·J·r code:
```julia
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.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [June 9, 2024, 5:34pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/6 "2024-06-09T17:34:39Z")
> [@ddcampayo](#):
>
> `SVector{2}`
This is also abstract and should be `SVector{2,Float64}`
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [June 9, 2024, 5:39pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/7 "2024-06-09T17:39:25Z")
> [@ddcampayo](#):
>
> `Jr_i.= J_i * r_i`
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!`).
---
### Author:  [@ddcampayo](https://discourse.julialang.org/u/ddcampayo)
#### Post date: [June 9, 2024, 6:02pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/8 "2024-06-09T18:02:37Z")
Hi, and thanks again for the two suggestions. The suggestion about turning `SVector{2}` into `SVector{2,Float64}` has a great effect:
```julia
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.
---
### Author:  [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)
#### Post date: [June 9, 2024, 6:07pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/9 "2024-06-09T18:07:01Z")
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).
---
### Author:  [@ddcampayo](https://discourse.julialang.org/u/ddcampayo)
#### Post date: [June 9, 2024, 7:20pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/10 "2024-06-09T19:20:43Z")
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
```julia
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.
---
### Author:  [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)
#### Post date: [June 9, 2024, 8:03pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/11 "2024-06-09T20:03:20Z")
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.
---
### Author:  [@ddcampayo](https://discourse.julialang.org/u/ddcampayo)
#### Post date: [June 9, 2024, 9:42pm UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/12 "2024-06-09T21:42:41Z")
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.
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [June 10, 2024, 12:18am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/13 "2024-06-10T00:18:17Z")
> [@ddcampayo](#):
>
> ```julia
> 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]
>
> ```
These implementations are wrong.
---
### Author:  [@Eben60](https://discourse.julialang.org/u/Eben60)
#### Post date: [June 10, 2024, 11:27am UTC](https://discourse.julialang.org/t/efficient-operations-on-vectors-of-static-things/115387/14 "2024-06-10T11:27:36Z")
> [@IlianPihlajamaa](#):
>
> But (as far as i know) there is no macro or syntax to prohibit a function from using any (untyped) global variables.
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:
```julia
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.