如何提高双重for循环的效率

综合讨论区 性能优化
1

for 循环提速问题

本人是Julia新手，下面代码中position和pos是两个非常大的数组，size分别是(300000, 3)和(9000, 3)，我需要对两个数组中的元素进行一次两两对比，来进行距离筛选。希望大家能帮助一下，提高计算效率。现在的版本计算一次我估计得20多分钟，我就没有跑完就关闭啦。希望能够多线程，并行啥的？计算机有很多核，希望可以充分使用起来，感谢大家！

Julia版本 1.5

运行环境 Jupyter

  1. 代码如下
    '''
此处position和pos都是空间三维点坐标，测试性能用随机数组应该也是一样的。比如假设：
position = rand(300000, 3).*500
pos = rand(9000, 3).*500
xlo, xhi = minimum(position[:, 1]), maximum(position[:, 1])
ylo, yhi = minimum(position[:, 2]), maximum(position[:, 2])
zlo, zhi = minimum(position[:, 3]), maximum(position[:, 3])
'''
function test(position, pos, xlo, xhi, ylo, yhi, zlo, zhi)
    pingheng = 3.2
    cellx, celly, cellz = xhi-xlo, yhi-ylo, zhi-zlo #都是一些浮点数
    marks = [] #用作标记筛选
    for i in 1:size(position, 1)
        for j in 1:size(pos, 1)
            L = position[i, :] - pos[j, :]

            Lx = L[1] - floor(L[1]/cellx+0.5)*cellx #考虑周期性边界
            Ly = L[2] - floor(L[2]/celly+0.5)*celly
            Lz = L[3] - floor(L[3]/cellz+0.5)*cellz

            r = (Lx^2 + Ly^2 + Lz^2)^0.5

            if r < pingheng
                push!(marks, false) #用于筛选
                break
            end
        end
    end
end

再次感谢各位！

2

给一个测试用例啊。能给一个手算的案例就更好了，帮助别人理解你的代码

3

你好，我重新编辑了一下问题，现在这个函数可以去测试计算啦，谢谢！

4

我不是很懂你这个的原理，另外没有输出，我也不敢乱改（不知道如何判断结果是否正确，所以你一定要有return 语句），所以你至少要提供一段可以运行（并且要输出正确结果，提供运行时间以及内存消耗给我们做参考）的代码。还有，你这个测试的数据太大了，建议改小一点（一般10秒钟以内运行完吧,你这个我消耗了一个小时的cpu 时间都没运行完。）

5

for 循环提速问题

本人是Julia新手，下面代码中position和pos是两个非常大的数组，size分别是(300000, 3)和(9000, 3)，我需要对两个数组中的元素进行一次两两对比，来进行距离筛选。希望大家能帮助一下，提高计算效率，现在我估计一下，我估计得20多分钟。多线程，并行啥的？计算机有很多核，希望可以以充分使用起来，感谢大家！

Julia版本 1.5

运行环境 Jupyter

  1. 代码如下
    """
此处position和pos都是空间三维点坐标，测试性能用随机数组应该也是一样的。比如假设：
position = rand(1000, 3).*100
pos = rand(1000, 3).*100
"""
function test(position::Array{Float64,2}, pos::Array{Float64,2})::Array{Float64,2}
    
    xlo, xhi = minimum(position[:, 1]), maximum(position[:, 1])
    ylo, yhi = minimum(position[:, 2]), maximum(position[:, 2])
    zlo, zhi = minimum(position[:, 3]), maximum(position[:, 3])
    rp = 3.0825 ::Float64 #设置的平衡距离
    cellx, celly, cellz = xhi-xlo, yhi-ylo, zhi-zlo ::Float64
    marks = [] 
    for i in 1:size(position, 1)
        for j in 1:size(pos, 1)
            L = position[i, :] - pos[j, :]

            Lx = L[1] - floor(L[1]/cellx+0.5)*cellx #考虑周期性边界
            Ly = L[2] - floor(L[2]/celly+0.5)*celly
            Lz = L[3] - floor(L[3]/cellz+0.5)*cellz

            r = (Lx^2 + Ly^2 + Lz^2)^0.5

            if r < rp
                push!(marks, i) #用作标记，那些坐标点距离小于平横距离 rp
                break
            end
        end
    end
    return position[marks, :]
end
  2. 程序测试如下
    @time test(position::Array{Float64,2}, pos::Array{Float64,2})::Array{Float64,2};
0.394665 seconds (2.83 M allocations: 302.081 MiB, 18.21% gc time)

感谢各位！

6

不好意思，昨天有点事情，现在算例变小啦，可以直接运行。我在考虑除了提高for循环的效率，就是我的写法应该有问题，没有发挥出julia的性能，我主要是想把这个结构，如何写成多线程计算这种，当然论坛里面之前的帖子讨论Cell-List的算法，我也在学习着写。您要是能帮助并行这一块就最好了，感谢。帖子我编辑不了，所以只能直接回复您啦！

7

给你改了，速度大概提升了一倍吧，我想不出其他提升的方法了。我做的主要改变是使用了 @views 宏避免了一部分内存分配，又使用了一个 BitVector 来避免push 操作。至于多线程和分布式计算，我只会用 Thread.@threads宏，具体用法看文档，这里我觉得提升不大。另外建议读一下文档的性能建议（觉得自己写的代码慢的时候，我也常常会去读）。

function test3(position=position, pos=pos)
    N=size(position,1)
    
    @views xlo, xhi = minimum(position[:, 1]), maximum(position[:, 1])
    @views ylo, yhi = minimum(position[:, 2]), maximum(position[:, 2])
    @views zlo, zhi = minimum(position[:, 3]), maximum(position[:, 3])
    rp = 3.0825 ::Float64 #设置的平衡距离
#     cellx, celly, cellz = xhi-xlo, yhi-ylo, zhi-zlo ::Float64
    cells=Float64[xhi-xlo, yhi-ylo, zhi-zlo] #cells=[ cellx, celly, cellz]
    Ls=rand(3) #Ls=[Lx,Ly,Lz]
    marks = falses(N)
    for i in 1:N
        for j in 1:N
            @views L = position[i, :] - pos[j, :]

            for k in 1:3
                Ls[k]=L[k] -floor(L[k]/cells[k]+0.5)*cells[k]
            end
#             Lx = L[1] - floor(L[1]/cellx+0.5)*cellx #考虑周期性边界
#             Ly = L[2] - floor(L[2]/celly+0.5)*celly
#             Lz = L[3] - floor(L[3]/cellz+0.5)*cellz

#             r2 = (Lx^2 + Ly^2 + Lz^2)
            r2=sum(a^2 for a in Ls)

            if r2 < rp^2
                marks[i]=true
                break
            end
        end
    end
    return position[marks, :]
end

加入多线程

function test4(position=position, pos=pos)
    N=size(position,1)
    
    @views xlo, xhi = minimum(position[:, 1]), maximum(position[:, 1])
    @views ylo, yhi = minimum(position[:, 2]), maximum(position[:, 2])
    @views zlo, zhi = minimum(position[:, 3]), maximum(position[:, 3])
    rp = 3.0825 ::Float64 #设置的平衡距离
    cells=Float64[xhi-xlo, yhi-ylo, zhi-zlo] #cells=[ cellx, celly, cellz]
- Ls=rand(3) 
    marks = falses(N)
    Threads.@threads for i in 1:N
+        Ls=rand(3) #移到这里，保证每个线程有一个自己的数组用于写入。
        for j in 1:N
            @views L = position[i, :] - pos[j, :]
            for k in 1:3
                Ls[k]=L[k] -floor(L[k]/cells[k]+0.5)*cells[k]
            end
            r2=sum(a^2 for a in Ls)

            if r2 < rp^2
                marks[i]=true
                break
            end
        end
    end
    return position[marks, :]
end

放一下benchmark.

using BenchmarkTools

你原始代码的benchmark

julia>@benchmark test(position::Array{Float64,2}, pos::Array{Float64,2})::Array{Float64,2}
BenchmarkTools.Trial: 
  memory estimate:  303.19 MiB
  allocs estimate:  2838215
  --------------
  minimum time:     89.933 ms (5.71% GC)
  median time:      91.152 ms (6.75% GC)
  mean time:        91.693 ms (6.70% GC)
  maximum time:     100.230 ms (6.65% GC)
  --------------
  samples:          55
  evals/sample:     1

提速后的 benchmark.

julia> @benchmark test3()
BenchmarkTools.Trial: 
  memory estimate:  101.05 MiB
  allocs estimate:  946017
  --------------
  minimum time:     41.323 ms (5.88% GC)
  median time:      42.436 ms (5.74% GC)
  mean time:        42.685 ms (6.50% GC)
  maximum time:     50.683 ms (4.46% GC)
  --------------
  samples:          118
  evals/sample:     1

julia> @benchmark test4()
BenchmarkTools.Trial: 
  memory estimate:  101.24 MiB
  allocs estimate:  947789
  --------------
  minimum time:     9.178 ms (0.00% GC)
  median time:      10.537 ms (0.00% GC)
  mean time:        35.980 ms (68.08% GC)
  maximum time:     228.600 ms (93.30% GC)
  --------------
  samples:          139
  evals/sample:     1

最后，新年快乐！

4 个赞
8

哇，非常感谢您的如此详细的回答，我下去再测试一下，有问题看文档不懂再来请教您。新年快乐！！！ :yum: