multithreading - 使用 Julia 1.3.1 的多线程速度较慢

Question

我最近开始使用 1.3.1 版本中的新多线程接口(interface)。在我尝试斐波那契示例之后 in this blog post为了获得显着的加速，我开始尝试我的一些旧算法。

我有一个使用梯形方法来计算曲线下方或上方积分的函数:

function trapezoid( x    :: AbstractVector       ,
                    y    :: AbstractVector       ;
                    y0   :: Number         = 0.0 ,
                    inv  :: Number         = NaN )
   int = zeros(length(x)-1)
   for i = 2:length(x)
      if isnan(inv) == true
         int[i-1] = (y[i]+y[i-1]-2y0) * (x[i]-x[i-1]) / 2
      else
         int[i-1] = (2inv-(y[i]+y[i-1])-2y0) * (x[i]-x[i-1]) / 2
      end # if
   end # for
   integral = sum(int) ;
   return integral
end

然后我有一个非常低效的算法，通过比较曲线下方和上方的面积来确定曲线的中点索引:

function EAM_without_threads( x     :: Vector{Float64} ,
                              y     :: Vector{Float64} ,
                              y0    :: Real            ,
                              ymean :: Real            )

    approx = Vector{Float64}(undef,length(x)-1)
    for i in 1:length(x)-1
        x1 = @view(x[1:i  ])
        x2 = @view(x[i:end])
        y1 = @view(y[1:i  ])
        y2 = @view(y[i:end])
        Al = trapezoid( x1 , y1 , y0=y0     )
        Au = trapezoid( x2 , y2 , inv=ymean )

        approx[i] = abs(Al-Au)
    end
    minind = findmin(approx)[2]

    return x[minind]
end

还有:

function EAM_with_threads( x     :: Vector{Float64} ,
                           y     :: Vector{Float64} ,
                           y0    :: Real            ,
                           ymean :: Real            )

    approx = Vector{Float64}(undef,length(x)-1)
    for i in 1:length(x)-1
        x1 = @view(x[1:i  ])
        x2 = @view(x[i:end])
        y1 = @view(y[1:i  ])
        y2 = @view(y[i:end])
        Al = @spawn trapezoid( x1 , y1 , y0=y0     )
        Au = @spawn trapezoid( x2 , y2 , inv=ymean )

        approx[i] = abs(fetch(Al)-fetch(Au))
    end
    minind = findmin(approx)[2]

    return x[minind]
end

这是我用来尝试这两个功能的:

using SpecialFunctions
using BenchmarkTools

x = collect(-10.0:5e-4:10.0)
y = erf.(x)

然后得到这些结果:

julia> @btime EAM_without_threads(x,y,-1.0,1.0)
  7.515 s (315905 allocations: 11.94 GiB)

julia> @btime EAM_with_threads(x,y,-1.0,1.0)
  10.295 s (1274131 allocations: 12.00 GiB)

我不明白...使用htop我可以看到我的所有8个线程几乎都在满负荷工作。这是我的机器:

julia> versioninfo()
Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-4712MQ CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIA_NUM_THREADS = 8

我知道处理多个线程的开销，并且在小问题中我知道它是否较慢，但为什么在这种情况下？

我也在寻找多线程“良好实践”，因为我猜不是每段代码都会从并行性中受益。

提前谢谢大家。

Answer 1

您的代码在这里做了一些非常多余的工作。它对每个步骤进行完整的梯形积分，而不是仅增量更新 Al 和 Au。在这里，我重写了代码，使其执行零分配，并且我的 EAM 版本在我的计算机上比原始版本快 5 个数量级，且不使用任何线程。

一般来说:在开始研究线程之类的事情之前，请考虑您的算法是否高效。与线程相比，快速算法可以获得更大的加速。

function trapz(x, y; y0=0.0, inv=NaN)
    length(x) != length(y) && error("Input arrays cannot have different lengths")
    s = zero(eltype(x))
    if isnan(inv)
        @inbounds for i in eachindex(x, y)[1:end-1]
            s += (y[i+1] + y[i] - 2y0) * (x[i+1] - x[i])
        end
    else
        @inbounds for i in eachindex(x, y)[1:end-1]
            s += (2inv - (y[i+1] + y[i]) - 2y0) * (x[i+1] - x[i])
        end
    end
    return s / 2
end



function eam(x, y, y0, ymean)
    length(x) != length(y) && error("Input arrays cannot have different lengths")
    Au = trapz(x, y; inv=ymean)
    Al = zero(Au)
    amin = abs(Al - Au)
    ind = firstindex(x)
    @inbounds for i in eachindex(x, y)[2:end-1] # 2:length(x)-1
        Al += (y[i] + y[i-1] - 2y0) * (x[i] - x[i-1]) / 2
        Au -= (2ymean - (y[i] + y[i-1])) * (x[i] - x[i-1]) / 2
        aval = abs(Al - Au)
        if aval < amin
            (amin, ind) = (aval, i)
        end
    end
    return x[ind]
end

此处的基准(我对您的代码使用 @time ，对我自己的代码使用 @btime ，因为使用 @btime< 太耗时了 非常慢的代码):

julia> x = collect(-10.0:5e-4:10.0);

julia> y = erf.(x);

julia> @time EAM_without_threads(x, y, -1.0, 1.0)
 15.611004 seconds (421.72 k allocations: 11.942 GiB, 11.73% gc time)
0.0

julia> @btime eam($x, $y, -1.0, 1.0)
  181.991 μs (0 allocations: 0 bytes)
0.0

一点额外的说明:你不应该写 if isnan(inv) == true，这是多余的。只需编写if isnan(inv)即可。