BLAS v. Równoległe aktualizacje dla obiektów Julia SharedArray

Question

Jestem zainteresowany wykorzystaniem Julia SharedArray s dla projektu naukowego informatyki. Moja obecna implementacja odwołuje się do BLAS dla wszystkich operacji macierzy-wektorowych, ale pomyślałem, że być może SharedArray zaoferuje pewne przyspieszenie na maszynach wielordzeniowych. Moim pomysłem jest po prostu zaktualizować indeks wyjściowy po indeksie, aktualizując indeksy do procesów roboczych.

Poprzednie dyskusje here o SharedArray s oraz here o obiektach pamięci współdzielonej nie zawierały jasnych wskazówek dotyczących tego problemu. Wydaje się to intuicyjnie proste, ale po testach jestem nieco zdezorientowany, dlaczego to podejście działa tak źle (zobacz poniższy kod). Na początek wygląda na to, że @parallel for przydziela dużo pamięci. I jeśli przedwzmacniam pętlę z @sync, co wydaje się mądrą rzeczą, jeśli cały wektor wyjściowy jest wymagany później, to pętla równoległa jest znacznie wolniejsza (ale bez @sync, pętla jest potężna szybko).

Czy źle zinterpretowałem właściwe użycie obiektu SharedArray? A może nieefektywnie przydzieliłem obliczenia?

### test for speed gain w/ SharedArray vs. Array ### 

# problem dimensions 
n = 10000; p = 25000 

# set BLAS threads; 64 seems reasonable in testing 
blas_set_num_threads(64) 

# make normal Arrays 
x = randn(n,p) 
y = ones(p) 
z = zeros(n) 

# make SharedArrays 
X = convert(SharedArray{Float64,2}, x) 
Y = convert(SharedArray{Float64,1}, y) 
Z = convert(SharedArray{Float64,1}, z) 

# run BLAS.gemv! on Arrays twice, time second case 
BLAS.gemv!('N', 1.0, x, y, 0.0, z) 
@time BLAS.gemv!('N', 1.0, x, y, 0.0, z) 

# does BLAS work equally well for SharedArrays? 
# check timing result and ensure same answer 
BLAS.gemv!('N', 1.0, X, Y, 0.0, Z) 
@time BLAS.gemv!('N', 1.0, X, Y, 0.0, Z) 
println("$(isequal(z,Z))") # should be true 

# SharedArrays can be updated in parallel 
# code a loop to farm updates to worker nodes 
# use transposed X to place rows of X in columnar format 
# should (hopefully) help with performance issues from stride 
Xt = X' 
@parallel for i = 1:n 
    Z[i] = dot(Y, Xt[:,i]) 
end 

# now time the synchronized copy of this 
@time @sync @parallel for i = 1:n 
    Z[i] = dot(Y, Xt[:,i]) 
end 

# still get same result? 
println("$(isequal(z,Z))") # should be true 

Wyjście z testu z 4 pracowników + węzłów 1 master:

elapsed time: 0.109010169 seconds (80 bytes allocated) 
elapsed time: 0.110858551 seconds (80 bytes allocated) 
true 
elapsed time: 1.726231048 seconds (119936 bytes allocated) 
true 

Answer 1

używasz do kilku zagadnień, z których najważniejszym jest to, że Xt[:,i] tworzy nową tablicę (przydzielania pamięci). Oto wykazanie, że dostaje się bliżej, co chcesz:

n = 10000; p = 25000 

# make normal Arrays 
x = randn(n,p) 
y = ones(p) 
z = zeros(n) 

# make SharedArrays 
X = convert(SharedArray, x) 
Y = convert(SharedArray, y) 
Z = convert(SharedArray, z) 

Xt = X' 

@everywhere function dotcol(a, B, j) 
    length(a) == size(B,1) || throw(DimensionMismatch("a and B must have the same number of rows")) 
    s = 0.0 
    @inbounds @simd for i = 1:length(a) 
     s += a[i]*B[i,j] 
    end 
    s 
end 

function run1!(Z, Y, Xt) 
    for j = 1:size(Xt, 2) 
     Z[j] = dotcol(Y, Xt, j) 
    end 
    Z 
end 

function runp!(Z, Y, Xt) 
    @sync @parallel for j = 1:size(Xt, 2) 
     Z[j] = dotcol(Y, Xt, j) 
    end 
    Z 
end 

run1!(Z, Y, Xt) 
runp!(Z, Y, Xt) 
@time run1!(Z, Y, Xt) 
zc = copy(sdata(Z)) 
fill!(Z, -1) 
@time runp!(Z, Y, Xt) 

@show sdata(Z) == zc 

Wyniki (gdy zaczynające julia -p 8):

julia> include("/tmp/paralleldot.jl") 
elapsed time: 0.465755791 seconds (80 bytes allocated) 
elapsed time: 0.076751406 seconds (282 kB allocated) 
sdata(Z) == zc = true 

Dla porównania, gdy działa na tej samej maszynie:

julia> blas_set_num_threads(8) 

julia> @time A_mul_B!(Z, X, Y); 
elapsed time: 0.067611858 seconds (80 bytes allocated) 

Tak więc surowa implementacja Julii jest co najmniej konkurencyjna w stosunku do BLAS.