Bounds checking in dot

[wraith1995]

The dot product of two static arrays has a bounds check in the inner loop as

using StaticArrays
using LinearAlgebra
a = @SVector rand(4)
b = @SVector rand(4)
@code_typed dot(a,b)

produces

CodeInfo(
1 ──       goto #18 if not true
2 ──       goto #17 if not true
3 ──       nothing::Nothing
4 ┄─ %4  = φ (#3 => 0.0, #16 => %32)::Float64
│    %5  = φ (#3 => 0, #16 => %33)::Int64
│    %6  = (Base.slt_int)(%5, 4)::Bool
└───       goto #17 if not %6
5 ── %8  = (Base.add_int)(%5, 1)::Int64
│    %9  = (Base.sub_int)(%8, 1)::Int64
│    %10 = (Base.add_int)(1, %9)::Int64
└───       goto #14 if not false
6 ── %12 = (Base.slt_int)(0, %8)::Bool
└───       goto #10 if not %12
7 ── %14 = (Base.sle_int)(%10, 4)::Bool
└───       goto #9 if not %14
8 ── %16 = (Base.sle_int)(1, %10)::Bool
└───       goto #11
9 ──       goto #11
10 ─       goto #11
11 ┄ %20 = φ (#8 => %16, #9 => false, #10 => false)::Bool
└───       goto #13 if not %20
12 ─       goto #14
13 ─       invoke Base.throw_boundserror($(QuoteNode(1:4))::UnitRange{Int64}, %8::Int64)::Union{}
└───       $(Expr(:unreachable))::Union{}
14 ┄       goto #15
15 ─       goto #16
16 ─ %27 = (Base.getfield)(a, :data)::NTuple{4,Float64}
│    %28 = (Base.getfield)(%27, %10, false)::Float64
│    %29 = (Base.getfield)(b, :data)::NTuple{4,Float64}
│    %30 = (Base.getfield)(%29, %10, false)::Float64
│    %31 = (Base.mul_float)(%28, %30)::Float64
│    %32 = (Base.add_float)(%4, %31)::Float64
│    %33 = (Base.add_int)(%5, 1)::Int64
│          $(Expr(:simdloop, false))::Any
└───       goto #4
17 ┄ %36 = φ (#4 => %4, #2 => 0.0)::Float64
18 ┄ %37 = φ (#17 => %36, #1 => 0.0)::Float64
└───       goto #19
19 ─       return %37
) => Float64

This seems surely wrong to me, but the dot function uses @inbounds correctly and as far as I can tell, all the needed @propagate_inbounds are in place! What on earth is going wrong?

For reference, I'm on julia v1.1.1 and StaticArrays v0.11.0.

[c42f]

I'm not sure this is an actual problem in practice? The type inferred code looks a bit mysterious here, but the native code looks good:

julia> @code_native debuginfo=:none dot(a, b)
	.text
	vmovupd	(%rdi), %ymm0
	vmulpd	(%rsi), %ymm0, %ymm0
	vextractf128	$1, %ymm0, %xmm1
	vaddpd	%ymm1, %ymm0, %ymm0
	vhaddpd	%ymm0, %ymm0, %ymm0
	vzeroupper
	retq
	nopw	(%rax,%rax)

Note that in #603 @tkoolen changed dot so that loop unrolling is done at the LLVM level which had performance benefits as you can see there.

The origin of the throw_boundserror call in the code_typed output seems to be the indexing of the UnitRange 1:prod(S) in the loop index computation. I do think that's a bit mysterious and it seems to be due to the use of the @simd macro (perhaps interfering somehow with inference). Consider the following MWE with and without the @simd:

function foo(r)
    s = zero(eltype(r))
    @inbounds @simd for i in r
        s += i
    end
    s
end

@code_warntype foo(1:4)

[c42f]

Possibly related:

• make simd loop over ranges perform better JuliaLang/julia#28166
• Change simd loop lowering JuliaLang/julia#29717

[vchuravy]

We don't do DCE elimination on typed IR. The line goto #14 if not false (L11) is the boundscheck elimination