モンテカルロ法で円周率の(並行)計算

先日Matlabで書いたものと同じアルゴリズムを、Juliaで並行計算できるようにしました。

tic()

N = 1000000;
K = 1000;
est_pi_total = @parallel (+) for k=1:K
  x = rand(N);
  y = rand(N);
  4*( sum(x.*x + y.*y < 1) / N)
end
printf("K=%d N=%d ==> pi=%f\n", K, N, est_pi_total / K);

toc()

MacBook Pro 13 (Mid 2009)はコア数が2個なので、使用コアが1つの時と2つの時で比較してみました。まずは標準の1コア版。

$ julia
julia> load("montecarlo_pi.jl")
K=1000 N=1000000 ==> pi=3.141652
elapsed time: 54.607524156570435 seconds

そしてオプションで使用するコアを2に指定します。

$ julia -p 2
julia> load("montecarlo_pi.jl")
K=1000 N=1000000 ==> pi=3.141499
elapsed time: 43.42981815338135 seconds

若干、速くなったかな。


【2017-08-04 追記】

上記コードがJulia 0.6で動かなかったので、以下のように修正しました。

N = 1000000;
K = 1000;
@time est_pi_total = @parallel (+) for k=1:K
  x = rand(N);
  y = rand(N);
  4.0 * ( sum(x.*x + y.*y .< 1.0) / N)
end
@printf("K=%d N=%d ==> pi=%f\n", K, N, est_pi_total / K);

起動時オプション「julia -p auto」で、手持ちのCPUコアを全て使ってくれます。ローカルファイルの読み込みと実行は「load」ではなく「include」を使います。MacBook Pro 15(Late 2013, 2.6 GHz Core i7)だと、1コア1スレッドで13.5秒だったのが、4コア8スレッドを使って6.0秒になりました。