CUDA：「超並列コンピューティング」はじめました（その2）
～高速化TipsとOpenGLでのグラフィクスへの応用

5年ぶりのマンデルブロ集合

　ゴキゲンなビデオ・カードが手に入ったことだし、最後にグラフィクスへのちょっとした応用を試してみましょうか、マンデルブロ集合の計算と描画です。このお題は5年ほど前、AMD製GPUを使ってOpenCLと.NETで実装しました。これを元ネタに、CUDAとOpenGLで再実装します。

　マンデルブロ集合は複素数の集合です。とある複素数pに対し複素数z（初期値 0+0i）を用意して、z = z * z + pを無限に繰り返してもzが発散しないとき、pはマンデルブロ集合の要素です。無限に繰り返すのではいつまで待っても結果が得られないので、繰り返し回数の上限：limitを与え、点pが発散を始める（zのノルムが4を超える）までの繰り返し回数を求めます。

list06

// 点pがマンデルブロ集合の要素であるかを調べる。
// limitを上限として、発散を始めるまでの演算回数を返す
//   名前空間cuにコーディングを楽にするためのヘルパ関数を用意した
__host__ __device__ unsigned int calculate(cu::fcomplex p, unsigned int limit) {
  cu::fcomplex z = cu::make_fcomplex(0.0f, 0.0f);
  unsigned int count;
  for ( count = 0U; count < limit && cu::norm(z) <= 4.0f; ++count ) {
    z = z*z + p;
  }
  return count;
}

　2つの複素数：begin（左下隅）とend（右上隅）で定まる矩形を碁盤目状に分割し、マス目の各点に対して上記calculateを呼んでマンデルブロ集合をoutに求めます。

list07

// 複素平面上のbegin(左下隅)とend(右上隅)で定まる矩形領域をsizeで刻み、
// 各点に対して上記calculateを呼んでマンデルブロ集合を求める
__global__ void kernel_mandelbrot(cu::fcomplex begin, cu::fcomplex end, uint2 size,
                                  unsigned int limit, unsigned int* out) {
  cu::fcomplex diff  = end - begin;
  cu::fcomplex step  = cu::make_fcomplex(cu::real(diff)/size.x, cu::imag(diff)/size.y);
  dim3 pos = blockDim * blockIdx + threadIdx;
  if ( pos.x < size.x && pos.y < size.y ) {
    cu::fcomplex p = begin + cu::make_fcomplex(cu::real(step)*pos.x, cu::imag(step)*pos.y);
    unsigned int count = calculate(p, limit);
    cu::ary2mtx<unsigned int>(out, size.x)[pos.y][pos.x] = count;

  }
}

　さらにもう一つ、得られたunsigned int値をRGBA（4byte）に変換します。0～limitの範囲を0～255に正規化した値をR,G,Bに、A（アルファ値）は255にします（α-brendを使ってないので不要なんですけどね）。

list08

// 各点のカウント値をグレイ・スケールRGBAに変換する
__global__ void kernel_count2RGBA(unsigned int* in, uchar4* out, unsigned int size, unsigned int limit, unsigned int bias) {
  unsigned int i = cu::width(blockDim * blockIdx + threadIdx);
  if ( i < size ) {
    unsigned int val = in[i];
    val = (unsigned int)(256 * rsqrtf((float)val / limit));
    val = (val + bias) % 256;
    out[i] = make_uchar4(val,val,val,255);
  }
}

　左下隅（-2.0 -1.5i）、右上隅（1.0 +1.5i）でマンデルブロ集合を求め、RGBA変換すると：

list09

void device_mandelbrot(unsigned int limit, uint2 size, unsigned int* buf, uchar4* out) {
  // 開始点(begin)/終了点(end) 
  cu::fcomplex begin = cu::make_fcomplex(-2.0f, -1.5f);
  cu::fcomplex end   = cu::make_fcomplex( 1.0f,  1.5f);
 
  // グレイ・スケール変換時のバイアス
  static unsigned int bias = 0;
  bias = ++bias % 256;

  // マンデルブロ集合を求め、
  {
  dim3 desired(size.x, size.y);
  dim3 block(16U,16U); // 16*16 = 256 threads per block 
  dim3 grid = cu::roundup_div(desired, block);
  kernel_mandelbrot<<<grid, block>>>(begin, end, size, limit, buf);
  }
  // グレイ・スケールRGBAに変換する
  {
  unsigned int desired = size.x * size.y;
  unsigned int block = 256;
  unsigned int grid = (desired + block - 1) / block;
  kernel_count2RGBA<<<grid, block>>>(buf, reinterpret_cast<uchar4*>(buf), desired, limit, bias);
  }
  cudaMemcpy(out, buf, size.x*size.y*sizeof(uchar4), cudaMemcpyDeviceToHost);
}

　RGBA列ができれば、こいつをOpenGLのテクスチャとしてそのまま使えます。OpenGLコードは今回本スジではないので割愛、僕のトモダチ西山さんのkindle本『チュートリアル形式で始めるOpenGL[2D編]』にあったサンプルをほとんどそのまま使わせていただきました（西山さんありがとです）。

　こんなのが現れましたよ♪

　計算／テクスチャ変換／描画コミコミで毎秒400フレーム処理できてます。計算とテクスチャ変換をホスト（CPU）で行った場合6FPS足らずでしたから、65倍速くなってます。GeForce700番台のローエンドGPU：GT720の乗ったビデオ・カードをアキバで手に入れ同じコードを実行したところ140FPSでした。GT720はCUDAコア192個積んでいてGTX750の3/8、FPS比はCUDAコア数の比とほぼ一致しました（偶然かしら?）。

　CUDA ToolkitのインストールでVisual Studioにビルトインされるデバッガ／プロファイラ：Nsightで実行中のタイムラインを観察すると、

　……あちゃー、赤く囲った部分がデバイス内の処理本体（kernel）、デバイスからホストへのデータ転送（cudaMemcpy）が思いっきり足引っ張ってますしくしく。

　実はCUDAにはOpenGL（とDirect3D）が利用しているビデオ・メモリ側の領域を直接読み書きできる手段が用意されていて、それを使えば時間のかかるホスト-デバイス間のデータ転送を端折って格段に速くなるらしいのですが……なんせOpenGLは触り始めてわずか3日、シロートにはいささか手強い代物です。NVIDIAのCUDAエンジニア：森野さんに助けていただきました。

　概略は以下のとおり。

　まずOpenGL側にPBO（Pixel Buffer Object）を用意します。CUDAが書いてOpenGLが読む領域ですね。

// Pixel Buffer Objectを作成: CUDAからの出力が書き込まれる。
GLuint pbo_id;
glGenBuffers(1, &pbo_id);

　つぎにCUDAが提供するOpenGL Interoperability-APIを用い、OpenGL側の上記PBOとCUDA側のデバイス・メモリとの仲介役となる cudaGraphicResource_tを作ります。この部分はかなりコスト高なのでモロモロの初期化時に済ませておきます。

// OpenGL interoperability
// PBOは、CUDAから、Write-onlyでアクセスされる。
cudaGraphicsResource_t cudaPBO = nullptr;
cudaGraphicsGLRegisterBuffer(&cudaPBO, pbo_id, cudaGraphicsRegisterFlagsWriteDiscard);

　これで準備完了。

　CUDAによる処理（kernel呼び出し）に先立ってPBOをデバイス・メモリにマッピング（紐づけ）を行い、マップされた領域のポインタと大きさを手に入れます。

// PBOをCUDAにマッピング
cudaGraphicsMapResources(1, &cudaPBO);
// マップされた領域のポインタ（と大きさ）を取得
size_t pboSize = 0;
uchar4 *d_pbo = nullptr;
cudaGraphicsResourceGetMappedPointer((void**)&d_pbo, &pboSize, cudaPBO);

　手に入れたポインタと大きさを引数にkernelを実行し、cudaGraphicsUnmapResources()でマップを解除した時点で結果がOpenGL側のPBOに書き込まれています。あとはOpenGLにおまかせ、PBOからテクスチャを作って描画。

　コードの全容はサンプルファイルを眺めていただくとして、さて少しは速くなったでしょうか。

　650FPS！ 1フレームを1.5msで処理しています、すばらしい。CPU（single-thread）と比べると114倍、ふたケタ違いのパフォーマンス。CUDAすげー。