Odd behavior of threaded matrix multiplication in Julia

Question

I am trying to make some linear algebra in a loop Threads.@threads for, and it returns some weird results. It seems that matrices are not being multiplied properly in a loop. Is it safe to do this in a threading loop?

Below is a minimal working example for creating a table TxRin matrices NxN. For each Riteration, the (t + 1) th matrix is ​​the product of the (t) th matrix with another random matrix. Multiplications are performed by different threads, and then checked for correctness by a single thread. The function should return a matrix of zeros. This is done for N<=3, however, as a result, N>=4there are several.

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() for me 4. Tested on a stable 64-bit version of Julia 0.5.2, 0.5.3, 0.6.0 on 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