を正月休みにしていた。CUDAに関しては雑な知識しかなかったので、『CUDA BY EXAMPLE』なる書籍を読んだ。『CUDA C プロフェッショナルプログラミング』なる本が名著であると勧められてはいたが、まずモチベーションを殺さないために短めのものを一読してから重い本で知識をアップデートするやり方に慣れていたので、手元にあった薄い方の本を手にとった次第だ。
CUDAの文法や基礎知識は解説できるほど理解していないので、本や何やで確認してほしい。
mandelbrot
恐らく最も並列化が容易でビジュアル的にもよいであろうマンデルブロ集合が最初の目標としては最適だろうと考えた。初めてのCUDAプログラムなので、__shared__
メモリだとか__constant__
メモリだとかそういうのは気にせずまったりやることにする。
naiveな実装
std::complex
をCUDAで使えるのかわからなかったので、ひどい再発明から入った。実際のところ、__device__
関数で何も考えずにstd::complex z, c;
z * z + c;
みたいなのを書くと、__host__
関数を__device__
からは呼べないと怒られるし、std::complex
のコンストラクタはconstexpr
なので--expt-relaxed-constexpr
オプションをつけないとそこでも怒られる。
とりあえず必要なものだけ定義する。__device__
が付いていること以外は通常の関数と変わらない。template
な__device__
関数の例は見たことがなかったが(CUDA C
の資料を見ているので当然である)、できるだろうと思って書いたらできた。CUDAのC++対応がどの程度進んでいるかまとまっているところはないだろうか。
template<typename realT> struct Complex { realT r, i; }; template<typename realT> __device__ __forceinline__ Complex<realT> add(const Complex<realT> a, const Complex<realT> b) { return Complex<realT>{a.r + b.r, a.i + b.i}; } template<typename realT> __device__ __forceinline__ Complex<realT> mul(const Complex<realT> a, const Complex<realT> b) { return Complex<realT>{a.r * b.r - a.i * b.i, a.r * b.i + b.r * a.i}; } template<typename realT> __device__ __forceinline__ Complex<realT> mandelbrot_step(const Complex<realT> z, const Complex<realT> c) { return add(mul(z, z), c); } template<typename realT> __device__ __forceinline__ realT abs(const Complex<realT> z) { return z.r * z.r + z.i * z.i; }
まあ、素直にthrust::complex
を使えば済む話なので、これを読んでいる人はこんなクソッタレな再発明をするべきではない。
で、各block
の各thread
でマンデルブロ集合内であるかを判定する関数を書く。threadIdx
とblockIdx
で担当する場所を決め、それぞれが判定する。
template<typename realT> __device__ bool is_mandelbrot(const Complex<realT> z, const std::size_t threshold) { Complex<realT> x{0., 0.}; for(std::size_t i=0; i<threshold; ++i) { x = mandelbrot_step(x, z); if(abs(x) > 2.) return false; } return true; } template<typename realT> __global__ void write_mandelbrot( bool* result, const realT rmin, const realT rmax, const realT imin, const realT imax, const std::size_t th, const std::size_t x_dim, const std::size_t y_dim) { const std::size_t x = threadIdx.x + blockIdx.x * blockDim.x; const std::size_t y = threadIdx.y + blockIdx.y * blockDim.y; const std::size_t offset = x + gridDim.x * blockDim.x * y; if(offset > x_dim * y_dim) return; const realT r = rmin + (rmax - rmin) / x_dim * x; const realT i = imin + (imax - imin) / y_dim * y; const Complex<realT> z{r, i}; result[offset] = is_mandelbrot(z, th); return; }
これを__host__
側から呼ぶ。template
な__global__
関数を呼ぶときはwrite_mandelbrot<double><<<blocks, threads>>>(args...)
のようにすればできる。
#include <png++/png.hpp> /*上記の関数達*/ int main() { const double rmin = -2.0; const double rmax = 1.0; const double imin = -1.5; const double imax = 1.5; const std::size_t threshold = 400; const std::size_t x_dim = 1024; const std::size_t y_dim = 1024; bool* mandelbrot; cudaError_t ermalloc = cudaMalloc((void**)&mandelbrot, x_dim * y_dim); assert(ermalloc == 0); dim3 blocks(x_dim/32, y_dim/32); dim3 threads(32, 32); write_mandelbrot<double><<<blocks, threads>>>( mandelbrot, rmin, rmax, imin, imax, threshold, x_dim, y_dim); bool *result = new bool[x_dim * y_dim]; cudaError_t ercpy = cudaMemcpy( result, mandelbrot, x_dim * y_dim, cudaMemcpyDeviceToHost); assert(ercpy == 0); cudaFree(mandelbrot); png::image<png::rgb_pixel> image(x_dim, y_dim); for(std::size_t i=0; i<x_dim; ++i) { for(std::size_t j=0; j<y_dim; ++j) { std::size_t offset = i + x_dim * j; if(result[offset]) image[i][j] = png::rgb_pixel(255, 0, 0); else image[i][j] = png::rgb_pixel(0,0,0); } } image.write("mandelbrot.png"); return 0; }
少々雑だが、これで動く。
thrust
さて、C++erとしては生のCUDAコードを書くよりも、thrust
の使用を検討したほうがよいという話は聞く。生CUDAを書いてしまってから検討するのは遅きに失している気もするが、thrust
を使ってみよう。
実はだいぶ以前に一度CUDAを使ってみようとthrust
のクイックスタートガイドを和訳したことがあるので(著作権などがどうなっているか厳しい気がしたので公開していないが)、thrust
に関してその程度の知識はある。結局その後今に至るまで実際にコードを書くことはなかったのだが。
というわけで、thrust::transform
によって、複素平面上の点を格納したthrust::device_vector<thrust::complex<double>>
を、マンデルブロ集合内かどうかの判定結果が格納されたthrust::device_vector<bool>
に変換することにしよう。
その判定用の関数オブジェクトを作る。なんのことはない、thrust::complex
を受け取ってthreshold
回のイテレーションの間に絶対値が2より大きくなるか調べるだけだ。
template<typename realT> struct judge_mandelbrot { const std::size_t threshold; judge_mandelbrot(const std::size_t th): threshold(th){} __host__ __device__ bool operator()(const thrust::complex<realT>& z) const { thrust::complex<realT> x(0., 0.); for(std::size_t i=0; i<threshold; ++i) { x = x * x + z; if(thrust::abs(x) > 2.) return false; } return true; } };
最初にthrust::device_vector<thrust::complex<double>>
を複素平面で埋めるやり方がCUDA初心者にはよくわからなかったので、諦めてhost_vector
を作ってコピーすることにした。
// generate complex plane thrust::host_vector<thrust::complex<double>> host_plane(size); for(std::size_t i=0; i<r_resolution; ++i) { for(std::size_t j=0; j<i_resolution; ++j) { const thrust::complex<double> val = min + thrust::complex<double>(dr * i, di * j); const std::size_t offset = i + r_resolution * j; host_plane[offset] = val; } } thrust::device_vector<thrust::complex<double>> complex_plane = host_plane;
で、これをthrust::transform
する。
judge_mandelbrot<double> judge(threshold); thrust::device_vector<bool> is_mandelbrot(size); thrust::transform(complex_plane.begin(), complex_plane.end(), is_mandelbrot.begin(), judge);
これで計算は終了だ。こんなに簡単でいいのだろうか。あとは出力するだけである。
ところで、出力するときにいちいちdevice_vector
にアクセスしていると、一度host_vector
にコピーしたときと比較して、信じられないくらい時間がかかることがわかった。要素一つごとにcudaMemcpy
しているのだろう。
thrust
つきのプログラムをコンパイルする方法に一瞬迷ったが、
$ nvcc -O3 -std=c++11 mandelbrot_thrust.cu -o mandelbrot -lpng12
でよいようだ。最後のはpng++ライブラリのリンク。
ところで、雑でもいいのでパフォーマンスを計測してみようとwrite_mandelbrot
関数にかかった時間とthrust::transform
にかかった時間をstd::chrono
で計測してみたのだが、素人が書いたwrite_mandelbrot
で5回の平均が18[us]、対してthrust
は平均10[us]であった。thrust
強すぎでは? 以前にArrayfireの紹介記事で素直な実装だと生CUDAコードを書いてもArrayfireの方が早いというのがあったが、thrustも(その記事とは比較にならないほど雑なコードだろうが)素直な生CUDAよりも早いとは。
この内容は
GitHub - ToruNiina/cuda_practice
で公開されている。
次回はisingモデルを書いたときのことをまとめる。