プロセス並列化＋リダクション演算による配列の和の計算 [Julia]

配列要素の総和を求めるような演算は，数値計算においてしばしば表れる．たとえば Julia で書かれた次のようなコードがあるとする．

# 配列 a, b の確保
sum_a = zero(eltype(a));
sum_b = zero(eltype(b));
for j in 1:N
    for i in 1:M
        sum_a += a(i,j);
        sum_b += b(i,j);
    end
end

このようなコードを並列化する場合，リソース（メモリ）の競合を避けるための対策が必要となる．たとえば，C++/OpenMP でスレッド並列化する場合は，reduction指示子があるので比較的簡単に並列化できる．

// 配列 a, b の確保
double sum_a = 0.0;
double sum_b = 0.0;
#pragma omp parallel for reduction(+:sum_a, sum_b)
for ( int i = 0; i < M; i++) {
    for ( int j = 0; j < N; j++ ) {
        sum_a += a[i][j];
        sum_b += b[i][j];
    }
}

残念ながら，現時点で Julia のスレッド並列にはリダクションが実装されていない．一時変数やアトミック操作を使う手もあるが，なるべく簡単に並列化したい．幸いプロセス並列にはリダクションが実装されているので，上のコードは次のようにすることでプロセス並列化できる．

using Distributed
addprocs(); # 複数プロセス立ち上げ
# 配列 a, b の確保
sum_a, sum_a = @distributed (+) for j in 1:N
    sum_a = zero(eltype(a));
    sum_b = zero(eltype(b));
    for i in 1:M
        sum_a += a(i,j);
        sum_b += b(i,j);
    end
    [sum_a, sum_b];
end