mul! appears not to be thread safe

[jeffreyesun]

When I run the following code, the output of foo! is sometimes incorrect, disagreeing with the output of foo_nothread!. The problem only seems to occur when the rmats are larger than 5x1_667_090 or so.

using Random
using Base.Threads: @threads
using LinearAlgebra: mul!
using Statistics: mean

Random.seed!(1234)
lmats = [rand(Float32, (5, 5)) for _=1:2]
rmats = [rand(Float32, (5, 5_000_000)) for _=1:2]
outmats = [zeros(Float32, (5, 5_000_000)) for _=1:2]

function foo!(outmats, lmats, rmats)
    @threads for i in 1:2
        bar!(outmats[i], lmats[i], rmats[i])
    end
end

function foo_nothread!(outmats, lmats, rmats)
    for i in 1:2
        bar!(outmats[i], lmats[i], rmats[i])
    end
end

function bar!(outmat, lmat, rmat)
    mul!(outmat, lmat, rmat)
    @show hash(outmat)
    @show mean(outmat)
end

foo!(outmats, lmats, rmats)
foo_nothread!(outmats, lmats, rmats)

Version info:

Julia Version 1.12.1
Commit ba1e628ee49 (2025-10-17 13:02 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × AMD Ryzen 9 7950X3D 16-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 6 default, 1 interactive, 6 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1
  JULIA_NUM_THREADS = 6

[yuyichao]

How different are the results? I don't imagine blas have much guarantee on reproducibility of the result, especially when there's competitions between threads.

[jeffreyesun]

The results are substantially different, to the point that the code in my project breaks when I introduce multithreading. The calculations are correct some of the time, producing exactly the same hash, but sometimes incorrect, producing not only different hashes but significantly different means.

My guess is that there is some cache reuse problem, perhaps due to the reworking of interactive vs. non-interactive threads in Julia. It looks like a similar issue used to exist but was fixed by implementing thread-local buffers in LinearAlgebra.jl.

I'm currently working around the problem by re-implementing matrix multiplication with nested loops, but that's obviously not ideal.

[mcabbott]

significantly different means

Can you post some actual numbers? Or re-work the example to try N times & pick the worst case?

I see only slight differences, this is fairly typical:

julia> foo!(outmats, lmats, rmats)
hash(outmat) = 0xe00f6ff569150891
hash(outmat) = 0x49a9c059f6efae19
mean(outmat) = 1.2612231f0
mean(outmat) = 1.333006f0

julia> foo_nothread!(outmats, lmats, rmats)
hash(outmat) = 0x49a9c059f6efae19
mean(outmat) = 1.2612238f0
hash(outmat) = 0x5baf39933702562e
mean(outmat) = 1.3330061f0

julia> 1.2612238 / 1.2612231
1.0000005550167927

julia> versioninfo()
Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 11 × Apple M3 Pro

[jeffreyesun]

Sure thing, here is an example:

julia> foo!(outmats, lmats, rmats)
hash(outmat) = 0x49a9c059f6efae19
hash(outmat) = 0xd5bf3168a3281f74
mean(outmat) = 1.3330413f0
mean(outmat) = 1.2612238f0

julia> foo_nothread!(outmats, lmats, rmats)
hash(outmat) = 0x49a9c059f6efae19
mean(outmat) = 1.2612238f0
hash(outmat) = 0x5baf39933702562e
mean(outmat) = 1.3330061f0

If I examine the resulting matrices, I get

julia> maximum(abs.(mat_1 ./ mat_2 .- 1))
2.6345253f-5

It's ends up being unsuitable for my application for the results to differ by as much as this. What's more, the fact that the threaded version often produces exactly the same output and hash as the non-threaded version leads me to believe that some incorrect state is being reached sometimes, stochastically. Maybe the result of race conditions or cache reuse? Unfortunately, I don't understand the internals here very well.

[yuyichao]

This amount of difference is more or less expected and not a bug. It’s also nothing specific to threading or any threading related bugs. This is simply floating point math being not reassociative so different computation order gives you different result.