ジュリア集合をOpenCLで実装

ジュリア集合自体は、以下の本に記載されている。

CUDA by Example 汎用GPUプログラミング入門

そして、それをOpenCLに移植してみよう。参考にするのは、例のフィックスターズの本:

改訂新版 OpenCL入門 1.2対応 マルチコアCPU・GPUのための並列プログラミング

github.com

とりあえず実装はできたものの、OpenCLについてまだ理解できていないところが多く、微妙な実装になっている。

• OpenCLで画像データを取り扱うには、もっと良い方法があるはず (今はCUDA by Exampleのglutのライブラリを利用している)
• CUDAコアのサイズについてはあまり考慮せず作ってある

カーネルの実装

CUDAの例では、カーネルの記述の中に構造体(実際にはクラス)が記述してあるが、OpenCLのカーネルではそれが記述できず、 しょうがないので、全てクラスの部分を配列と関数に置き換えた。

float magnitude2(float x, float y)
{
        return x*x + y*y;
}


void array_mult(float a[2], float b[2], float c[2])
{
        c[0] = a[0] * b[0] - a[1] * b[1];
        c[1] = a[1] * b[0] + a[0] * b[1];
        return;
}

void array_add(float a[2], float b[2], float c[2])
{
        c[0] = a[0] + b[0];
        c[1] = a[1] + b[1];
        return;
}


int julia(int x, int y)
{
        const float scale = 1.5;
        float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);
        float jy = scale * (float)(DIM / 2 - y) / (DIM / 2);

        float c[2] = { -0.8, 0.156 };
        float a[2] = { jx, jy };

        int i = 0;
        for (i = 0; i < 200; i++) {
                // a = a * a + c;
                float tmp[2];
                array_mult(a, a, tmp);
                array_add(tmp, c, a);
                if (magnitude2(a[0], a[1]) > 1000) {
                        return 0;
                }
        }
        return 1;
}

しかし、それ以外は基本的にCUDAの場合と一緒だ。

__kernel void julia_kernel(__global unsigned char *ptr)
{
        int x = get_group_id(0);
        int y = get_local_id(0);

        int offset = y * DIM + x;

        int juliaValue = julia(x, y);

        ptr[offset * 4 + 0] = 255 * juliaValue;
        ptr[offset * 4 + 1] = 0;
        ptr[offset * 4 + 2] = 0;
        ptr[offset * 4 + 3] = 255;
}

ホストの実装

ホストは基本的なOpenCLの呼び出しを利用している。引数として、入力値はまったく設定しておらず、出力値だけReadBufferで読み込むようにしている。

        /* Generate Program */
        program = clCreateProgramWithSource(context, 1, (const char **)&source_str, (const size_t *)&source_size, &ret);
        ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
        kernel = clCreateKernel(program, "julia_kernel", &ret);
        printf("err:%d\n", ret);

        if (ret != CL_SUCCESS) {
                size_t len;
                clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, NULL, NULL, &len);
                char *log = new char[len]; //or whatever you use
                clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, len, log, NULL);
                printf("Error was occurred.%d\n", len);
                printf("%s\n", log);
                exit(EXIT_FAILURE);
        }

        ret = clSetKernelArg(kernel, 0, sizeof(memobj), (void *)&memobj);
...
        CPUBitmap bitmap(DIM, DIM);
        ret = clEnqueueReadBuffer(command_queue, memobj, CL_TRUE, 0, DIM * DIM * 4, bitmap.get_ptr(), 0, NULL, NULL);
        uint8_t *mem = bitmap.get_ptr();

とりあえずジュリア集合は表示されるようになったが、OpenCLの呼び出しと、CUDAコアがどのように割り付けられて実行されているのかいまいち関係がまだ掴めていない。 ちょっとずつ探っていこう。 今使っているGPUはCUDAコアは2880個なので、1ピクセル1個ずつ処理してるわけじゃないんだよな。global_work_itemとlocal_work_itemにより、どのように配置されるのかが いまいちまだ分かっていない。。。

www.nvidia.co.jp