Slow in-place broadcasted assignment (.=) of mutable FieldVectors

[stillyslalom]

As noted on Discourse (https://discourse.julialang.org/t/julia-style-arrays-matrix-vs-structs/32424/10?u=stillyslalom), broadcasted assignment seems to be many times slower for mutable FieldVectors than for MVectors.

using StaticArrays

struct Point3D{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end

mutable struct MPoint3D{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end

function f1!(v)
    for i in eachindex(v)
        x, y, z = v[i]
        v[i] = Point3D(2x, x*y, x*z)
    end
end

function f2!(v)
    for i in eachindex(v)
        x, y, z = v[i]
        v[i] .= (2x, x*y, x*z)
    end
end

using BenchmarkTools
v1 = [@SVector(rand(3)) for i = 1:1000]
v2 = [@MVector(rand(3)) for i = 1:1000]
v3 = [Point3D(rand(3))  for i = 1:1000]
v4 = [MPoint3D(rand(3)) for i = 1:1000]

@btime f1!($v1) # 853.985 ns (0 allocations: 0 bytes)
@btime f2!($v2) # 1.775 μs (0 allocations: 0 bytes)

@btime f1!($v3) # 852.917 ns (0 allocations: 0 bytes)
@btime f2!($v4) # 60.091 μs (3000 allocations: 46.88 KiB)

[KristofferC]

How is it for a normal MVector?

[stillyslalom]

That's v2 in the benchmark above:

@btime f2!($v2) # 1.775 μs (0 allocations: 0 bytes)

[KristofferC]

Ah, sorry, missed that.

I think it would be interesting to also add this style of assignment to see how that compares with broadcasting for the different types, https://discourse.julialang.org/t/julia-style-arrays-matrix-vs-structs/32424/10?u=kristoffer.carlsson.

[stillyslalom]

Your field-assignment function is still faster than broadcasting over MVectors, though:

function f3!(v)
   for vv in v
       x, y, z = vv
       vv.x = 2x
       vv.y = x*y
       vv.z = x*z
   end
end

@btime f3!($v4) # 1.158 μs (0 allocations: 0 bytes)

[mateuszbaran]

As far as I can tell indices used in broadcasting over a FieldVector are not constant-propagated to setfield! somehow. We can fix that by passing the index in type of an argument using Val but it's quite ugly. Do you have a better idea?

[c42f]

mutable struct MPoint3D{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end

It's true that this is a performance bug, but I'm also tempted to say "don't do that" :-) Even though we do allow mutable field vectors I personally can't think of a use for them in numerical work. In practice I find that the immutable versions encourage a more functional style in geometric computation which is much more naturally mathematical, far more memory efficient when stored in a large array, and faster.

@mateuszbaran thanks for investigating. Did you have a more minimal piece of code you were looking at in coming to the conclusion about setfield!?

[mateuszbaran]

@mateuszbaran thanks for investigating. Did you have a more minimal piece of code you were looking at in coming to the conclusion about setfield!?

I've inspected the code generated for the problematic case and compared it with this:

julia> using StaticArrays

julia> mutable struct MPoint3D{T} <: FieldVector{3, T}
           x::T
           y::T
           z::T
       end

julia> f(a, i) = a[i] = 100
f (generic function with 1 method)

julia> a = MPoint3D(1,2,3)
3-element MPoint3D{Int64} with indices SOneTo(3):
 1
 2
 3

julia> @code_native f(a, 2)
	.text
; ┌ @ REPL[3]:1 within `f'
	pushq	%r14
	pushq	%rbx
	subq	$56, %rsp
	vxorps	%xmm0, %xmm0, %xmm0
	vmovaps	%xmm0, (%rsp)
	movq	$0, 16(%rsp)
	movq	%fs:0, %rbx
; │┌ @ FieldArray.jl:107 within `setindex!'
	movq	$2, (%rsp)
	movq	-15712(%rbx), %rax
	movq	%rax, 8(%rsp)
	movq	%rdi, %r14
	movq	%rsp, %rax
	movq	%rax, -15712(%rbx)
	movabsq	$jl_box_int64, %rax
	movq	%rsi, %rdi
	callq	*%rax
	movq	%rax, 16(%rsp)
	movq	%r14, 32(%rsp)
	movq	%rax, 40(%rsp)
	movabsq	$139887313344800, %rax  # imm = 0x7F3A0D9ED920
	movq	%rax, 48(%rsp)
	movabsq	$jl_f_setfield, %rax
	leaq	32(%rsp), %rsi
	xorl	%edi, %edi
	movl	$3, %edx
	callq	*%rax
	movq	8(%rsp), %rax
	movq	%rax, -15712(%rbx)
; │└
	movl	$100, %eax
	addq	$56, %rsp
	popq	%rbx
	popq	%r14
	retq
	nopw	%cs:(%rax,%rax)
; └

Although:

julia> g(a) = f(a, 2)
g (generic function with 1 method)

julia> @code_native g(a)
	.text
; ┌ @ REPL[6]:1 within `g'
; │┌ @ REPL[3]:1 within `f'
; ││┌ @ REPL[6]:1 within `setindex!'
	movq	$100, 8(%rdi)
; │└└
	movl	$100, %eax
	retq
	nop
; └

so the constant propagation sometimes works and sometimes it doesn't.

[kylebeggs]

It's true that this is a performance bug, but I'm also tempted to say "don't do that" :-) Even though we do allow mutable field vectors I personally can't think of a use for them in numerical work. In practice I find that the immutable versions encourage a more functional style in geometric computation which is much more naturally mathematical, far more memory efficient when stored in a large array, and faster.

@c42f Solving PDEs with moving boundaries is a use case for me.