Matrix-vector product with StaticArrays yields different results in each iteration of an `@batch` loop

[efaulhaber]

MWE:

using StaticArrays: SVector, SMatrix
using Polyester: @batch

function test()
    A = SMatrix{1, 2}(-160.12486538933146, 39.88422637822968)
    x = SVector(-1.000000018394, -1.000000013424)

    result_batch = zeros(2)
    result_serial = zeros(2)

    @batch for i in 1:2
        result = A * x
        result_batch[i] = result[1]
    end

    for i in 1:2
        result = A * x
        result_serial[i] = result[1]
    end

    @info "" result_batch.-result_serial
end

Tested on different architectures with Polyester v0.7.16:

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD Ryzen Threadripper 3990X 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver2)
Threads: 32 default, 0 interactive, 16 GC (on 128 virtual cores)

julia> test()
┌ Info: 
│   result_batch .- result_serial =
│    2-element Vector{Float64}:
│     -1.4210854715202004e-14
└      0.0

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 10 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 6 default, 0 interactive, 3 GC (on 6 virtual cores)

julia> test()
┌ Info: 
│   result_batch .- result_serial =
│    2-element Vector{Float64}:
│     -1.4210854715202004e-14
└      0.0

Happens both with --check-bounds=auto and --check-bounds=yes.

[efaulhaber]

Probably unrelated to Polyester.jl and caused by muladd.

[chriselrod]

Probably unrelated to Polyester.jl and caused by muladd.

Seems like a bug to get different results from one iteration to another.

[efaulhaber]

You're right. That doesn't happen with Threads.@threads.

[efaulhaber]

Here is a more minimal example:

using Polyester: @batch

function test(result)
    x = 1.2

    @batch for i in 1:2
        result[i] = muladd(100.0, -1.0, 100.0 * x)
    end

    return result[1] - result[2]
end

println(test(zeros(2)))

$ julia -t 2 test.jl
-3.552713678800501e-15

This does not happen with Threads.@threads and it doesn't happen when I pass x to the function instead of defining it inside.

[chriselrod]

julia> function test_batch(result)
           x = 1.2
       
           @batch for i in 1:2
               result[i] = muladd(100.0, -1.0, 100.0 * x)
           end
       
           return result[1] - result[2]
       end
test_batch (generic function with 1 method)

julia> function test_serial(result)
           x = 1.2
       
           for i in 1:2
               result[i] = muladd(100.0, -1.0, 100.0 * x)
           end
       
           return result[1] - result[2]
       end
test_serial (generic function with 1 method)

julia> function test_threads(result)
           x = 1.2
       
           Threads.@threads for i in 1:2
               result[i] = muladd(100.0, -1.0, 100.0 * x)
           end
       
           return result[1] - result[2]
       end
test_threads (generic function with 1 method)

julia> t = zeros(2); s = zeros(2); b = zeros(2);

julia> test_serial(s); test_threads(t); test_batch(b);

julia> b .== s
2-element BitVector:
 0
 1

julia> b .== t
2-element BitVector:
 1
 0

When running on N threads, the thread running @batch launches N-1 threads, and then starts the Nth work item itself.
This is why the last item matches serial execution, while the former matches Threads.@threads.

Looks like const-prop results in different values.

Ideally, every iteration would match test_serial/it'd be able to const-prop through the closure.