Julia@Parallel与SharedArray的并行计算

Question

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

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]]

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)

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矩阵。我怀疑朱莉娅在悄悄地为你改变类型。

ShareDarRay对象将不允许这种松散的类型转换行为。ShareDarRay的组件必须是单个基元类型，如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向量。

你有什么选择？

null

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