StaticArrays.jl
StaticArrays.jl copied to clipboard
Inconsistent matrix multiplication on Julia 1.10
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
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.