StaticArrays.jl icon indicating copy to clipboard operation
StaticArrays.jl copied to clipboard

Inconsistent matrix multiplication on Julia 1.10

Open danielmatz opened this issue 1 year ago • 1 comments

I get subtly different matrix multiplication results in different contexts in Julia 1.10. This is the simplest example I've come up with so far:

using StaticArrays

struct Foo{A}
    array::A
end

Base.:*(x::Foo, y::Foo) = Foo(x.array * y.array)
Base.:(==)(x::Foo, y::Foo) = x.array == y.array
Base.hash(x::Foo, h::UInt64) = hash((Foo, x.array), h)
Base.isequal(x::Foo, y::Foo) = isequal(x.array, y.array)

function example(x, y)
    a = Foo(x)
    b = Foo(y)
    c = a * b
    d = x * y
    e = Foo(d)
    (; x, y, a, b, c, d, e, failure = c != e)
end

results = map(1:20) do _
    example(@SArray(rand(3, 3)), @SArray(rand(3, 3)))
end

failure_results = filter(results) do result
    result.failure
end

@show length(failure_results) # roughly half the cases fail

f = first(failure_results)
@show f.x * f.y == f.d # this is false!

nothing

I would expect c == e, but roughly half of the cases I run fail that comparison. The example function returns its internal values, and if I try to reproduce the computation of d = x * y, it gives a different result! The difference is only in the last digit of one or two entries in the matrix.

I tried running this example on Julia 1.9, and I didn't see any cases give inconsistent results. So this seems to be new in Julia 1.10.

Here's my version info:

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores
Environment:
  JULIA_ERROR_COLOR = red

danielmatz avatar Feb 05 '24 20:02 danielmatz

The exact result of matrix multiplication depends on the exact code emitted by LLVM. StaticArrays.jl relies on its auto-vectorization features, and the exact code may also depend on what exactly Julia decides to inline in a particular context. So I don't know if we can do anything about this inconsistency without compromising performance.

mateuszbaran avatar Feb 06 '24 10:02 mateuszbaran