使用 Julia @parallel 和 SharedArray 进行并行计算

Question

我一直在尝试使用 @parallel 和 SharedArrays 在 Julia 中实现一些并行编程。

Xi = Array{Float64}([0.0, 450.0, 450.0, 0.0, 0.0, 450.0, 450.0, 0.0])
Yi = Array{Float64}([0.0, 0.0, 600.0, 600.0, 0.0, 0.0, 600.0, 600.0])
Zi = Array{Float64}([0.0, 0.0, 0.0, 0.0, 400.0, 400.0, 400.0, 400.0])
Xj = Array{Float64}([0.0, 450.0, 450.0, 0.0, 0.0, 450.0, 450.0, 0.0])
Yj = Array{Float64}([0.0, 0.0, 600.0, 600.0, 0.0, 0.0, 600.0, 600.0])
Zj = Array{Float64}([0.0, 0.0, 0.0, 0.0, 400.0, 400.0, 400.0, 400.0])
L = Array{Float64}([400.0, 400.0, 400.0, 400.0, 450.0, 600.0, 450.0, 600.0])
Rot = Array{Float64}([90.0, 90.0, 90.0, 90.0, 0.0, 0.0, 0.0, 0.0])

显然，这些向量会很大，但为简单起见，我只放了这个有限的大小。

这是没有并行计算的操作：

function jt_transcoord(Xi, Yi, Zi, Xj, Yj, Zj, Rot, L)
r = Vector(length(Xi))
for i in 1:length(Xi)
    rxX = (Xj[i] - Xi[i]) / L[i]
    rxY = (Yj[i] - Yi[i]) / L[i]
    rxZ = (Zj[i] - Zi[i]) / L[i]
        if rxX == 0 && rxY == 0
            r[i] = [0 0 rxZ; cosd(Rot[i]) -rxZ*sind(Rot[i]) 0; sind(Rot[i]) rxZ*cosd(Rot[i]) 0]
        else
            R=sqrt(rxX^2+rxY^2)
            r21=(-rxX*rxZ*cosd(Rot[i])+rxY*sind(Rot[i]))/R
            r22=(-rxY*rxZ*cosd(Rot[i])-rxX*sind(Rot[i]))/R
            r23=R*cosd(Rot[i])
            r31=(rxX*rxZ*sind(Rot[i])+rxY*cosd(Rot[i]))/R
            r32=(rxY*rxZ*sind(Rot[i])-rxX*cosd(Rot[i]))/R
            r33=-R*sind(Rot[i])
            r[i] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]
        end
end
return r
end

返回值基本上是一个数组，每个向量行中都包含一个矩阵。看起来像这样：

r = 
[[0.0 0.0 1.0; 0.0 -1.0 0.0; 1.0 0.0 0.0], 
[0.0 0.0 1.0; 0.0 -1.0 0.0; 1.0 0.0 0.0], 
[0.0 0.0 1.0; 0.0 -1.0 0.0; 1.0 0.0 0.0], 
[0.0 0.0 1.0; 0.0 -1.0 0.0; 1.0 0.0 0.0], 
[1.0 0.0 0.0; 0.0 -0.0 1.0; 0.0 -1.0 -0.0], 
[0.0 1.0 0.0; 0.0 -0.0 1.0; 1.0 0.0 -0.0], 
[-1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 -0.0], 
[0.0 -1.0 0.0; -0.0 0.0 1.0; -1.0 -0.0 -0.0]]

这是我使用@parallel 的函数。首先我需要将向量转换为SharedArrays：

Xi = convert(SharedArray, Xi)
Yi = convert(SharedArray, Yi)
Zi = convert(SharedArray, Zi)
Xj = convert(SharedArray, Xj)
Yj = convert(SharedArray, Yj)
Zj = convert(SharedArray, Zj)
L = convert(SharedArray, L)
Rot = convert(SharedArray, Rot)

这是相同的代码，但使用@parallel

function jt_transcoord_parallel(Xi, Yi, Zi, Xj, Yj, Zj, Rot, L)
r = SharedArray{Float64}(zeros((length(Xi),1)))
@parallel for i in 1:length(Xi)
    rxX = (Xj[i] - Xi[i]) / L[i]
    rxY = (Yj[i] - Yi[i]) / L[i]
    rxZ = (Zj[i] - Zi[i]) / L[i]
        if rxX == 0 && rxY == 0
            r[i] = [0 0 rxZ; cosd(Rot[i]) -rxZ*sind(Rot[i]) 0; sind(Rot[i]) rxZ*cosd(Rot[i]) 0]
        else
            R=sqrt(rxX^2+rxY^2)
            r21=(-rxX*rxZ*cosd(Rot[i])+rxY*sind(Rot[i]))/R
            r22=(-rxY*rxZ*cosd(Rot[i])-rxX*sind(Rot[i]))/R
            r23=R*cosd(Rot[i])
            r31=(rxX*rxZ*sind(Rot[i])+rxY*cosd(Rot[i]))/R
            r32=(rxY*rxZ*sind(Rot[i])-rxX*cosd(Rot[i]))/R
            r33=-R*sind(Rot[i])
            r[i] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]
        end
end
return r
end

我刚得到一个零向量。我的问题是：有没有办法在 Julia 中使用 @parallel 来实现这个函数并获得与我在原始函数中得到的结果相同的结果？

Answer 1

函数jt_transcoord 和jt_transcoord_parallel 存在重大编码缺陷。

在jt_transcoord 中，您将数组分配给向量元素位置。例如，您写r = Vector(length(Xi))，然后分配r[i] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]。但是r[i] 应该是一个数字，而你应该为它分配一个3x3 矩阵。我怀疑 Julia 正在悄悄地为你改变类型。

SharedArray 对象不会承认这种松散的类型转换行为。 SharedArray 的组件必须是单个primitive type，例如Float64，而Vector{Matrix} 不是原始类型。打开 Julia v0.6 REPL 并复制/粘贴以下代码：

r = SharedArray{Float64}(length(Xi))

for i in 1:length(Xi)
    rxX = (Xj[i] - Xi[i]) / L[i]
    rxY = (Yj[i] - Yi[i]) / L[i]
    rxZ = (Zj[i] - Zi[i]) / L[i]
        if rxX == 0 && rxY == 0
            r[i] = [0 0 rxZ; cosd(Rot[i]) -rxZ*sind(Rot[i]) 0; sind(Rot[i]) rxZ*cosd(Rot[i]) 0]
        else
            R    = sqrt(rxX^2+rxY^2)
            r21  = (-rxX*rxZ*cosd(Rot[i])+rxY*sind(Rot[i]))/R
            r22  = (-rxY*rxZ*cosd(Rot[i])-rxX*sind(Rot[i]))/R
            r23  = R*cosd(Rot[i])
            r31  = (rxX*rxZ*sind(Rot[i])+rxY*cosd(Rot[i]))/R
            r32  = (rxY*rxZ*sind(Rot[i])-rxX*cosd(Rot[i]))/R
            r33  = -R*sind(Rot[i])
            r[i] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]
        end
end

最后，我得到：

ERROR: MethodError: Cannot `convert` an object of type Array{Float64,2} to an object of type Float64
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] setindex!(::SharedArray{Float64,2}, ::Array{Float64,2}, ::Int64) at ./sharedarray.jl:483
 [2] macro expansion at ./REPL[26]:6 [inlined]
 [3] anonymous at ./<missing>:?

基本上，Julia 是在告诉您它不能将矩阵分配给 SharedArray 向量。

你有什么选择？

• 如果您坚持使用Vector{Matrix} 返回类型，则在jt_transcoord 中使用r = Vector{Matrix{Float64}}(length(Xi))。但是您不能为此使用 SharedArrays，因为 Vector{Matrix} 不是可接受的原始类型。
• 或者，如果您愿意使用张量（即 3 路数组）进行操作，则可以使用下面的伪代码 A。但是SharedArray 计算只有在您仔细考虑哪个进程拥有张量的哪个部分时才会对您有所帮助。否则，进程将需要相互通信，并且您的并行化函数可能执行得非常缓慢。
• 如果您愿意以 3n x 3 列方式排列 3x3 矩阵，则可以使用下面的伪代码 B。

伪代码A

function jt_transcoord_tensor(Xi, Yi, Zi, Xj, Yj, Zj, Rot, L)
    # initialize array
    r = Array{Float64}(3,3,length(Xi))
    # r = SharedArray{Float64,3}((3,3,length(Xi))) # for SharedArrays

    for i in 1:length(Xi)
    # @parallel for i in 1:length(Xi) # for SharedArrays
        # other code...
        r[:,:,i] = [0 0 rxZ; cosd(Rot[i]) -rxZ*sind(Rot[i]) 0; sind(Rot[i]) rxZ*cosd(Rot[i]) 0]
        # other code...
        r[:,:,i] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]
        end
    end
    return r
end

伪代码 B

function jt_transcoord_parallel(Xi, Yi, Zi, Xj, Yj, Zj, Rot, L)
    n = length(Xi)
    r = SharedArray{Float64}((3*n,3))
    @parallel for i in 1:length(Xi)
            # other code...
            r[(3*(i-1)+1):3*(i),:] = [0 0 rxZ; cosd(Rot[i]) -rxZ*sind(Rot[i]) 0; sind(Rot[i]) rxZ*cosd(Rot[i]) 0]
            # other code...
            r[(3*(i-1)+1):3*(i),:] = [rxX rxY rxZ;r21 r22 r23;r31 r32 r33]
        end
    end
    return r
end