JuliaでGPGPUして爆速でBuddhabrotを描き出す
はじめに
Scrapbox から失礼します。電気科4年生のToTTiです。
この記事ではプログラミング言語の Julia で GPGPU をする方法を紹介します。
実行環境は
OS: Windows10
Julia: 1.5.3
GPU: GeForce GTX 1050 Ti
です。
追記(2020/12/13): Buddhabrot のアルゴリズムのセクションで
2. $ z_0=0として、漸化式$ z_{n+1} = z_{n}を$ nが決めた値になるまで計算する。
となっていたところを
2. $ z_0=0として、漸化式$ z_{n+1}=z_n^2+cを$ nが決めた値になるまで計算する。
と修正しました。ご指摘ありがとうございました。
Juliaとは？
「そもそも、Julia ってなんだ？」
という方のために、簡単に紹介させていただきます。
Julia は Python と同じ "動的言語" です。さらに "動的型付け" でもあります。
動的言語というのは、その記述の柔軟性の変わりに動作速度が遅い傾向にあります。ところが、Julia は LLVM コンパイラーを用いることによって十分な柔軟性を保ちながら C言語相当の実行速度を得ることができます。
その他、科学計算用の豊富な外部ライブラリ(パッケージ)や多重ディスパッチ、他言語 (例: C, Python, FORTRAN)のコードの呼び出しなどの便利機能が多くあります。私的には Python の numpy を踏襲した ブロードキャスト という機能が大好きです。
詳しくは HP から Documentation を覗いてみてください。 ちなみにインストールは HP の緑色のDownload v1.5.3をクリックしてご自身の環境にあったインストーラをダウンロードして実行すればOKです。とっても簡単。 実行環境としては対話形式の REPL やみんな大好き Jupyter 、その他 Atom などのエディタなんかもあります。
ご自由にどうぞ。
もっぱら、私は Julia をプログラミング言語というよりは 超高級電卓ソフト として捉えています。
GPGPUってなに？
「 GPGPU ってなんやねん」
という方のために GPGPU についてもほんのり説明します。
GPGPU とは、グラフィック処理用のデバイスである GPU をグラフィック処理以外に使っちゃおうということです。たとえば、「RTX3080 を2枚使って 機械学習を行う」なんかは GPGPU のよくある例です。
CUDA.jlを使う
低レベルな領域で GPGPU を行うならnVidia製のGPUに使える CUDA C/C++ や OpenCL などがありますが、わざわざメモリ管理の確保や解放を記述したがるマゾはいないでしょう。
そこで今回は、Julia のパッケージ CUDA.jl を使ってなんともコメントの困る微妙な性能の GTX1050Ti を馬車馬のように働かせてみたいと思います。
そのために、まずはNVIDIAドライバーが導入します。NVIDIAのGPUを使っていれば基本的にはインストールされているはずですが、以下の手順でつまずいた場合などにはココからダウンロードしてインストールしてみてください。 ここまでくれば後は簡単、Julia の REPL (Julia 1.5.3.exe) を起動して]add CUDA と入力して Enter キーを押せばオッケーです。
ちなみに、CUDA.jlは名前の通り CUDA というものを利用しているため、CUDA が対応していない GPU は扱えません。ナンテコッタイ。対応している GPU はココから確認できます。 GPGPUプログラミング
CUDA.jl を使った GPGPU プログラミングには2通りの方法があります。
ブロードキャストを利用して配列同士の演算だけで記述する方法
関数をGPUに実行させる方法
まずは簡単な前者から。Julia ではブロードキャストによって、変数を入れる箇所に無理やり配列を入れることができます。
例えば
code: test1.jl
# 普通の関数の使い方
foge(2)
このように、関数名と()の間に.を入れることによって配列の要素それぞれを関数の引数に指定できます。
もちろん、ブロードキャストは演算子にも適用することができて、例えば
code: test2.jl
C = A .+ B
というように書くことができます。便利ですね。
この記法をそのままに GPGPU することができます。
では、GPUで計算するにはどうすればいいのか。
普通、Juliaの配列はArrayという型です。これはCPU用の配列なのでこれを GPU が使うことはできません。そこで、CuArrayという型を持ったGPU専用の配列を用意してこれを使ってやります。つまり
code: test3.jl
using CUDA # パッケージの読み込み
# CuArray型の配列を確保
C = A .+ B
と書けちゃうわけです。
ただ、この記法は配列の要素によって処理が変わるようなものを実装したい場合には使いにくいです。
そこで、今度は関数を記述する方法について紹介します。
ということで早速コードを紹介します。
code: test4.jl
using CUDA
# 関数を記述
function parallel_add(X, Y)
index = threadIdx().x
stride = blockDim().x
for i in index:stride:length(X)
end
return nothing
end
A = CUDA.zeros(10) # CuArray型の 0 が 10 個の配列
B = CUDA.rand(10) # CuArray型の 乱数 が 10 個の配列
@cuda threads=10 parallel_add(A, B)
println(Array(A))
上のコードで一番需要なのは @cuda threads=10, parallel_add(A, B) です。threads=10 によって関数 parallel_add を 10 個のGPUのコア(thread)で実行させています。
続いて、関数 parallel_add(X, Y) の中身について見る前に、GPUがどのように thread を管理しているかについて説明します。
GPUは複数の thread を block にまとめています。例えば、
GPU「block1 の thread1 さん、これやっといて」
GPU「block210 の thread500さんはこれね」
GPU「block2371 の ... 」
みたいに。(実際に指示を出すのはCPUなのですが...)
それを表したのが下の図
https://juliagpu.github.io/CUDA.jl/stable/tutorials/intro1.png
１つのブロックで並列に処理できる thraed の数は決まっているので、thread と block の数は適切に選んであげたほうが効率的です。
ただ、今回は簡単のためにblockの数を1つとしましたので @cuda threads=10, block=1 parallel_add(A, B) の blocks=1 は省略してます。
ではparallel_add(A, B) の中身を見ていきます。
threadIdx().xというのは実行している thread の id を返す関数です。
blockDim().x は１つの block に入っている thread の数を返す関数です。
先ほどの図にも似たようなものが書いてあるかと思います。これによって、配列 X, Y にアクセスするインデックスをスレッドごとに割り振ってあげて、並列処理を行います。
1:stride:length(X) は Python の range(1, length(X), stride)と一緒です。1 から length(X) までの数を stride ずつ飛ばして列挙します。
これによって、thread が少ない場合でも配列 X および Y の全ての要素に対して処理するようにします。
これ以降は大丈夫でしょう。
Buddhabrot のアルゴリズム
CUDA.jl によるGPGPUの方法も説明したところで早速、本題のコードをお見せしたいのですが、そもそも、何をやっているのかを説明しないと何が何だか分からないと思いますので、そもそものアルゴリズムについて説明します。
今回、ご紹介するアルゴリズムは Buddhabrot という、Mandelbrot 集合の特殊な可視化方法の一つです。
ざっくりと流れを説明します。
1. ランダムに複素数$ c をとる。
2. $ z_0=0として、漸化式$ z_{n+1}=z_n^2+cを$ nが決めた値になるまで計算する。
3. $ z_nの絶対値が途中で$ 2以上になった場合は、そこまでの$ z_nを複素平面上にプロットする。
4. 1~3 を思う存分 (10万回~100万回程度) 繰り返す。
ざっとこんな感じです。繰り返せば繰り返すほどきれいな画像が得られるのですが、その分計算量が大変なことになって、並みの CPU ではあまり太刀打ちできません。
そこで GPU の出番というわけです。
上記の手順の 1~3 は10個並列にやろうが100個並列にやろうが影響がないので、GPGPU にピッタリです。
Buddhabrot を CUAD.jl で描く
ということでお待たせしました。
本題のコードがこちらになります。
code:buddha.jl
using Random, CUDA
function gpu_buddha(Grid, c_array, max_iterate::Int64)
index = threadIdx().x
stride = blockDim().x
width, height = size(Grid)
for i in index:stride:length(c_array)
# パラメータの初期化
z = 0f0 + 0f0im
is_out = false
out_iterate = 1
# zが発散するかどうかのチェック
for j in 1:max_iterate
z = z*z + c
if !(real(z)^2 + imag(z)^2 < 4f0)
is_out = true
out_iterate = j
break
end
end
# 発散するzの軌道を記録
if is_out
z = 0f0 + 0f0im
for j in 1:out_iterate-1
z = z*z + c
# Grid 上に対応する z のインデックスに 1 足す
m = round(Int64, (2f0+real(z))*(width-1)/4f0 + 1f0)
n = round(Int64, (2f0+imag(z))*(height-1)/4f0 + 1f0)
end
end
end
return nothing
end
# 描画設定
Grid = CUDA.zeros(Int32, 512, 512) # 画像サイズ
max_iterate = 2^10 # 最大でどれくらい漸化式を計算するか
loop_iterate = 2^14 # どれくらい点を記録するか
for t in 1:loop_iterate
# 本計算
c_array = CUDA.rand(Complex{Float32}, 256) .* 4 .- (2f0 + 2f0im)
@cuda threads=256 gpu_buddha(Grid, c_array, max_iterate)
end
synchronize()
本体は関数 gpu_buddha です。コイツで$ -2 - 2iから$ 2 + 2iまでの正方領域上の複素平面に対応した Grid に$ z_nをプロットして画像にします。
引数 max_iterate が Buddhabrot のアルゴリズムの「$ nが決めた値になるまで計算する 」の決められた値です。
この値によって得られる画像がかなり変化して面白いです。
gpu_buddha 内は基本的に Float32 で書いています。というのも、GPU は CPU と比べると高度な計算能力を持っていないため、Float64 で計算すると遅いのです。
@inbounds は配列のインデックスがその配列をはみ出していないかのチェックをしないようにするマクロです。
後は Grid を Plots の heatmap などで画像にすると次のようなものが得られます。
https://gyazo.com/5082d109331be74fada4158576d9ae61
どうでしょう。ざらつきは多いですがなかなか面白い画像じゃないでしょうか。
複素数$ cの取り方をちょっと制限したり、適切なノイズ除去を行うことによって、もう少し綺麗な画像を得ることもできるのですが
それはまた別の話。
まとめ
そんなわけで、Julia で GPGPU する方法について紹介しました。
結構お手軽にできるんじゃないかなと思います。
Julia は本当に科学計算との相性がいいので、流体のシミュレーションや機械学習とかを GPGPU で行うのには本当に向いていると思います。
面白いなと思った方は是非、プログラミング言語 兼 超高級電卓 Julia で遊んでみてください。
明日は kagitaharuyo さんが 簡易電気回路シュミレータの作成記録 について書いてくれるそうです。楽しみですね。