在 Julia 中使用比一个慢的线程

Question

我在脚本foo.jl 中有以下代码来执行一个简单的可并行化任务（它可能看起来很普通，但实际上类似于我研究中的一项关键任务）：

using LinearAlgebra
using Statistics

N = 60
length = 200

function matOps(mat1, mat2)
    mat = mat1'mat2 - mat2'mat1
    eye = Matrix{Float64}(I,N,N)
    for i in 0:10000
        mat = eye*transpose(mat)*eye
    end
    return mat
end

y = zeros(length)
println("Number of threads: $(Threads.nthreads())")
@time Threads.@threads for i in 1:length
    mat1 = rand(N,N)
    mat2 = rand(N,N)
    y[i] = mean(matOps(mat1, mat2))
end

所以这一切只是在 60x60 矩阵上执行一些（无意义的）操作 10,000 次，然后取矩阵元素的平均值并将其存储在长度为 200 的数组 y 中。所以我们循环这个 200 长度数组和每个元素我们做这个矩阵运算链，并存储单个平均值。

这是非常可并行化的，因为每次迭代都是彼此完全独立的，但是使用更多线程运行循环实际上会更慢。下面是一些输出：

$:> julia foo.jl

Number of threads: 1
184.268869 seconds (9.85 M allocations: 107.718 GiB, 5.16% gc time)

$:> $env:JULIA_NUM_THREADS=2
$:> julia foo.jl

Number of threads: 2
377.960229 seconds (9.85 M allocations: 107.583 GiB, 71.04% gc time)

$:> $env:JULIA_NUM_THREADS=4
$:> julia foo.jl

Number of threads: 4
1121.259542 seconds (9.84 M allocations: 107.190 GiB, 94.31% gc time)

这实际上给了我使用Distributed 和SharedArrays 而不是Threads 的加速。我也理解Threads 是实验性的，所以我理解这是否只是正在处理的事情。只是想我会检查我是否做错了什么，因为只看时间是如何缩放的，似乎每个线程都在按顺序执行整个循环。

我在 Windows 10 上使用 Powershell，Julia 版本 1.0.3

编辑：我将矩阵运算从eye*mat'eye 更改为eye*transpose(mat)*eye，因为它弄乱了这里的语法。因为垫子是真实的，所以它是一样的。

Answer 1

也许试试这个（注意我们为每个线程有一个单独的 RNG）：

using LinearAlgebra
using Statistics
using Random 

println("Number of threads: $(Threads.nthreads())")

function matOps(mat1, mat2, N)
    mat = mat1'mat2 - mat2'mat1
    eye = Matrix{Float64}(I,N,N)
    for i in 0:10000
        mat = eye*transpose(mat)*eye
    end
    return mat
end
matOps(rand(2,2), rand(2,2), 2) # on rare occassions 
                        # multithreaded compiling fails

function get_y(N, L, rands)
    y = Vector{Float64}(undef, L)    
    Threads.@threads for i in 1:L
        mat1 = rand(@inbounds(rands[Threads.threadid()]),N,N)
        mat2 = rand(@inbounds(rands[Threads.threadid()]),N,N)
        y[i] = mean(matOps(mat1, mat2, N))
    end
    y
end

const N = 60
const L = 20

const rands = [MersenneTwister(i) for i in 1:Threads.nthreads()]
@time get_y(N, L, rands);
@time get_y(N, L, rands);

一个线程：

> julia foo.jl
Number of threads: 1
 15.137731 seconds (958.56 k allocations: 10.771 GiB, 9.40% gc time)
 12.788418 seconds (800.34 k allocations: 10.763 GiB, 7.24% gc time)

两个线程：

> julia foo.jl
Number of threads: 2
  8.573475 seconds (958.52 k allocations: 10.767 GiB, 14.13% gc time)
  7.947657 seconds (800.25 k allocations: 10.756 GiB, 12.73% gc time)