2015-02-06 7 views
5

Sono interessato all'utilizzo di Julia SharedArray s per un progetto di calcolo scientifico. La mia attuale implementazione fa appello a BLAS per tutte le operazioni con la matrice vettoriale, ma ho pensato che forse uno SharedArray avrebbe offerto una certa velocità su macchine multicore. La mia idea è semplicemente aggiornare un vettore di output indice per indice, coltivando gli aggiornamenti dell'indice ai processi di lavoro.BLAS v. Aggiornamenti paralleli per oggetti Julia SharedArray

discussioni precedenti here su SharedArray s e here su oggetti di memoria condivisa non offrono indicazioni chiare su questo problema. Sembra intuitivamente abbastanza semplice, ma dopo i test, sono un po 'confuso sul perché questo approccio funzioni così male (vedi il codice sotto). Per i principianti, sembra che @parallel for allochi molta memoria. E se prefisso il ciclo con @sync, che sembra una cosa intelligente da fare se l'intero vettore di output è richiesto in seguito, il ciclo parallelo è sostanzialmente più lento (anche se senza @sync, il loop è molto veloce).

Ho interpretato erroneamente l'uso corretto dell'oggetto SharedArray? O forse ho assegnato in modo inefficiente i calcoli?

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

uscita dal test con 4 operai + nodo 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 

risposta

4

si sta eseguendo in diverse questioni, di cui la più importante è che crea un Xt[:,i] (memoria allocazione) nuovo array. Ecco una dimostrazione che si avvicina a ciò che si vuole:

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 

Risultati (quando partire 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 

Per confronto, se in esecuzione su questa stessa macchina:

julia> blas_set_num_threads(8) 

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

Quindi l'implementazione raw di Julia è almeno competitiva con BLAS.

+0

ooh, non avevo pensato alle allocazioni di memoria per 'Xt [:, i]'. E il tuo codice di esempio è facile da seguire. Grazie! –

+0

Perché 'dotcol' è più veloce quando' B [i, j] 'è usato (come nel tuo codice) piuttosto che' B [j, i] '(senza richiedere che la matrice' X' sia trasposta)? – Taiki

+2

@Taiki è questo un problema di accesso alla memoria a grandi passi? Julia memorizza gli array in formato colonna principale. Vedere la [documentazione] (http://docs.julialang.org/en/release-0.4/manual/performance-tips/#access-arrays-in-memory-order-along-columns) per i dettagli. –