标量优先似乎要求挺高？

XuJingye · 2023 年3 月 21 日 09:19

# 准备初始相空间坐标
zvec = randn(10000)*0.12
δvec = randn(10000)*1.06e-3

const restenergy = 0.511e6            # [eV]

nturns = 10001   # 圈数

struct Par
    v1::Float64
    r::Float64
    ϕs::Float64
    ϕ2s::Float64
    h1::Int64
    h::Int64
    circum::Float64
    centerenergy::Float64
    αc::Float64
end
Base.broadcastable(x::Par) = Ref(x)

function revolution(z, δ, p::Par)::Tuple{Float64,Float64}
    δ = δ + p.v1 / p.centerenergy * (sin(p.ϕs - 2 * p.h1 * π / p.circum * z ) - sin(p.ϕs) + p.r * sin(p.ϕ2s - 2 * p.h1 * p.h * π / p.circum * z) - p.r * sin(p.ϕ2s))
    η = p.αc - 1 / (p.centerenergy * (1 + δ) / restenergy)^2
    z = z - η * δ * p.circum
    z, δ
end

function revolution(z, δ, k1, k2, c1, c2, c3, p)::Tuple{Float64,Float64}
    δ = δ + c1 * (sin(p.ϕs - k1 * z ) + p.r * sin(p.ϕ2s - k2 * z) - c2)
    η = p.αc - c3 / (1+δ)^2
    z = z - η * δ * p.circum
    z, δ
end

function main!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
    k1 = 2 * p.h1 * π / p.circum
    k2 = 2 * p.h1 * p.h * π / p.circum
    γ = p.centerenergy/restenergy
    c1 = p.v1 / p.centerenergy
    c2 = sin(p.ϕs) + p.r * sin(p.ϕ2s)
    c3 = 1/γ^2
    @inbounds for i = 1:size(zvec, 1)
        z, δ = zvec[i], δvec[i]
        @inbounds for j = 1:nturns
            z, δ = revolution(z, δ, k1, k2, c1, c2, c3, p)   # 这个减少了参数的重复计算
            # z, δ = revolution(z, δ, p)    # 这个参数重复计算很费时间
        end
        zvec[i], δvec[i] = z, δ
    end
end

p = Par(3.6395053705870048e+06, 0.1801760487622639, 2.1575815005097385, 6.0620746573056303, 756, 3, 1360.4, 6e9, 0.0000156)
@btime main!(zvec, δvec, nturns, p)

我按照物理公式，定义了函数：

revolution(z, δ, p)

运行时间：

  8.505 s (0 allocations: 0 bytes)

  8.486 s (0 allocations: 0 bytes)

  8.603 s (0 allocations: 0 bytes)

离MATLAB的最快3.7s差了很多，甚至比矢量化的3.973/3.991Julia写法慢很多。

我猜测主要是“标量优先”写法中，重复的参数计算太多。于是定义了

revolution(z, δ, k1, k2, c1, c2, c3, p)

这个函数，主要是把能先算的都算好了。然而时间为

  4.777 s (0 allocations: 0 bytes)

  4.880 s (0 allocations: 0 bytes)

  4.889 s (0 allocations: 0 bytes)

  4.887 s (0 allocations: 0 bytes)

可以看到，重复的参数计算确实影响很大，但是离矢量化写法还是有不小的差距。

感觉，在写标量函数的时候，“参数的重复计算”和“内存的调用”之间存在一种取舍。部分情况下，你“标量优先”做出的无用重复计算开销，会大于“加载内存”开销。

有没有开发者大佬能把Julia这个问题给解决一下的……

另一个问题，我用@code_warntype看到

julia> @code_warntype main!(zvec, δvec, nturns, p)
Variables
  #self#::Core.Const(main!)
  zvec::Vector{Float64}
  δvec::Vector{Float64}
  nturns::Int64
  p::Par
  @_6::Union{Nothing, Tuple{Int64, Int64}}    # 看这儿
  ηvec::Vector{Float64}
  i::Int64

Body::Nothing
# 以下忽略

上面出现了类似警告的颜色，其中有个Nothing是为什么呢？如何消除呢？

gxliu · 2023 年3 月 22 日 14:02

关于重复计算的问题，我觉得从设计代码的人这边更好解决，而不是指望编译器解决所有（毕竟 Julia 编译器不是550 W）。还有就是，如果不重新写 revolution() 的话，代码改成向量化还是有很多重复计算。

另外，从未贴出的另一部分 @code_llvm 输出可以看到 @_6 就是对循环 for i = 1:size(zvec, 1) 生成的迭代器。这个是黄色提示，基本不用管。

│   %29 = Main.size(zvec, 1)::Int64
│   %30 = (1:%29)::Core.PartialStruct(UnitRange{Int64}, Any[Core.Const(1), Int64])    
│         (@_6 = Base.iterate(%30))
│   %32 = (@_6 === nothing)::Bool
│   %33 = Base.not_int(%32)::Bool
└──       goto #7 if not %33

XuJingye · 2023 年3 月 23 日 02:12

XuJingye:

function revolution(z, δ, p::Par)::Tuple{Float64,Float64}
    δ = δ + p.v1 / p.centerenergy * (sin(p.ϕs - 2 * p.h1 * π / p.circum * z ) - sin(p.ϕs) + p.r * sin(p.ϕ2s - 2 * p.h1 * p.h * π / p.circum * z) - p.r * sin(p.ϕ2s))
    η = p.αc - 1 / (p.centerenergy * (1 + δ) / restenergy)^2
    z = z - η * δ * p.circum
    z, δ
end

function revolution(z, δ, k1, k2, c1, c2, c3, p)::Tuple{Float64,Float64}
    δ = δ + c1 * (sin(p.ϕs - k1 * z ) + p.r * sin(p.ϕ2s - k2 * z) - c2)
    η = p.αc - c3 / (1+δ)^2
    z = z - η * δ * p.circum
    z, δ
end

主要是这两个比起来，上面那个，我们这个专业的人一看就知道物理量是什么。下面的那个，我自己下次调用的话，都要瞪半天。

然后又有两个新的问题。

问题一、大佬能更改我目前for循环写法第二个revolution，以达到和下面main1!()矢量化一样的速度，或者示范更改下面矢量化写法的main2!()吗？这两个在我看来，计算次数已经压榨到极限了，如果可以进一步提升速度的话，我学习学习是怎么处理的（不使用多线程）。

附上我矢量化的代码：

function main1!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
    ηvec = zeros(size(zvec, 1))
    for i = 1:nturns
        δvec .= δvec .+ p.v1 / p.centerenergy .* (sin.(p.ϕs .- 2 * p.h1 * π / p.circum .* zvec) .- sin(p.ϕs) .+ p.r .* sin.(p.ϕ2s .- 2 * p.h1 * p.h * π / p.circum .* zvec) .- p.r * sin(p.ϕ2s))
        ηvec .= p.αc .- 1 ./ (p.centerenergy/restenergy .* (1 .+ δvec)) .^ 2
        zvec .= zvec .- p.circum .* ηvec .* δvec
    end
end

function main2!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
    k1 = 2 * p.h1 * π / p.circum
    k2 = 2 * p.h1 * p.h * π / p.circum
    c1 = p.v1 / p.centerenergy
    c2 = sin(p.ϕs) + p.r * sin(p.ϕ2s)
    c3 = 1 / (p.centerenergy / restenergy)^2

    ηvec = zeros(size(zvec, 1))
    for i = 1:nturns
        δvec .= δvec .+ c1 .* (sin.(p.ϕs .- k1 .* zvec) .+ p.r .* sin.(p.ϕ2s .- k2 .* zvec) .- c2)
        ηvec .= p.αc .- c3 ./ (1 .+ δvec) .^ 2
        zvec .= zvec .- p.circum .* ηvec .* δvec
    end
end

问题二、我上面main2!()代码，是我尝试去减少循环中的计算，但是速度居然还慢一丢丢。我不理解。

# 第一个
  1.958 s (2 allocations: 78.17 KiB)
  1.956 s (2 allocations: 78.17 KiB)
  1.959 s (2 allocations: 78.17 KiB)
  1.958 s (2 allocations: 78.17 KiB)
  1.961 s (2 allocations: 78.17 KiB)
# 第二个看起来更快的实际更慢
  1.989 s (2 allocations: 78.17 KiB)
  1.983 s (2 allocations: 78.17 KiB)
  1.982 s (2 allocations: 78.17 KiB)
  1.993 s (2 allocations: 78.17 KiB)
  1.988 s (2 allocations: 78.17 KiB)

johnnychen94 · 2023 年3 月 23 日 14:26

重复的参数计算确实影响很大，但是离矢量化写法还是有不小的差距 。

也许你可以考虑交换一下循环顺序。

@inbounds for j = 1:nturns
    @inbounds for i = 1:size(zvec, 1)
        zvec[i], δvec[i] = revolution(zvec[i], δvec[i], k1, k2, c1, c2, c3, p)   # 这个减少了参数的重复计算
    end
end

简单尝试了一下，快了大约一倍。

离MATLAB的最快3.7s 差了很多

我自己的经验，Julia 的 sin/cos/exp 这些函数的向量化比 MATLAB/numpy 的版本要慢一些。没有去深究为什么，我的猜测是它们调用了更好的向量化实现，而 Julia 的向量化实现似乎比较朴素。

@_6::Union{Nothing, Tuple{Int64, Int64}}

Julia 的 for 循环采用的迭代器协议是 iterate(x)::Union{Nothing, T} 的形式，其中返回 Nothing 表示迭代结束（类似于 Python 的 StopIteration）。对于这种类型可能性有限，并且可能的类型都是具体类型的情况，Julia 有所谓的 union-splitting 技术来优化它，所以你看到的警告颜色是黄色（也就是不太需要关心）

johnnychen94 · 2023 年3 月 23 日 14:34

高性能的实现本身可能会和纸面的推导有一些差异导致可读性变差，比如说 ImageMorphology.extreme_filter。我个人的倾向是引入辅助函数和注释说明：

function _revolution_cache(p)
    k1 = 2 * p.h1 * π / p.circum
    k2 = 2 * p.h1 * p.h * π / p.circum
    γ = p.centerenergy/restenergy
    c1 = p.v1 / p.centerenergy
    c2 = sin(p.ϕs) + p.r * sin(p.ϕ2s)
    c3 = 1/γ^2
    return (; k1, k2, c1, c2, c3)
end

function revolution(z, δ, p::Par, c)::Tuple{Float64,Float64}
    # c is the precomputed intermediate result that can be reused

    # δ = δ + p.v1 / p.centerenergy * (sin(p.ϕs - 2 * p.h1 * π / p.circum * z ) - sin(p.ϕs) + p.r * sin(p.ϕ2s - 2 * p.h1 * p.h * π / p.circum * z) - p.r * sin(p.ϕ2s))
    δ = δ + c.c1 * (sin(p.ϕs - c.k1 * z ) + p.r * sin(p.ϕ2s - c.k2 * z) - c.c2)
    # η = p.αc - 1 / (p.centerenergy * (1 + δ) / restenergy)^2
    η = p.αc - c.c3 / (1+δ)^2
    z = z - η * δ * p.circum
    z, δ
end

N5N3 · 2023 年3 月 25 日 09:57

Base.sin/cos/exp 默认都是不内联，因此无法向量化。
这种时候个人建议是能上LV就上LV。如果LV不支持，只能调SLEEF.jl了（记得内联）

XuJingye · 2023 年3 月 28 日 00:03

你好，前阵子一直没有机会测试。刚刚试了下，按您所说，确实能提速：

function main1!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
    ...略...
    @inbounds for i = 1:size(zvec, 1)
        z, δ = zvec[i], δvec[i]
        @inbounds for j = 1:nturns
            z, δ = revolution(z, δ, k1, k2, c1, c2, c3, p)
        end
        zvec[i], δvec[i] = z, δ
    end
end

function main2!(zvec::Vector{Float64}, δvec::Vector{Float64}, nturns, p::Par)
    ...略...
    @inbounds for j = 1:nturns
        @inbounds for i = 1:size(zvec, 1)
            zvec[i], δvec[i] = revolution(zvec[i], δvec[i], k1, k2, c1, c2, c3, p)
        end
    end
end

julia> @time main1!(zvec, δvec, nturns, p)
  3.056780 seconds (371.13 k allocations: 19.812 MiB, 0.60% gc time, 3.34% compilation time)
  2.963028 seconds
  2.976086 seconds
  2.973847 seconds
  2.975546 seconds

julia> @time main2!(zvec, δvec, nturns, p)
  2.165210 seconds (44.41 k allocations: 2.447 MiB, 0.63% compilation time)
  2.160676 seconds
  2.139684 seconds
  2.137427 seconds

大概从3.0 s到2.15 s。

========================================

但我完全不能理解：我特意为了避免反复读写数组，特意先从数组中拿一个粒子坐标出来，跑到头再算下一个粒子。你这里是每一圈逐粒子计算，每圈都要读写，怎么会比我的快的？

julia> sizeof(zvec)
80000

可以看到zvec都有80kB，加上δvec有160kB，这俩肯定只能放在下图L1 Cache中；而如果只拿一个z,δ，就直接存在寄存器中，速度最快。我这么想哪儿错了吗？

XuJingye · 2023 年3 月 28 日 00:27

这个写法是看的哪里的？在官方文档里随便搜了下，没搜到

XuJingye · 2023 年3 月 28 日 00:40

请问一下内联是什么东西？我直接将SLEEF.sin去做矢量化运算，速度比原版慢得多，估计问题就是这个内联没弄

johnnychen94 · 2023 年3 月 28 日 04:22

你这里是每一圈逐粒子计算，每圈都要读写，怎么会比我的快的？

高性能计算中存在缓存局部性的概念 (Locality): CPU 在读取数据进入缓存的时候是批量读取的 (cacheline)，因此如果你在较短的时间内对矩阵进行迭代操作的话，数据依然在缓存中的可能性比较高。所以你特意调整循环的顺序可能并没有按照预期产生性能上的优化。

我的理解是这类优化的有效性大多数需要经验和尝试来做判断，应该不太有一个绝对不变的答案（比如 for 循环内部计算变复杂/简单之后，情况可能也会跟着变化）。

有兴趣的话还可以阅读 https://biojulia.github.io/post/hardware/

(; k1, k2, c1, c2, c3)

这个写法是看的哪里的？在官方文档里随便搜了下，没搜到

额，这个是创建 NamedTuple 的语法糖；我也不记得在哪看到的了。

N5N3 · 2023 年3 月 28 日 07:33

内联 = @inline

julia> a = randn(1024); b = sin.(a);

julia> import Tools.MyFastFloat: SLEEF

julia> a = randn(1024); b = sin.(a);

julia> @btime $b .= sin.($a);
  9.000 μs (0 allocations: 0 bytes)

julia> @btime $b .= SLEEF.sin.($a);
  29.900 μs (0 allocations: 0 bytes)

julia> @inline xsin(x) = @inline SLEEF.sin(x);

julia> @btime $b .= xsin.($a);
  7.300 μs (0 allocations: 0 bytes)

julia> @inline xsin_fast(x) = @inline SLEEF.sin_fast(x);

julia> @btime $b .= xsin_fast.($a);
  2.311 μs (0 allocations: 0 bytes)

以上计算我们实际上是依靠LLVM的Vectorlizer进行向量化, 因此性能相比LV还是有比较大的差距的。
此外，对于仅使用实数的场景，如果LV不能正常运行，那LLVM大概率也无法向量化。

gxliu · 2023 年3 月 30 日 03:37

这篇博客我翻译完了，也可阅读中文版。

photor · 2023 年5 月 2 日 13:23

我用的Julia 1.7.3版提示错误，似乎不支持这样的@inline语法

N5N3 · 2023 年5 月 2 日 13:45

这个是1.8的新特性，先前的版本应该没有很方便的办法在调用端强制inline