OpenCLでマンデルブロ集合を書いてみた

元ネタ　http://codezine.jp/article/detail/5439
この記事ではRadeonだが、僕のノートについているのはGeForceM8400GSなので、nVidiaの開発者ページから、
GPU Computing SDK code samples and more
を落としてきてインストール
ドライバはCUDA対応のものが必要だが、最近は標準でCUDA対応になっているようなので割愛。
VisualC++2010Expressをインストール。
includeパスと、ライブラリパスを加えて、以下のプログラムをコンパイルした。

#include 
#include 
#include 

#define __CL_ENABLE_EXCEPTIONS
#include 

using namespace std;

int main(){
  try{
    vector platforms;
    cl::Platform::get(&platforms);
    cl::Platform platform = platforms[0];
    cl_context_properties prop[3] = {CL_CONTEXT_PLATFORM,reinterpret_cast(platform()),0 };
    cl::Context  context(CL_DEVICE_TYPE_GPU, prop);
    vector devices = 
      context.getInfo();
    cl::Device device = devices[0];
    cl::CommandQueue queue(context,device);
    cl::Buffer buffer(context,CL_MEM_READ_ONLY, (cl_int)(sizeof(int)*128*128));

    static string src_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(src_str.c_str(),src_str.size()));
    cl::Program program(context,sources);
    program.build(devices);

    cl::Kernel kernel(program,"mandelbrot");
    kernel.setArg(0,10000);
    kernel.setArg(1,-2.0f);
    kernel.setArg(2,1.5f);
    kernel.setArg(3,0.0234375f);
    kernel.setArg(4,-0.0234375f);
    kernel.setArg(5,buffer);
    int* result = new int[128*128];
    queue.enqueueNDRangeKernel(kernel,cl::NullRange,cl::NDRange(128,128),cl::NullRange);
    queue.enqueueReadBuffer(buffer,CL_TRUE,0,sizeof(int)*128*128,result);
    
    for(int i = 0;i < 128;i++){
      for(int j = 0;j < 128;j++){
        cout << result[j+i*128] << ",";
      }
      cout << endl;
    }
 
  }catch(const cl::Error& ex){
  }
}

汚いプログラムだけどご勘弁。結果はコンソールにCSVで出力されるので、それをファイルに書き出してExcelでグラフにしてみた。