# Performance of custom \`Vec\` type versus \`SVector{3, Float64}\`

**URL:** <https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594>\
**Category:** Performance\
**Tags:** staticarrays\
**Created:** [February 1, 2022, 6:48pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594 "2022-02-01T18:48:04Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![CameronBieganek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cameronbieganek/32/6915_2.png) [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)\
**Post date:** [February 1, 2022, 6:48pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594/1 "2022-02-01T18:48:04Z")

</div>

I’m having trouble replicating the performance of `SVector{3, Float64}` with my own custom `Vec` type. Here’s the code:

```julia
using LinearAlgebra
using StaticArrays
using BenchmarkTools

struct Vec
    x::Float64
    y::Float64
    z::Float64
end

Vec(u) = Vec(u[1], u[2], u[3])

function Base.iterate(v::Vec, state=1)
    if state == 1
        v.x, 2
    elseif state == 2
        v.y, 3
    elseif state == 3
        v.z, 4
    else
        nothing
    end
end

Base.length(v::Vec) = 3
Base.eltype(::Type{Vec}) = Float64
Base.IteratorSize(::Type{Vec}) = Base.HasLength()
Base.IteratorEltype(::Type{Vec}) = Base.HasEltype()

vs = [SVector{3}(rand(3)) for _ in 1:10000]
ps = [Vec(svec) for svec in vs]

x = SVector(2.0, 3.0, 4.0)
y = Vec(2, 3, 4)

function foo(vs, x)
    dot.(vs, Ref(x))
end

@btime foo($vs, $x);
@btime foo($ps, $y);

```

And here are the benchmark results:

```julia
julia> @btime foo($vs, $x);
  2.674 μs (2 allocations: 78.17 KiB)

julia> @btime foo($ps, $y);
  16.559 μs (2 allocations: 78.17 KiB)

```

Am I doing something wrong? Is there anything easy I can do to improve the performance of my `Vec` type? One of the guiding principles of the Julia language is that users should be able to get the same performance with custom types that they would get from Base types. However, that doesn’t seem to be the case with StaticArrays (well, I know that StaticArrays is not in Base, but it’s the same idea). They seem to have some extra magic to squeeze out the maximum performance from `SVector`. Or maybe it’s just that the compiler is better at optimizing `NTuple`s?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 1, 2022, 7:03pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594/2 "2022-02-01T19:03:58Z")

</div>

Why not make it an AbstractVector?

```julia

struct Vec <: AbstractVector{Float64}
    x::Float64
    y::Float64
    z::Float64
end

Vec(u) = Vec(u[1], u[2], u[3])

function Base.getindex(v::Vec, i)
    if i == 1
        v.x
    elseif i == 2
        v.y
    elseif i == 3
        v.z
    else
        throw(BoundsError())
    end
end

Base.length(::Vec) = 3
Base.size(::Vec) = (3,)
Base.eltype(::Type{Vec}) = Float64
Base.IteratorSize(::Type{Vec}) = Base.HasLength()
Base.IteratorEltype(::Type{Vec}) = Base.HasEltype()
Base.IndexStyle(::Type{Vec}) = Base.IndexLinear()

function mydot(x,y)
    s = zero(promote_type(eltype(x),eltype(y)))
    @inbounds @simd for i ∈ eachindex(x,y)
        s += x[i]*y[i]
    end
    s
end

LinearAlgebra.dot(x::Vec,y::Vec) = mydot(x,y)

```

Compare just the `@code_typed dot(y,y)` and `@code_typed dot(x,x)`.

---

<div class="post-metadata">

**Author:** ![CameronBieganek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cameronbieganek/32/6915_2.png) [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)\
**Post date:** [February 1, 2022, 7:25pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594/3 "2022-02-01T19:25:58Z")

</div>

Nice, that makes a big difference. The native code appears to be exactly the same now:

```julia
julia> @code_native dot(y, y)
	.text
; ┌ @ REPL[13]:1 within `dot`
	vmovsd	(%rdi), %xmm0 # xmm0 = mem[0],zero
	vmovupd	8(%rdi), %xmm1
; │┌ @ REPL[10]:3 within `mydot`
; ││┌ @ simdloop.jl:77 within `macro expansion` @ REPL[10]:4
; │││┌ @ float.jl:405 within `*`
	vmulpd	8(%rsi), %xmm1, %xmm1
; │││└
; │││┌ @ float.jl:399 within `+`
	vfmadd132sd	(%rsi), %xmm1, %xmm0 # xmm0 = (xmm0 * mem) + xmm1
	vpermilpd	$1, %xmm1, %xmm1 # xmm1 = xmm1[1,0]
	vaddsd	%xmm1, %xmm0, %xmm0
; │└└└
	retq
	nop
; └

julia> @code_native dot(x, x)
	.text
; ┌ @ linalg.jl:205 within `dot`
; │┌ @ linalg.jl:218 within `_vecdot`
; ││┌ @ simdloop.jl:77 within `macro expansion` @ linalg.jl:219
; │││┌ @ generic.jl:905 within `dot`
; ││││┌ @ float.jl:405 within `*`
	vmovsd	(%rdi), %xmm0 # xmm0 = mem[0],zero
	vmovupd	8(%rdi), %xmm1
	vmulpd	8(%rsi), %xmm1, %xmm1
; │││└└
; │││┌ @ float.jl:399 within `+`
	vfmadd132sd	(%rsi), %xmm1, %xmm0 # xmm0 = (xmm0 * mem) + xmm1
	vpermilpd	$1, %xmm1, %xmm1 # xmm1 = xmm1[1,0]
	vaddsd	%xmm1, %xmm0, %xmm0
; │└└└
	retq
	nop
; └

```

And the performance is better, although oddly enough the performance with `Vec` is still slightly slower on my machine, even though the native code appears to be the same:

```julia
julia> @benchmark foo($vs, $x)
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max): 3.014 μs … 55.250 μs ┊ GC (min … max): 0.00% … 72.93%
 Time (median): 5.828 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 6.344 μs ± 3.844 μs ┊ GC (mean ± σ): 6.16% ± 9.12%

      ▃█▁                                                     
  ▂▃▄▅███▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▂ ▂
  3.01 μs Histogram: frequency by time 33.8 μs <

 Memory estimate: 78.17 KiB, allocs estimate: 2.

julia> @benchmark foo($ps, $y)
BenchmarkTools.Trial: 10000 samples with 5 evaluations.
 Range (min … max): 5.768 μs … 112.114 μs ┊ GC (min … max): 0.00% … 82.16%
 Time (median): 6.993 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 7.903 μs ± 6.630 μs ┊ GC (mean ± σ): 6.83% ± 7.68%

   ▃▅▇██▇▅▃ ▁▁▁ ▂
  ██████████▇▅▄▃▁▁▁▄▃▄▄▁▁▁▁▁▃▃▁▁▁▅▄▅▃▃▃▁▁▁▁▁▁▁▁▁▃▄▃▆▇▇█████▇▇ █
  5.77 μs Histogram: log(frequency) by time 19.6 μs <

 Memory estimate: 78.17 KiB, allocs estimate: 2.

```

Version info:

```julia
julia> versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: 11th Gen Intel(R) Core(TM) i9-11900H @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, tigerlake)
Environment:
  JULIA_EDITOR = vim

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 1, 2022, 7:32pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594/4 "2022-02-01T19:32:36Z")

</div>

Odd. FWIW, this fixes it for me:

```julia
julia> @btime foo($ps, $y);
  8.176 μs (2 allocations: 78.17 KiB)

julia> LinearAlgebra.dot(x::Vec,y::Vec) = muladd(x.z,y.z,muladd(x.y,y.y,x.x*y.x))

julia> @btime foo($ps, $y);
  3.958 μs (2 allocations: 78.17 KiB)

julia> @btime foo($vs, $x);
  3.932 μs (2 allocations: 78.17 KiB)

```

Julia+LLVM really do not like vectorizing outer loops, whenever an inner loop is present.  
Thus manually unrolling the inner loop can make a big difference.  
Thus, `foo`, which loops over `vs`, is SIMD if we manually unroll, but not otherwise. =/

---

<div class="post-metadata">

**Author:** ![CameronBieganek](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cameronbieganek/32/6915_2.png) [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)\
**Post date:** [February 1, 2022, 7:47pm UTC](https://discourse.julialang.org/t/performance-of-custom-vec-type-versus-svector-3-float64/75594/5 "2022-02-01T19:47:05Z")

</div>

Thanks! That works for me too:

```julia
julia> @btime foo($vs, $x);
  2.899 μs (2 allocations: 78.17 KiB)

julia> @btime foo($ps, $y);
  2.664 μs (2 allocations: 78.17 KiB)

```

Although I wonder how StaticArrays avoids that SIMD problem, because it looks like their implementation is pretty similar to your first implementation:

[https://github.com/JuliaArrays/StaticArrays.jl/blob/d336cfb672a1d5a427b402c8a3e3026cc75d83e3/src/linalg.jl#L208](https://github.com/JuliaArrays/StaticArrays.jl/blob/d336cfb672a1d5a427b402c8a3e3026cc75d83e3/src/linalg.jl#L208)