我有一些代码可以计算一个矩阵中的每个笛卡尔坐标与另一个矩阵中的每个其他坐标之间的距离。对于每个坐标，将返回最小距离以及产生最小值的坐标的索引位置。
function MED3D(m1, m2)
n1::Int = size(m1,1)
Dist = SharedArray{Float64}((n1,3))
@sync @distributed for k in 1:n1
Dist[k,:] = MD3D(m1[k,:], m2, k)
end
return Dist
end
@everywhere function MD3D(v1, m2, k)
dsum::Float64 = Inf
dtemp::Float64 = Inf
i = 0
for j in 1:size(m2,1)
@inbounds dtemp = sqrt((v1[1] - m2[j,1]) * (v1[1] - m2[j,1]) + (v1[2] - m2[j,2]) * (v1[2] - m2[j,2]) + (v1[3] - m2[j,3]) * (v1[3] - m2[j,3]))
if dtemp < dsum
dsum = dtemp
i = j
end
end
return [dsum, k, i]
end
m1 = rand(10,3)
m2 = rand(15,3)
results = MED3D(m1,m2)
虽然这适用于较小的 3D 点云，但我希望通过使用基于 GPU 的分析来提高大型点云的性能。但是，在 Julia 中使用更典型的方法进行矩阵运算似乎是不可能的，因为我必须返回索引位置和最小距离。我已经尝试了几种不同的方法来为这项任务采用 CUarray，但到目前为止，在没有使用实际 for 循环的情况下，它们都失败了。此外，由于将距离矩阵存储在内存中，许多实现它的方法似乎效率极低，对于我的特定数据集，这很快超过了 128gb 的内存。
有人可以帮助我如何在 Julia 中正确实现它以在 GPU 上运行吗？CUarrays 甚至是正确的方法，还是因为除了距离之外我还返回索引，它是否过于抽象？我尝试使用 product 和 dot 计算 L2 范数，但它并没有完全满足我的要求。
更新：
这是我使用广播对内部循环进行 GPUify 的失败尝试。
using CuArrays
function difff(m1,m2)
n1 = size(m1,1)
Dist = Array{Float64}(undef, n1,3)
m2 = CuArray(m2)
m1 = CuArray(m1)
for z in 1:size(m1)
v1 = transpose(m1[z,:])
i = 0
dsum::Float64 = Inf
mi = v1 .- m2
mi = mi .* mi
mi = sum(mi, dims=2)
mi = mi .^ 0.5
mi = findmin(mi)
i = mi[2][1]
dsum = mi[1]
@inbounds Dist[z,:] = [dsum,z,i]
end
end
更新：
尝试 #2 失败。我试图计算最小距离并忘记了索引。这对我的应用程序来说并不理想，但我可以接受它。但是，这只有在第一个数组只有一行时才能正常工作。我试图通过使用 maplices 来解决这个问题，但它不起作用。
using CuArray
a = rand(1,3)
b = rand(3,3)
a = CuArray(a)
b = CuArray(b)
function GK(m1, m2)
reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end
mapslices(GK(b), a, 2)
更新：
通过使用外循环取得进展，但肯定有更好的方法来做到这一点？
using CuArray
using BenchmarkTools
aa = rand(2,3)
bb = rand(5000000,3)
a = CuArray(aa)
b = CuArray(bb)
function GK(m1, m2)
reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end
function D(a,b)
Dist = Array{Float64}(undef,size(a,1),1)
for i in 1:size(a,1)
Dist[i] = GK(a[i,:]',b)
end
return Dist
end
@benchmark test = D(a,b)
@benchmark test = D(aa,bb)
更新：
我之前的分布式版本、修改后的分布式版本、GPU 版本和串行版本之间的一些基准测试。编辑：在扩展到 1000 亿次比较之后，GPU 版本的性能不再优于我之前的分布式版本......关于这是为什么的任何想法？ ??
using Distributed
using SharedArrays
using CuArrays
using BenchmarkTools
aa = rand(4,3)
bb = rand(500000,3)
a = CuArray(aa)
b = CuArray(bb)
function MED3D(m1, m2)
n1::Int = size(m1,1)
Dist = SharedArray{Float64}((n1,1))
@sync @distributed for k in 1:n1
Dist[k] = MD3D(m1[k,:]', m2)
end
return Dist
end
@everywhere function MD3D(v1, m2)
dsum::Float64 = Inf
dtemp::Float64 = Inf
for j in 1:size(m2,1)
@inbounds dtemp = sqrt((v1[1] - m2[j,1]) * (v1[1] - m2[j,1]) + (v1[2] - m2[j,2]) * (v1[2] - m2[j,2]) + (v1[3] - m2[j,3]) * (v1[3] - m2[j,3]))
if dtemp < dsum
dsum = dtemp
end
end
return dsum
end
function MED3DGK(m1, m2)
n1::Int = size(m1,1)
Dist = SharedArray{Float64}((n1,1))
@sync @distributed for k in 1:n1
@inbounds Dist[k] = GK(m1[k,:]',m2)
end
return Dist
end
@everywhere function GK(m1, m2)
reduce(min, sum((m1 .- m2) .^ 2,dims=2) .^ 0.5)
end
function D(a,b)
Dist = Array{Float64}(undef,size(a,1),1)
for i in 1:size(a,1)
@inbounds Dist[i] = GK(a[i,:]',b)
end
return Dist
end
@benchmark test = D(a,b)
@benchmark test = D(aa,bb)
@benchmark test = MED3D(aa,bb)
@benchmark test = MED3DGK(aa,bb)
更新：
使用 NearestNeighbors.jl 和分布式处理实现。关于如何使这更快的任何想法？：
function MED3D(m1, m2)
m2 = Matrix(m2')
kdtree = KDTree(m2)
n1::Int = size(m1,1)
Dist = SharedArray{Float64}((n1,1))
Ind = SharedArray{Float64}((n1,1))
@sync @distributed for k in 1:n1
Ind[k,:], Dist[k,:] = knn(kdtree, m1[k,:], 1)
end
return [Ind,Dist]
end