使用 julia 并行简单数组操作
Simple array operation in parallel using julia
我正在做一个项目，它在一个巨大的数组中包含一些简单的数组操作。
即这里的一个例子
function singleoperation!(A::Array,B::Array,C::Array)
@simd for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
我尝试将其并行化以获得更快的速度。为了并行化它，我使用了 distirbuded 和 share 数组函数，它只是对我刚刚展示的函数做了一点修改：
@everywhere function paralleloperation(A::SharedArray,B::SharedArray,C::SharedArray)
@sync @distributed for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
然而，即使我使用 4 个线程（在 R7-5800x 和 I7-9750H CPU 上试用），两个函数之间也没有时间差异。我能知道我可以在此代码中改进什么吗？非常感谢！我将post下面的完整测试代码：
using Distributed
addprocs(4)
@everywhere begin
using SharedArrays
using BenchmarkTools
end
@everywhere function paralleloperation!(A::SharedArray,B::SharedArray,C::SharedArray)
@sync @distributed for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
function singleoperation!(A::Array,B::Array,C::Array)
@simd for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
N = 128;
A,B,C = fill(0,N,N,N),fill(.2,N,N,N),fill(.3,N,N,N);
AN,BN,CN = SharedArray(fill(0,N,N,N)),SharedArray(fill(.2,N,N,N)),SharedArray(fill(.3,N,N,N));
@benchmark singleoperation!(A,B,C);
BenchmarkTools.Trial: 1612 samples with 1 evaluation.
Range (min … max): 2.582 ms … 9.358 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.796 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.086 ms ± 790.997 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
@benchmark paralleloperation!(AN,BN,CN);
BenchmarkTools.Trial: 1404 samples with 1 evaluation.
Range (min … max): 2.538 ms … 17.651 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.154 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.548 ms ± 1.238 ms ┊ GC (mean ± σ): 0.08% ± 1.65%
如评论所述，这看起来更像是多线程而不是多处理的工作。详细的最佳方法通常取决于您是 CPU-bound 还是 memory-bandwith-bound。通过示例中如此简单的计算，很可能是后者，在这种情况下，您将达到减少 returns 添加额外线程的程度，并且可能想要转向具有显式内存建模的东西， and/or 到 GPU。
但是，一种非常简单的通用方法是使用LoopVectorization.jl
内置的多线程
A = rand(10000,10000)
B = rand(10000,10000)
C = zeros(10000,10000)
# Base
function singleoperation!(A,B,C)
@inbounds @simd for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
using LoopVectorization
function singleoperation_lv!(A,B,C)
@turbo for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
# Multithreaded (make sure you've started Julia with multiple threads)
function threadedoperation_lv!(A,B,C)
@tturbo for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
这给了我们
julia> @benchmark singleoperation!(A,B,C)
BenchmarkTools.Trial: 31 samples with 1 evaluation.
Range (min … max): 163.484 ms … 164.022 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 163.664 ms ┊ GC (median): 0.00%
Time (mean ± σ): 163.701 ms ± 118.397 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▄ ▄▄ ▁ ▁ ▁
▆▁▁▁▁▁▁▁▁▁▁▁▁▆█▆▆█▆██▁█▆▁▆█▁▁▁▆▁▁▁▁▆▁▁▁▁▁▁▁▁█▁▁▁▆▁▁▁▁▁▁▆▁▁▁▁▆ ▁
163 ms Histogram: frequency by time 164 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark singleoperation_lv!(A,B,C)
BenchmarkTools.Trial: 31 samples with 1 evaluation.
Range (min … max): 163.252 ms … 163.754 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 163.408 ms ┊ GC (median): 0.00%
Time (mean ± σ): 163.453 ms ± 130.212 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▃ ▃█▃ █ ▃ ▃
▇▁▁▁▁▁▇▁▁▁▇██▇▁▇███▇▁▁█▇█▁▁▁▁▁▁▁▁█▁▇▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▇▇▁▇▁▇ ▁
163 ms Histogram: frequency by time 164 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark threadedoperation_lv!(A,B,C)
BenchmarkTools.Trial: 57 samples with 1 evaluation.
Range (min … max): 86.976 ms … 88.595 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 87.642 ms ┊ GC (median): 0.00%
Time (mean ± σ): 87.727 ms ± 439.427 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▅ █ ▂ ▂
▅▁▁▁▁██▁▅█▁█▁▁▅▅█▅█▁█▅██▅█▁██▁▅▁▁▅▅▁▅▁▅▁▅▁▅▁▅▁▁▁▅▁█▁▁█▅▁▅▁▅█ ▁
87 ms Histogram: frequency by time 88.5 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
现在，单线程 LoopVectorization @turbo 版本与单线程 @inbounds @simd 版本几乎完美绑定这一事实对我来说暗示我们可能在此处受到内存带宽限制（通常 @turbo 明显快于 @inbounds @simd，所以这个平局表明实际计算不是瓶颈）——在这种情况下，多线程版本只是通过让我们访问更多的内存带宽来帮助我们（虽然随着减少 returns，假设有一些主内存总线无论它可以与多少个内核通信都只能如此之快。
为了获得更深入的了解，让我们尝试使算法更难一些：
function singlemoremath!(A,B,C)
@inbounds @simd for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
using LoopVectorization
function singlemoremath_lv!(A,B,C)
@turbo for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
function threadedmoremath_lv!(A,B,C)
@tturbo for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
果然如此
julia> @benchmark singlemoremath!(A,B,C)
BenchmarkTools.Trial: 2 samples with 1 evaluation.
Range (min … max): 2.651 s … 2.652 s ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.651 s ┊ GC (median): 0.00%
Time (mean ± σ): 2.651 s ± 792.423 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ █
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
2.65 s Histogram: frequency by time 2.65 s <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark singlemoremath_lv!(A,B,C)
BenchmarkTools.Trial: 19 samples with 1 evaluation.
Range (min … max): 268.101 ms … 270.072 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 269.016 ms ┊ GC (median): 0.00%
Time (mean ± σ): 269.058 ms ± 467.744 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁ █ ▁ ▁ ▁▁█ ▁ ██ ▁ ▁ ▁ ▁ ▁
█▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█▁▁█▁▁███▁█▁▁██▁▁▁█▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁▁█ ▁
268 ms Histogram: frequency by time 270 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark threadedmoremath_lv!(A,B,C)
BenchmarkTools.Trial: 56 samples with 1 evaluation.
Range (min … max): 88.247 ms … 93.590 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 89.325 ms ┊ GC (median): 0.00%
Time (mean ± σ): 89.707 ms ± 1.200 ms ┊ GC (mean ± σ): 0.00% ± 0.00%
▄ ▁ ▄█ ▄▄▄ ▁ ▄ ▁ ▁ ▁ ▁ ▁▄
█▁█▆██▆▆▆███▆▁█▁█▁▆▁█▆▁█▆▆▁▁▆█▁▁▁▁▁▁█▆██▁▁▆▆▁▁▆▁▁▁▁▁▁▁▁▁▁▆▆ ▁
88.2 ms Histogram: frequency by time 92.4 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
现在我们更接近 CPU-bound，现在线程和 SIMD 向量化是 2.6 秒和 90 毫秒之间的差异！
如果您的实际问题将像示例问题一样受到内存限制，您可以考虑在 GPU 上工作，在针对内存带宽优化的服务器上，and/or 使用一个包含大量致力于内存建模。
您可能会签出的其他一些软件包可能包括 Octavian.jl (CPU), Tullio.jl (CPU or GPU), and GemmKernels.jl (GPU)。
我正在做一个项目，它在一个巨大的数组中包含一些简单的数组操作。 即这里的一个例子
function singleoperation!(A::Array,B::Array,C::Array)
@simd for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
我尝试将其并行化以获得更快的速度。为了并行化它，我使用了 distirbuded 和 share 数组函数，它只是对我刚刚展示的函数做了一点修改：
@everywhere function paralleloperation(A::SharedArray,B::SharedArray,C::SharedArray)
@sync @distributed for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
然而，即使我使用 4 个线程（在 R7-5800x 和 I7-9750H CPU 上试用），两个函数之间也没有时间差异。我能知道我可以在此代码中改进什么吗？非常感谢！我将post下面的完整测试代码：
using Distributed
addprocs(4)
@everywhere begin
using SharedArrays
using BenchmarkTools
end
@everywhere function paralleloperation!(A::SharedArray,B::SharedArray,C::SharedArray)
@sync @distributed for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
function singleoperation!(A::Array,B::Array,C::Array)
@simd for k in eachindex(A)
@inbounds C[k] = A[k] * B[k] / (A[k] +B[k]);
end
end
N = 128;
A,B,C = fill(0,N,N,N),fill(.2,N,N,N),fill(.3,N,N,N);
AN,BN,CN = SharedArray(fill(0,N,N,N)),SharedArray(fill(.2,N,N,N)),SharedArray(fill(.3,N,N,N));
@benchmark singleoperation!(A,B,C);
BenchmarkTools.Trial: 1612 samples with 1 evaluation.
Range (min … max): 2.582 ms … 9.358 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.796 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.086 ms ± 790.997 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
@benchmark paralleloperation!(AN,BN,CN);
BenchmarkTools.Trial: 1404 samples with 1 evaluation.
Range (min … max): 2.538 ms … 17.651 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.154 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.548 ms ± 1.238 ms ┊ GC (mean ± σ): 0.08% ± 1.65%
如评论所述，这看起来更像是多线程而不是多处理的工作。详细的最佳方法通常取决于您是 CPU-bound 还是 memory-bandwith-bound。通过示例中如此简单的计算，很可能是后者，在这种情况下，您将达到减少 returns 添加额外线程的程度，并且可能想要转向具有显式内存建模的东西， and/or 到 GPU。
但是，一种非常简单的通用方法是使用LoopVectorization.jl
内置的多线程A = rand(10000,10000)
B = rand(10000,10000)
C = zeros(10000,10000)
# Base
function singleoperation!(A,B,C)
@inbounds @simd for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
using LoopVectorization
function singleoperation_lv!(A,B,C)
@turbo for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
# Multithreaded (make sure you've started Julia with multiple threads)
function threadedoperation_lv!(A,B,C)
@tturbo for k in eachindex(A)
C[k] = A[k] * B[k] / (A[k] + B[k])
end
end
这给了我们
julia> @benchmark singleoperation!(A,B,C)
BenchmarkTools.Trial: 31 samples with 1 evaluation.
Range (min … max): 163.484 ms … 164.022 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 163.664 ms ┊ GC (median): 0.00%
Time (mean ± σ): 163.701 ms ± 118.397 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▄ ▄▄ ▁ ▁ ▁
▆▁▁▁▁▁▁▁▁▁▁▁▁▆█▆▆█▆██▁█▆▁▆█▁▁▁▆▁▁▁▁▆▁▁▁▁▁▁▁▁█▁▁▁▆▁▁▁▁▁▁▆▁▁▁▁▆ ▁
163 ms Histogram: frequency by time 164 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark singleoperation_lv!(A,B,C)
BenchmarkTools.Trial: 31 samples with 1 evaluation.
Range (min … max): 163.252 ms … 163.754 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 163.408 ms ┊ GC (median): 0.00%
Time (mean ± σ): 163.453 ms ± 130.212 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▃▃ ▃█▃ █ ▃ ▃
▇▁▁▁▁▁▇▁▁▁▇██▇▁▇███▇▁▁█▇█▁▁▁▁▁▁▁▁█▁▇▁▁▁▁▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▇▇▁▇▁▇ ▁
163 ms Histogram: frequency by time 164 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark threadedoperation_lv!(A,B,C)
BenchmarkTools.Trial: 57 samples with 1 evaluation.
Range (min … max): 86.976 ms … 88.595 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 87.642 ms ┊ GC (median): 0.00%
Time (mean ± σ): 87.727 ms ± 439.427 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▅ █ ▂ ▂
▅▁▁▁▁██▁▅█▁█▁▁▅▅█▅█▁█▅██▅█▁██▁▅▁▁▅▅▁▅▁▅▁▅▁▅▁▅▁▁▁▅▁█▁▁█▅▁▅▁▅█ ▁
87 ms Histogram: frequency by time 88.5 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
现在，单线程 LoopVectorization @turbo 版本与单线程 @inbounds @simd 版本几乎完美绑定这一事实对我来说暗示我们可能在此处受到内存带宽限制（通常 @turbo 明显快于 @inbounds @simd，所以这个平局表明实际计算不是瓶颈）——在这种情况下，多线程版本只是通过让我们访问更多的内存带宽来帮助我们（虽然随着减少 returns，假设有一些主内存总线无论它可以与多少个内核通信都只能如此之快。
为了获得更深入的了解，让我们尝试使算法更难一些：
function singlemoremath!(A,B,C)
@inbounds @simd for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
using LoopVectorization
function singlemoremath_lv!(A,B,C)
@turbo for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
function threadedmoremath_lv!(A,B,C)
@tturbo for k in eachindex(A)
C[k] = cos(log(sqrt(A[k] * B[k] / (A[k] + B[k]))))
end
end
果然如此
julia> @benchmark singlemoremath!(A,B,C)
BenchmarkTools.Trial: 2 samples with 1 evaluation.
Range (min … max): 2.651 s … 2.652 s ┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.651 s ┊ GC (median): 0.00%
Time (mean ± σ): 2.651 s ± 792.423 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
█ █
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
2.65 s Histogram: frequency by time 2.65 s <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark singlemoremath_lv!(A,B,C)
BenchmarkTools.Trial: 19 samples with 1 evaluation.
Range (min … max): 268.101 ms … 270.072 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 269.016 ms ┊ GC (median): 0.00%
Time (mean ± σ): 269.058 ms ± 467.744 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁ █ ▁ ▁ ▁▁█ ▁ ██ ▁ ▁ ▁ ▁ ▁
█▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█▁▁█▁▁███▁█▁▁██▁▁▁█▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁▁█ ▁
268 ms Histogram: frequency by time 270 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark threadedmoremath_lv!(A,B,C)
BenchmarkTools.Trial: 56 samples with 1 evaluation.
Range (min … max): 88.247 ms … 93.590 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 89.325 ms ┊ GC (median): 0.00%
Time (mean ± σ): 89.707 ms ± 1.200 ms ┊ GC (mean ± σ): 0.00% ± 0.00%
▄ ▁ ▄█ ▄▄▄ ▁ ▄ ▁ ▁ ▁ ▁ ▁▄
█▁█▆██▆▆▆███▆▁█▁█▁▆▁█▆▁█▆▆▁▁▆█▁▁▁▁▁▁█▆██▁▁▆▆▁▁▆▁▁▁▁▁▁▁▁▁▁▆▆ ▁
88.2 ms Histogram: frequency by time 92.4 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
现在我们更接近 CPU-bound，现在线程和 SIMD 向量化是 2.6 秒和 90 毫秒之间的差异！
如果您的实际问题将像示例问题一样受到内存限制，您可以考虑在 GPU 上工作，在针对内存带宽优化的服务器上，and/or 使用一个包含大量致力于内存建模。
您可能会签出的其他一些软件包可能包括 Octavian.jl (CPU), Tullio.jl (CPU or GPU), and GemmKernels.jl (GPU)。