Juliaでも並列計算ができる、と
https://docs.julialang.org/en/v0.6.4/manual/parallel-computing/#Parallel-Computing-1
とある。
スレッド並列については良記事
"Juliaのスレッド並列を使って簡単な並列計算を楽しむ"
https://qiita.com/goropikari/items/dc9ceddfb2d92a744e25
があったので、この記事ではプロセス並列(MPI並列のようなもの)を行う。
並列化経験について
私: Fortran90でスパコンでMPI+openmpでブンブン回している。
バージョン
Julia 0.6.4
(08/03/18追記)
@goropikari さんにコメントをいただいた。その結果、もっとコンパクトに書ける事がわかった。最後に追記しておく。
最終結論だけここに先に書いておくと、forループで回すindexをxとして、1からnmaxまで回した結果を配列として持ちたいときは、
A = pmap(x -> func(x), 1:nmax )
とすれば、サイズnmaxの配列Aが得られる。
(08/06/18追記)
バッチジョブなどでコントロールされたクラスターマシン上では、
julia --machinefile=$PBS_NODEFILE /path/to/your/home/dir/test_julia.jl
とすれば、割り当てられたCPUコアを使った並列計算が動く。実際にこのやり方でやると並列化が効き計算時間が減っていくことを確認した。ここで、-pオプションは使わないことに注意。
参考にしたURL:
https://discourse.julialang.org/t/help-setting-up-julia-on-a-cluster/5519/5
(04/28/20追記)
上のmachinefileはJulia 1.0以降では変更されていて、
julia --machine-file=$PBS_NODEFILE /path/to/your/home/dir/test_julia.jl
となることに注意。
コード
forループを分割して、分割したそれぞれを各CPUコアに割り当てて計算させるコードがあれば、かなり広範囲の種類の並列化ができる。
試行錯誤して得られたコードを先に載せておく。これが最適かどうかはわからない。
@everywhere function test(nm,nmax)
ssum = zeros(Float64,nm)
v = procs()
nv = length(v)-1
if nv == 0
nv = 1
end
a = Array{Any}(nv)
if nmax % nv != 0
println("Error nmax % nv should be 0. But, nmax= $(nmax) and nv = $(nv)")
end
for i in 1:nv
starti = 1+(i-1)*nmax/nv
endi = starti + nmax/nv -1
a[i] = @spawn svdsum(nm, starti,endi)
end
println("After nv loop")
for i in 1:nv
ssum += fetch(a[i])
end
return ssum
end
println(nprocs())
@time test(400,4000)
実行したいときは、
julia -p 4 test.jl
とすれば、4並列で動く。
forループで毎回処理したい関数は、ここでは、svdsumというものにしている。
svdsumでは、startiからendiまで、forループを実行する。このstartiやendiがプロセスごとに異なることで、全体のループを分割して実行することが可能になる。
@everywhere function svdsum(nm,starti,endi)
ssum = zeros(Float64,nm)
println("Start!")
for i = starti:endi
u,s,v = svd((3*i).*eye(nm)+cos(i)*ones(nm,nm))
ssum += s
end
println("End!")
return ssum
end
結局、@spawnを何回も繰り返して呼べば、次々と仕事を確保したプロセス分やってくれる。
そして、それを配列aに保存しておいて、forループを回りきった後に、改めてfetchでデータを回収する。
これによって、それぞれのプロセスが自分の仕事をこなし、データが最後に集まることになる。
この方法で、使うプロセスを倍に増やすと(物理CPUコアの数よりプロセス数を超えない場合に)、計算時間は半分になることを確かめた。
なお、@everywhereをつけないと1番以外のプロセスが関数を読み込んでくれないので、@everywhereは必要。
MPIと比べて超簡単。
(08/03/18追記)
@goropikari さんにコメントをいただき、さらに簡単に書ける事がわかった。
functionのtest部分は
nm = 100
nmax = 1000
@time ssum = sum( pmap(x -> svdi(nm,x), 1:nmax ) )
println(sum(ssum))
だけで良いようだ。ここで、要素iを引数にするfunctionとして
@everywhere function svdi(nm,i)
u,s,v = svd((3*i).*eye(nm)+cos(i)*ones(nm,nm))
return s
end
を定義した。この書き方を使えば、forループで様々な処理をしたあとに、それぞれの結果を配列に格納することも可能で、その場合には、
nm = 100
nmax = 1000
@time smat =pmap(x -> svdi2(nm,x), 1:nmax )
println(smat)
とする。svdi2は
@everywhere function svdi2(nm,i)
u,s,v = svd((3*i).*eye(nm)+cos(i)*ones(nm,nm))
return sum(s)
end
としている。