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

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
    

再次感谢各位!

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

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

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

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)
    

感谢各位!

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

给你改了,速度大概提升了一倍吧,我想不出其他提升的方法了。我做的主要改变是使用了 @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

最后,新年快乐!

3赞

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

京ICP备17009874号-2