.NETでマンデルブロ集合を描く（番外編）

OpenCLで超並列コンピューティング

Buffer

　今回実装するマンデルブロ集合描画アプリケーションは幅×高さが512×512のピクセル平面を複素平面に見立て、512×512=26万数千個の複素数C=x+yiに対し:

int calculate(const float x, const float y, const int limit) {
    int   count = 0;
    // Z = 0 + 0i                                             
    float prevX = 0.0f;                                       
    float prevY = 0.0f;                                       
    for (int i = 0; i < limit; i++) {                         
        // Z = Z^2 + C を |Z| > 2.0 となるまで繰り返す 
        count = i + 1;                                        
        float newX = (prevX * prevX) - (prevY * prevY) + x;   
        float newY = (2.0f * prevX * prevY) + y;              
        if (((newX * newX) + (newY * newY)) > 4.0f) {         
            break;                                            
        } else {                                              
            prevX = newX;                                     
            prevY = newY;                                     
        }                                                     
    }                                                         
    return count;                                             
}                                                             

を呼び出し、得られた値に応じた色でピクセルを着色します。この処理は与えられた引数のみに依存し、他に影響を及ぼすことがないので多数の点について同時に行うことができます。小さなCPUを多数持つGPUにはうってつけといえます。

　計算結果を保持する512×512×sizeof(int)の領域をDevice内に保持するBuffer（バッファ）を用意します。

Buffer buffer(context, CL_MEM_READ_ONLY, (cl_int)(sizeof(int)*512*512));

Program

　次にProgram（プログラム）。ProgramはDevice上で実行されるコードを集めたライブラリとして機能します。Device上で実行されるコードはOpenCL Cで記述し、コード文字列をビルドします。

tatic string source_str =
  "int calculate(const float x, const float y, const int limit) { \n"
  "    int   count = 0;                                           \n"
  "    float prevX = 0.0f;                                        \n"
  "    float prevY = 0.0f;                                        \n"
  "    for (int i = 0; i < limit; i++) {                          \n"
  "        count = i + 1;                                         \n"
  "        float newX = (prevX * prevX) - (prevY * prevY) + x;    \n"
  "        float newY = (2.0f * prevX * prevY) + y;               \n"
  "        if (((newX * newX) + (newY * newY)) > 4.0f) {          \n"
  "            break;                                             \n"
  "        } else {                                               \n"
  "            prevX = newX;                                      \n"
  "            prevY = newY;                                      \n"
  "        }                                                      \n"
  "    }                                                          \n"
  "    return count;                                              \n"
  "}                                                              \n"
  "                                                               \n"
  "__kernel void mandelbrot(                                      \n"
  "                const int max_calc,                            \n"
  "                const float start_x,                           \n"
  "                const float start_y,                           \n"
  "                const float step_x,                            \n"
  "                const float step_y,                            \n"
  "               __global int *out) {                            \n"
  "    size_t x     = get_global_id(0);                           \n"
  "    size_t y     = get_global_id(1);                           \n"
  "    size_t width = get_global_size(0);                         \n"
  "                                                               \n"
  "    int index = (width * y) + x;                               \n"
  "                                                               \n"
  "    float px = mad(step_x, (float)x, start_x);                 \n"
  "    float py = mad(step_y, (float)y, start_y);                 \n"
  "                                                               \n"
  "    out[index] = calculate(px, py, max_calc);                  \n"
  "}                                                              \n";
  cl::Program::Sources sources;
  sources.push_back(make_pair(source_str.c_str(), source_str.size()));
  cl::Program program(context, sources);
  program.build(devices);

Kernelと処理の実行/結果の読み出し

　Kernel（カーネル）はProgram内の関数から選ばれた関数ひとつを指し、これがDevice上での処理の開始点を表すことになります。

cl::Kernel kernel(program,"mandelbrot");

　kernelに引数を与え、CommandQueueに登録することでDevice内での処理が行われます。計算結果をbufferから読み出す処理もCommandQueueに与えておきましょう。bufferはデバイス内にありますから結果の読み出しにはそれと同じ（かそれ以上）の領域をホスト側に確保しておかなくてはなりません。

int* result = new int[512*512];
complex<float> start; // ピクセル(0,0)に対応する複素数点
complex<float> strp;  // 1ピクセルの（複素平面での）幅と高さ
kernel.setArg(0, 512);
kernel.setArg(1, start.real());
kernel.setArg(2, start.imag());
kernel.setArg(3, step.real());
kernel.setArg(4, step.imag());
kernel.setArg(5, buffer);
queue.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(512,512), cl::NullRange);
queue.enqueueReadBuffer(buffer, CL_TRUE, 0, sizeof(int)*512*512, result);

　これで処理終了。resultに求まった値を基にピクセルを着色してできあがりとなります。

アプリケーション

　アプリケーションは『.NETでマンデルブロ集合を描く（後日談）』をベースにしました。フォームはC++/CLIで実装し、OpenCLを使った処理本体はC++でDLLを作り、C++/CLIから呼び出します。

　実行結果はご覧のとおり。

　マンデルブロ集合の計算そのものに要する時間は8ms程度、native-C++とPPLを用いたdual-coreでの計算時間が100msでしたからさらにその10倍以上のパフォーマンスを叩き出しています。今回使用したRadeon HD5670はミドルレンジエントリーモデルです。ハイエンドのビデオカードであればピクセルシェーダー数が1000基を超えますから、さらに高速な処理が期待できます。