在 Julia 中使用线程进行矩阵乘法的奇怪行为

Question

我试图在Threads.@threads for循环中做一些线性代数，它返回一些奇怪的结果。似乎矩阵在循环中没有正确相乘。在线程 for 循环中执行是否安全？

下面是生成矩阵TxR表的最小工作示例NxN。对于每次R迭代，第 (t+1) 矩阵是第 (t) 矩阵与另一个随机矩阵的乘积。乘法由不同的线程执行，然后由一个线程检查正确性。该函数应该返回一个零矩阵。它这样做是为了N<=3，但是，在 的结果中有几个N>=4。

function testMM(T, R, N)
    m1 = zeros(Int64, (T,R,N,N))
    m2 = rand(0:2, (T-1,R,N,N))
    m1[1,:,:,:] = rand(0:1,(R,N,N))
    Threads.@threads for i=1:R
        for t=2:T
            m1[t,i,:,:] = m2[t-1,i,:,:] * m1[t-1,i,:,:]
        end
    end

    odds = zeros(Int64,(T-1,R))
    for i=1:R
        for t=2:T
            if m1[t,i,:,:] != m2[t-1,i,:,:] * m1[t-1,i,:,:]
                odds[t-1,i] = 1
            end
        end
    end
    return odds
end

Threads.nthreads()对我来说是 4。在稳定的 64 位 Julia 0.5.2、0.5.3、0.6.0 上在 Windows 上测试。

编辑：这个例子更简单：矩阵的几个副本独立地平方数次。所有副本的结果应该相同，但该函数通常为 返回 false N>=4。看起来好像不同线程的数据混合在 BLAS 内部的某处。

function testMM2(T, R, N)
m0 = rand(0:2, (N,N))
m = [deepcopy(m0) for i=1:R]
Threads.@threads for i=1:R
    for t=1:T
        m[i] = m[i]^2
    end
end
return all(x->(x==m[1]),m)
end