Design of N dimensional structure for Particle

General Usage
2

  1. I think those are exactly the same, because the representation in memory of a static vector of one dimension will be a number. Edit: not exactly true (see below), the lowered codes still have one call to getindex in addition in comparison to the version with Vector{Float64} (I thought that the compiler would inline that).

  2. The performance of a particle simulation is dependent on the performance of the algorithm computing the interactions, generally involving a cutoff and using cell lists to avoid looping over particles that are far from each other. That makes the access to the vector of positions discontinuous in memory anyway from particle to particle. Because of that, the best strategy is to store a vector of static vectors.

(If the access was always continuous, a matrix would be the best strategy to take the most of simd operations, and that is the case of the update of positions given the forces, but the costly part will be the calculation of the forces).

Adding to the answer 1. Check this example:

julia> using StaticArrays

julia> x = rand(1000);

julia> y = [ SVector{1,Float64}(x[i]) for i in 1:1000 ];

julia> function mysum1(x) # with Vector{Float64}
         s = 0.
         @inbounds @simd for i in 1:length(x)
           s += sin(x[i])
         end
         s
       end
mysum1 (generic function with 1 method)

julia> function mysum2(x) # with Vector{SVector{1,Float64}}
         s = 0.
         @inbounds @simd for i in 1:length(x)
           s += sin(x[i][1])
         end
         s
       end
mysum2 (generic function with 1 method)

Functions mysum1 and mysum2 get lowered to the following codes, which differ by:

mysum1:

│   %28 = s
│   %29 = Base.getindex(x, i)
│   %30 = Main.sin(%29)

mysum2:

│   %28 = s
│   %29 = Base.getindex(x, i)
│   %30 = Base.getindex(%29, 1)
│   %31 = Main.sin(%30)
Details 
julia> @code_lowered mysum1(x)
CodeInfo(
1 ─       Core.NewvarNode(:(val))
│         s = 0.0
│         $(Expr(:inbounds, true))
│   %4  = Main.length(x)
│         r#273 = 1:%4
│   %6  = Base.simd_outer_range
│   %7  = (%6)(r#273)
│         @_5 = Base.iterate(%7)
│   %9  = @_5 === nothing
│   %10 = Base.not_int(%9)
└──       goto #8 if not %10
2 ┄ %12 = @_5
│         i#274 = Core.getfield(%12, 1)
│   %14 = Core.getfield(%12, 2)
│   %15 = Base.simd_inner_length
│   %16 = r#273
│         n#275 = (%15)(%16, i#274)
│   %18 = Main.zero(n#275)
│   %19 = %18 < n#275
└──       goto #6 if not %19
3 ─       i#276 = Main.zero(n#275)
4 ┄ %22 = i#276 < n#275
└──       goto #6 if not %22
5 ─ %24 = Base.simd_index
│   %25 = r#273
│   %26 = i#274
│         i = (%24)(%25, %26, i#276)
│   %28 = s
│   %29 = Base.getindex(x, i)
│   %30 = Main.sin(%29)
│         s = %28 + %30
│         i#276 = i#276 + 1
│         $(Expr(:loopinfo, Symbol("julia.simdloop"), nothing))
└──       goto #4
6 ┄       @_5 = Base.iterate(%7, %14)
│   %36 = @_5 === nothing
│   %37 = Base.not_int(%36)
└──       goto #8 if not %37
7 ─       goto #2
8 ┄       val = Main.nothing
│         $(Expr(:inbounds, :pop))
│         val
└──       return s
)

julia> @code_lowered mysum2(y)
CodeInfo(
1 ─       Core.NewvarNode(:(val))
│         s = 0.0
│         $(Expr(:inbounds, true))
│   %4  = Main.length(x)
│         r#277 = 1:%4
│   %6  = Base.simd_outer_range
│   %7  = (%6)(r#277)
│         @_5 = Base.iterate(%7)
│   %9  = @_5 === nothing
│   %10 = Base.not_int(%9)
└──       goto #8 if not %10
2 ┄ %12 = @_5
│         i#278 = Core.getfield(%12, 1)
│   %14 = Core.getfield(%12, 2)
│   %15 = Base.simd_inner_length
│   %16 = r#277
│         n#279 = (%15)(%16, i#278)
│   %18 = Main.zero(n#279)
│   %19 = %18 < n#279
└──       goto #6 if not %19
3 ─       i#280 = Main.zero(n#279)
4 ┄ %22 = i#280 < n#279
└──       goto #6 if not %22
5 ─ %24 = Base.simd_index
│   %25 = r#277
│   %26 = i#278
│         i = (%24)(%25, %26, i#280)
│   %28 = s
│   %29 = Base.getindex(x, i)
│   %30 = Base.getindex(%29, 1)
│   %31 = Main.sin(%30)
│         s = %28 + %31
│         i#280 = i#280 + 1
│         $(Expr(:loopinfo, Symbol("julia.simdloop"), nothing))
└──       goto #4
6 ┄       @_5 = Base.iterate(%7, %14)
│   %37 = @_5 === nothing
│   %38 = Base.not_int(%37)
└──       goto #8 if not %38
7 ─       goto #2
8 ┄       val = Main.nothing
│         $(Expr(:inbounds, :pop))
│         val
└──       return s
)
2 Likes