D言語で始めるOpenCL(1) 導入編

目次

• D言語で始めるOpenCL(1) 導入編(今回)
• D言語で始めるOpenCL(2) ローカルメモリー編
• D言語で始めるOpenCL(3) 下準備＆行列パディング編
• D言語で始めるOpenCL(4) プライベートメモリー編

ソースコード

なぜD言語?

タイトルの通り、D言語でOpenCLやります。行列計算などを速くやってみたいと思います。

なぜわざわざD言語みたいなマイナー言語でやるかというと、一番の理由は好きだからですが、他にいくつか理由を挙げることは可能です。

• C言語のAPIを簡単に扱える。
• ポインタ操作なども含めたメモリの直接操作が簡単に行える。
• 上記の利点を阻害することなく、厳密な型による制約や実行時のチェックなど、豊富なチェック機能を好きなだけ活用できる。
• 言語としての表現力が高く、GPU側ソースコード生成などの文字列操作が簡単に行える。
• OpenCLの環境構築が比較的簡単に行える。

なぜOpenCL?

GPGPUのデファクトスタンダードとも言えるNVIDIAのGPUを使う限り、どうもCUDAを使った方が速く、しかもライブラリなどの資産も豊富なようです。

その上で、OpenCLを使用した場合の利点としては、下記が考えられると思います。

AMDなど他社のGPUでも動く

以前、NVIDIAがGeForceのデータセンターでの使用を禁じるよう利用規約を変更して騒動になりました。
ソフトウェアがNVIDIAのみにロックインされていると、同様の予期しない事態が発生して色々困る可能性はどうしてもあります。

OpenCLであれば、基本的にはGPUベンダーに依存しないので、別GPUに乗り換えることは比較的容易です……少なくとも不可能ではないはずです。

FPGAなどのGPUとは異なるプロセッサでも動かせるかもしれない

OpenCLは、GPGPU限定ではなく、FPGAやDSPで実行することも考えられている仕様です。
書いたコードがそれらの異なるターゲットで再利用できるかもしれません。

でも普通は、アーキテクチャに合わせた最適化をしないといけないから、そのままというのはいくらなんでも無理そう……。

それでも、GPGPUが苦手な問題の場合、FPGA等との組み合わせで対応できる可能性はあります。

Intel HDなど専用GPUのない環境でもGPGPUの恩恵に預かれる

これが結構大きい利点かもしれません。貧乏人にも優しいOpenCL。

CUDA（とNVIDIAのGPU）がないと動かせないプログラムは、気軽に手元のノートPCで試したり、会社のデスクトップPCでこっそり動かすことが難しいです。
しかしOpenCLのプログラムであれば、そういったマシンでも、ある程度実用的な速度で動く見込みがあります。
たとえIntel HDであっても、GPUを使うとCPUの何倍も速かったりします。(最適化をちゃんとやれば)

環境構築

D言語のインストール

標準的なD言語のコンパイラであるdmdとdubが使えます。まずはこれらをインストールします。

この記事を読むような人はもうインストールされていると思いますが……。

OpenCL ランタイムライブラリの用意

Intel HDしかない環境の場合

下記からOpenCL用のドライバやランタイムライブラリがダウンロードできるそうです。

Macの場合

標準でOpenCLが入っているようです。

私のMacbook Air(Mid 2011)はなぜかIntel HD 3000を認識してくれませんが……。（古すぎる）

ちゃんとしたGPUを持っている場合

羨ましいですね。

GPUベンダに従って頑張ってドライバとライブラリをインストールしてください。

OpenCL仕様

OpenCLのAPI仕様およびGPU側C/C++仕様は下記にまとまっています。

この記事ではOpenCL 1.2の仕様を参照したいと思います。

なんでそんな古い仕様を見るのか？

私のメインマシンがMacbook Air(Mid 2011)で1.2までしか対応していないからです……。
上位互換性はあるので、より新しい環境でもとりあえず動きます。

OpenCLの構成要素

OpenCLは、異種混合システムのための並列プログラミングのオープンな標準仕様です。
よって、異種混合するために各種計算資源を色々管理できるようになっています。

OpenCLを使うためには、そんな計算資源を検索して初期化する(そして後処理する)コードをたくさん書く必要があります。

私が使った範囲での主な登場人物は下記の通りです。

• Platform
  • GPUベンダーによるランタイムライブラリなどのプラットフォーム。
• Device
  • Platform上で使用できる計算デバイス。
• Context
  • Deviceを使用するコンテキスト。複数のデバイスを同時に使用できる。
• CommandQueue
  • Deviceに対してKernelの実行やBufferへの読み書きを命じるためのキュー。
• Buffer
  • Device上で使用できるメモリオブジェクト。
• Program
  • Device上で実行できるプログラム。
• Kernel
  • Device上で実行できる関数。Programの中で定義されていて、ホストマシンから引数等を指定して呼び出せる。

他にOpenGLにあるようなImageやSamplerもあるみたいですが、とりあえず割愛します。

簡単なラッパーライブラリ

たとえD言語であっても素のOpenCL C APIを叩き続けるのは辛いので、ごく薄いラッパーライブラリを用意しました。

まだまだ機能は足りませんが、適宜追加していくことにします。
基本的に行なっていることは下記の通りです。

• エラーコードをチェックしてエラー発生時はメッセージ付き例外(OpenClException)を飛ばす
• パラメーターの省略やポインタの配列化
• よく使いそうなフラグ・オプション指定でのAPI呼び出しを、独立した関数とする
  • デバイス情報の取得など

依存ライブラリは今のところDerelictCLのみです。

HelloWorld

とりあえずOpenCLでHello Worldしてみます。OpenCL Cではなぜかprintfが使えて標準出力にメッセージが出せるので、それを使ってみます。

なお、ソースコードはラッパーライブラリのexamplesに含まれています。

OpenCLのランタイムライブラリがインストールされていれば、git cloneしてdub run --config=examples-hello-worldするだけで、おそらくきっとすぐに動かせます。

ソースコード

import std.stdio : writefln;

// OpenCLのラッパーライブラリ等のimport
import cl = dcltk;
import derelict.opencl.cl : cl_event;

/// とりあえずのメイン関数
void main() {
    // OpenCLプラットフォームの取得。とりあえず最初のプラットフォームをロードしている。
    auto platformId = cl.loadOpenCl();

    // プラットフォーム情報の表示
    writefln("loaded: %s %s %s(%s) %s [%s]",
        cl.getPlatformProfile(platformId),     // プロファイル情報
        cl.getPlatformName(platformId),        // プラットフォーム名
        cl.getPlatformVersion(platformId),     // バージョン
        cl.getPlatformCLVersion(platformId),   // OpenCLバージョン(enum CLVersion型)
        cl.getPlatformVendor(platformId),      // ベンダー名
        cl.getPlatformExtensions(platformId)); // 対応している拡張機能

    // 存在するデバイスを列挙する
    auto deviceIds = cl.getAllDeviceIds(platformId);
    writefln("devices:");
    foreach(i, d; deviceIds) {
        // デバイス情報の表示
        writefln("%d: %s %s gmem: %d, lmem: %d, cu: %d, w: %d, dim: %d %s",
                i,                                     // 番号
                cl.getDeviceName(d),                   // デバイス名
                cl.getDeviceVersion(d),                // バージョン
                cl.getDeviceGlobalMemorySize(d),       // グローバルメモリサイズ
                cl.getDeviceLocalMemorySize(d),        // ローカルメモリサイズ
                cl.getDeviceMaxComputeUnits(d),        // 最大計算ユニット数
                cl.getDeviceMaxWorkGroupSize(d),       // 最大ワークグループ数
                cl.getDeviceMaxWorkItemDimensions(d),  // 最大ワークアイテム次元数
                cl.getDeviceMaxWorkItemSizes(d));      // 各次元の最大ワークアイテム数
    }

    // デフォルトデバイス向けのOpenCLコンテキスト生成
    auto context = cl.createDefaultContext(platformId);
    scope(exit) cl.releaseContext(context);

    // 最初のデバイス向けのコマンドキュー生成
    auto device = deviceIds[0];
    auto commandQueue = cl.createCommandQueue(context, device);
    scope(exit) cl.releaseCommandQueue(commandQueue);

    // プログラムをソースコード文字列から生成
    auto program = cl.createProgramFromSource(context, `
        __kernel void helloWorld(void) {
            printf("Hello,World!\n");
        }
    `);
    scope(exit) cl.releaseProgram(program);

    // プログラムをビルド
    cl.buildProgram(program, deviceIds);

    // プログラムに含まれるカーネル関数からカーネルを生成
    auto kernel = cl.createKernel(program, "helloWorld");
    scope(exit) cl.releaseKernel(kernel);

    // コマンドキューにカーネルの実行を投げる
    cl_event event;
    cl.enqueueKernel(commandQueue, kernel, [1], [1], event);
    cl.flushCommandQueue(commandQueue);

    // コマンド終了までイベントを待機する。
    cl.waitAndReleaseEvents(event);
    cl.finishCommandQueue(commandQueue);
}

実行結果

% dub run --config=examples-hello-world
Performing "debug" build using dmd for x86_64.
derelict-util 3.0.0-beta.2: target for configuration "library" is up to date.
derelict-cl 3.2.0: target for configuration "library" is up to date.
dcltk ~master: building configuration "examples-hello-world"...
Linking...
To force a rebuild of up-to-date targets, run again with --force.
Running ./dcltk 
loaded: FULL_PROFILE Apple OpenCL 1.2 (May 24 2018 20:07:03)(CL12) Apple [cl_APPLE_SetMemObjectDestructor cl_APPLE_ContextLoggingFunctions cl_APPLE_clut cl_APPLE_query_kernel_names cl_APPLE_gl_sharing cl_khr_gl_event]
devices:
0: Intel(R) Core(TM) i7-2677M CPU @ 1.80GHz OpenCL C 1.2  gmem: 4294967296, lmem: 32768, cu: 4, w: 1024, dim: 3 [1024, 1, 1]
Hello,World!

計算する

GPGPUでHello, Worldしていても虚しいだけなので、得意の計算をやらせてみます。

とりあえず行列の積を計算させてみます。これさえできれば大抵のことはできます。

ソースコード

import std.datetime.stopwatch : benchmark;
import std.math : approxEqual, sqrt, floor, ceil;
import std.random : uniform01;
import std.stdio : writefln;

import cl = dcltk;
import derelict.opencl.cl : cl_event;

/// テスト用にCPUで行列積を計算する関数も用意する。
void productCpu(
        const(float)[] lhs,
        const(float)[] rhs,
        float[] result,
        uint rows,
        uint cols,
        uint resultCols)
in {
    assert(lhs.length == rows * cols);
    assert(rhs.length == cols * resultCols);
    assert(result.length == rows * resultCols);
} body {
    // 色々なものを犠牲にした実装
    for(size_t i = 0; i < rows; ++i) {
        for(size_t j = 0; j < resultCols; ++j) {
            float value = 0.0f;
            for(size_t k = 0; k < cols; ++k) {
                value += lhs[i * cols + k] * rhs[k * resultCols + j];
            }
            result[i * resultCols + j] = value;
        }
    }
}

/// 念のため適当にユニットテスト
unittest {
    immutable(float)[] lhs = [1, 2, 3, 4];
    immutable(float)[] rhs = [5, 6, 7, 8];
    auto result = new float[2 * 2];

    productCpu(lhs, rhs, result, 2, 2, 2);

    assert(approxEqual(1 * 5 + 2 * 7, result[0 * 2 + 0]));
    assert(approxEqual(1 * 6 + 2 * 8, result[0 * 2 + 1]));
    assert(approxEqual(3 * 5 + 4 * 7, result[1 * 2 + 0]));
    assert(approxEqual(3 * 6 + 4 * 8, result[1 * 2 + 1]));

    auto resultScalar = new float[1];
    productCpu(lhs, rhs, resultScalar, 1, 4, 1);
    assert(approxEqual(1 * 5 + 2 * 6 + 3 * 7 + 4 * 8, resultScalar[0]));
}

/// ここからが本体
void main() {
    // 行列サイズ。100行200列 × 200行300列 = 100行300列
    enum {
        ROWS = 100,
        COLS = 200,
        RESULT_COLS = 300
    }

    // 両辺の行列を[0, 1)の一様乱数で適当に初期化する。
    auto lhs = new float[ROWS * COLS];
    foreach(ref e; lhs) {
        e = uniform01!float;
    }
    auto rhs = new float[COLS * RESULT_COLS];
    foreach(ref e; rhs) {
        e = uniform01!float;
    }

    // 結果格納用の配列を用意する。
    auto cpuResult = new float[ROWS * RESULT_COLS];
    auto gpuResult = new float[ROWS * RESULT_COLS];

    // OpenCL プラットフォーム初期化
    auto platformId = cl.loadOpenCl();
    writefln("loaded: %s %s %s(%s) %s [%s]",
        cl.getPlatformProfile(platformId),
        cl.getPlatformName(platformId),
        cl.getPlatformVersion(platformId),
        cl.getPlatformCLVersion(platformId),
        cl.getPlatformVendor(platformId),
        cl.getPlatformExtensions(platformId));

    // デフォルトのコンテキストを生成
    auto context = cl.createDefaultContext(platformId);
    scope(exit) cl.releaseContext(context);

    // デバイス取得。無ければエラー。
    auto deviceIds = cl.getContextDeviceIds(context);
    if(deviceIds.length <= 0) {
        throw new cl.OpenClException("device not found.");
    }

    // 最初のデバイスを使用する。関連情報を表示する。
    auto device = deviceIds[0];
    writefln("device: %s %s gmem: %d, lmem: %d, cu: %d, w: %d, dim: %d %s",
        cl.getDeviceName(device),
        cl.getDeviceVersion(device),
        cl.getDeviceGlobalMemorySize(device),
        cl.getDeviceLocalMemorySize(device),
        cl.getDeviceMaxComputeUnits(device),
        cl.getDeviceMaxWorkGroupSize(device),
        cl.getDeviceMaxWorkItemDimensions(device),
        cl.getDeviceMaxWorkItemSizes(device));
    if(!cl.isGpuDevice(device)) {
        // GPUでなければ警告表示
        writefln("WARNING: device is not GPU!");
    }

    // コマンドキュー生成
    auto commandQueue = cl.createCommandQueue(context, device);
    scope(exit) cl.releaseCommandQueue(commandQueue);

    // GPUで実行するプログラムをソースコードから生成する。
    auto program = cl.createProgramFromSource(context, `
        __kernel void product(
                __global const float *lhs,
                __global const float *rhs,
                __global float *result,
                uint rows,
                uint cols,
                uint resultCols) {
            for(size_t i = 0; i < rows; ++i) {
                for(size_t j = 0; j < resultCols; ++j) {
                    float value = 0.0f;
                    for(size_t k = 0; k < cols; ++k) {
                        value += lhs[i * cols + k] * rhs[k * resultCols + j];
                    }
                    result[i * resultCols + j] = value;
                }
            }
        }
    `);
    scope(exit) cl.releaseProgram(program);

    // プログラムをビルド
    cl.buildProgram(program, deviceIds);

    // ビルドしたプログラムからカーネルを生成
    auto kernel = cl.createKernel(program, "product");
    scope(exit) cl.releaseKernel(kernel);

    // カーネルに渡す引数のバッファを生成。
    auto lhsBuffer = cl.createReadBuffer(context, lhs);
    scope(exit) cl.releaseBuffer(lhsBuffer);
    auto rhsBuffer = cl.createReadBuffer(context, rhs);
    scope(exit) cl.releaseBuffer(rhsBuffer);

    // 結果格納用のバッファを生成。バイト単位のサイズのみ指定する。
    auto resultBuffer = cl.createWriteBuffer(context, gpuResult.length * float.sizeof);
    scope(exit) cl.releaseBuffer(resultBuffer);

    // カーネル引数を設定
    cl.setKernelArg(kernel, 0, lhsBuffer);
    cl.setKernelArg(kernel, 1, rhsBuffer);
    cl.setKernelArg(kernel, 2, resultBuffer);
    cl.setKernelArg(kernel, 3, ROWS);
    cl.setKernelArg(kernel, 4, COLS);
    cl.setKernelArg(kernel, 5, RESULT_COLS);

    // カーネルの最大ワークグループサイズ・推奨ワークグループサイズ係数を表示
    writefln("kernel w: %s, pw: %s",
        cl.getKernelWorkGroupSize(kernel, device),
        cl.getKernelPreferredWorkGroupSizeMultiple(kernel, device));

    // ベンチマーク対象となるGPU実行関数
    void productGpu() {
        cl_event event;
        cl.enqueueKernel(commandQueue, kernel, [1024], [1]);
        cl.enqueueReadBuffer(commandQueue, resultBuffer, 0, gpuResult, event);
        cl.flushCommandQueue(commandQueue);
        cl.waitAndReleaseEvents(event);
        cl.finishCommandQueue(commandQueue);
    }

    // CPUとGPUでそれぞれベンチマークを行う。
    immutable cpuMsecs = benchmark!(() => productCpu(
                lhs, rhs, cpuResult, ROWS, COLS, RESULT_COLS))(1)[0].total!"msecs";
    immutable gpuMsecs = benchmark!(() => productGpu())(1)[0].total!"msecs";
    writefln("cpu: %d msecs, gpu: %d msecs", cpuMsecs, gpuMsecs);

    //writefln("%s", cpuResult);
    //writefln("%s", gpuResult);

    // CPUの計算結果とGPUの計算結果を比較する。差があればエラー。
    foreach(i, e; cpuResult) {
        assert(approxEqual(e, gpuResult[i]));
    }
}

実行結果

実装ができたので早速動かしてみます。

GPUマシンとしてとりあえずGoogle Compute Engineで適当なインスタンスを作り、NVIDIA Tesla K80を挿して実行します。
見せてあげよう、GPUの雷を!

loaded: FULL_PROFILE NVIDIA CUDA OpenCL 1.2 CUDA 9.2.176(CL12) NVIDIA Corporation [cl_khr_global_int32_base_atomics
 cl_khr_global_int32_extended_atomics cl_khr_local_int32_base_atomics cl_khr_local_int32_extended_atomics cl_khr_fp
64 cl_khr_byte_addressable_store cl_khr_icd cl_khr_gl_sharing cl_nv_compiler_options cl_nv_device_attribute_query c
l_nv_pragma_unroll cl_nv_copy_opts cl_nv_create_buffer]
device: Tesla K80 OpenCL C 1.2  gmem: 11996954624, lmem: 49152, cu: 13, w: 1024, dim: 3 [1024, 1024, 64]
kernel w: 1024, pw: 32
cpu: 37 msecs, gpu: 711 msecs

うわっ…私のK80、遅すぎ…？

CPUと同じようなことをやらせた結果、CPUの20倍遅いという結果が出ました……。

この結果はやばいので、最低限GPUらしく動作するよう修正します。

OpenCLでの並列計算

先ほどの実装で、よりによってCPUよりも遅くなってしまったのは、1コア・1スレッドでしか計算していないのが原因です。

GPUの細かいアーキテクチャはよくわかりませんが、多数のコアを並列に動作させて速度を出しているらしいことは知っています。
それなのに1コアだけ使うようなコードを書いていたら、つまり普通のCPUを動かすのと同じようなコードでは、遅くて当たり前です。
CPUよりパワーのないコアをワンオペさせる深夜のブラック牛丼屋のようなものです。

ここから、OpenCLで多数のコアをどうすれば働かせられるのか、その一番基本的なところをやってみます。

ワークアイテムとワークグループ

OpenCLでは、計算を実行する単位のことをワークアイテムと呼びます。
これは普通のCPUにおけるスレッドのようなものらしいです。
GPUなどでは、このワークアイテムがたくさん同時に実行でき、結果的に速く計算できます。
逆に言えば、ワークアイテムをたくさん同時に実行させないと遅いです。

ワークアイテムの数を設定するのは、enqueueKernel(OpenCL APIではclEnqueueNDRangeKernel)です。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [1], // グローバルワークアイテム数
    [1]  // ローカルワークアイテム数
);

上記でコメントにより示した箇所で、カーネルを実行するワークアイテム数を指定できます。
ここでワークアイテム数を増やせば良いのですが、いくつにするべきでしょうか。
必要なパラメーターは、getDeviceInfo系の関数で取得できます。

// 計算ユニット数。K80では13。
cl.getDeviceMaxComputeUnits(device);

// カーネルの推奨ワークグループサイズの係数。
// つまりワークグループサイズはこの数の倍数にしろということ。
// K80では32。
cl.getKernelPreferredWorkGroupSizeMultiple(kernel, device);

// 各次元の最大ワークアイテム数。K80では[1024, 1024, 64](3次元)。
cl.getDeviceMaxWorkItemSizes(device);

数字の意味がよくわからないですね。私もよくわかりません。でもなんとか説明してみます。

ワークアイテムは空間に広げられる

ワークアイテムは、先ほども言ったとおり計算を実行する単位です。それを複数走らせることができます。
ワークアイテムは複数あるので、それぞれに0から始まるIDが付与されます。
付与されるのですが、それが実は1次元だけではなく2次元、3次元といったベクトル値とすることができます。

行列や画像といった空間的な広がりのあるデータを処理しやすいように、そういった機能が盛り込まれているようです。
GPUが何次元までのワークアイテムを扱えるかは、getDeviceMaxWorkItemDimensionsで取得できます。(基本的には3次元までのようですが……)

例として、8つのワークアイテムを3次元に展開して処理を行わせるようenqueueKernelを呼び出してみます。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [2, 2, 2], // グローバルワークアイテム数
    [2, 2, 2]  // ローカルワークアイテム数
);

これで、実行されるワークアイテムの数が2 * 2 * 2 = 8になり、カーネルが8つ並列で実行されることになります。
ところで、引数に2つワークアイテム数があるのが気になりますよね……。

ワークグループでワークアイテムを分割する

OpenCLでは、展開したワークアイテムをワークグループとしてまとめることができます。
ワークグループとしてまとめると何が良いかというと、以下の機能があるそうです。

• ワークグループ内でのみ共有できる高速なメモリーがある。うまく使うと速い。
• 1つのワークグループは1つの計算ユニット(Compute Unit)で実行される。つまりワークグループを分けておくと複数の計算ユニットを同時に使用でき、速い。

つまり頑張ると速くなるということです。

さて、ではワークアイテムをワークグループに分割するにはどうすればよいでしょうか。
そのためには、enqueueKernelの引数を再び修正することになります。
ワークアイテム数は8のままとし、それらをワークグループに分けてみます。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [2, 2, 2], // グローバルワークアイテム数
    [1, 1, 1]  // ローカルワークアイテム数
);

上記の指定でワークグループが8つになります。どういうことだ。

ワークグループ数は、グローバルワークアイテム数 ÷ ローカルワークアイテム数 で計算されます。

今回は3次元で指定しているので、体積で考えるとわかりやすいです。

全体のワークアイテムの体積は、2 * 2 * 2 = 8 です。(グローバルワークアイテム数)
いっぽう、ローカルのワークアイテム体積は、 1 * 1 * 1 = 1 です。(ローカルワークアイテム数)
全体のワークアイテムの体積をローカルのワークアイテム体積で埋めると、ローカルのワークアイテムがぴったり8つ入ることになります。

この、グローバルなワークアイテム体積にぴったりはまる、ローカルなワークアイテム体積のそれぞれが、1つのワークグループとなります。

ワークアイテムを2次元に広げる場合は、面積ということになります。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [2, 2], // グローバルワークアイテム数
    [1, 1]  // ローカルワークアイテム数
);

この場合、全体のワークアイテム数は 2 * 2 = 4、ワークグループごとのワークアイテム数は 1 * 1 = 1、ワークグループ数は 4 ÷ 1 = 4 になります。

なんとなく分かりましたね。きっとそうですね。
では早速ワークアイテムを増やしていきましょう。
計算ユニット数13ということなので、ワークグループが13になるようにします。(縁起の悪い数字だ……)
また、推奨ワークグループサイズ係数(getKernelPreferredWorkGroupSizeMultiple)が32ということで、ワークグループサイズは32にします。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [13 * 32], // グローバルワークアイテム数
    [32]  // ローカルワークアイテム数
);

実行結果(ワークアイテム増加版)

cpu: 36 msecs, gpu: 712 msecs

うわっ…私のK80、遅すぎ…？

さて……。

ワークアイテム別に処理をさせる

全然速くならなかったですね。

なぜかというと、肝心のカーネルのコードの方を全然いじっていないからです。

        __kernel void product(
                __global const float *lhs,
                __global const float *rhs,
                __global float *result,
                uint rows,
                uint cols,
                uint resultCols) {
            for(size_t i = 0; i < rows; ++i) {
                for(size_t j = 0; j < resultCols; ++j) {
                    float value = 0.0f;
                    for(size_t k = 0; k < cols; ++k) {
                        value += lhs[i * cols + k] * rhs[k * resultCols + j];
                    }
                    result[i * resultCols + j] = value;
                }
            }
        }

現状では、このコードが全ワークアイテムで同時実行されます。
同じ計算をするだけなので、結果は正しくなりますが、完全に無駄です。

どれだけ並列に計算処理が動いても、全てが同じ計算であれば無意味です。
意味のある並列処理をさせるには、それぞれが異なる処理を行う必要があります。

このコードを修正して、ワークアイテムごとに異なる処理をさせるようにします。

__kernel void product(
        __global const float *lhs,
        __global const float *rhs,
        __global float *result,
        uint rows,
        uint cols,
        uint resultCols) {
    // ワークアイテムのグローバルIDを取得する。
    const size_t groupIndex = get_global_id(0);
    const size_t groupSize = get_global_size(0);

    for(size_t i = groupIndex; i < rows; i += groupSize) {
        for(size_t j = 0; j < resultCols; ++j) {
            float value = 0.0f;
            for(size_t k = 0; k < cols; ++k) {
                value += lhs[i * rows + k] * rhs[k * resultCols + j];
            }
            result[i * resultCols + j] = value;
        }
    }
}

get_global_idは、引数で指定した次元のワークアイテムのIDを取得します。
get_global_sizeは、引数で指定した次元のワークアイテムID数を取得します。
get_global_idの呼び出し結果が各ワークアイテムで異なるので、その値を利用することで、ワークアイテムごとに異なる処理を行わせることができます。

上記コードでは、ワークアイテムIDを使用し、それぞれのワークアイテムが別々の行を計算するようにしています。
担当の行の計算が終わったら、他のワークアイテムが担当している行はスキップし、次の担当の行を処理するようになっています。(i += groupSizeの部分)

計算結果

cpu: 35 msecs, gpu: 28 msecs

なんとか速くなりましたがまだ遅い……。

GPUを走らせるために、コマンドをエンキューしたりイベントを待ったりするので、ある程度オーバーヘッドがあります。
今の行列のサイズ程度ではオーバーヘッドの方が大きく出ているのかもしれません。
試しに、行数を10倍にして様子をみてみましょう。

    enum {
        ROWS = 1000,
        COLS = 200,
        RESULT_COLS = 300
    }

cpu: 372 msecs, gpu: 85 msecs

大きなサイズのデータであれば、CPUと比較してだいたい4倍速い感じが見えてきました。

列方向にも分割する

ここまでのコードでは、行方向にしか並列計算をさせていません。
極端な話、1行しかない巨大ベクトルが引数に来ると1ワークアイテムしか動作せず、超遅くなります。
そこで、列方向にもワークアイテムを分割するようにしてみます。

__kernel void product(
        __global const float *lhs,
        __global const float *rhs,
        __global float *result,
        uint rows,
        uint cols,
        uint resultCols) {
    // ワークアイテムのグローバルIDを取得する。
    const size_t groupI = get_global_id(0);
    const size_t groupRows = get_global_size(0);
    const size_t groupJ = get_global_id(1);
    const size_t groupCols = get_global_size(1);

    for(size_t i = groupI; i < rows; i += groupRows) {
        for(size_t j = groupJ; j < resultCols; j += groupCols) {
            float value = 0.0f;
            for(size_t k = 0; k < cols; ++k) {
                value += lhs[i * rows + k] * rhs[k * resultCols + j];
            }
            result[i * resultCols + j] = value;
        }
    }
}

さて、こちらの列分割版コードのためにワークアイテムも2次元にする必要があります。
ワークグループ数13程度、1ワークグループあたり32ワークアイテム程度という制約を保ったまま2次元に展開する必要があるので、ちょっと面倒です。

ワークグループ数の13という数(素数!)は掛け算で作れないので、妥協してワークグループ数を16とします。
16であれば、縦横4つずつワークグループができれば良いことになります。
また、推奨ワークグループサイズの32も平方根が整数にならないので、いっそのこと32 * 32を1ワークグループとします。最大ワークグループサイズが1024なので収まります。

cl.enqueueKernel(
    commandQueue,
    kernel,
    [4 * 32, 4 * 32], // グローバルワークアイテム数
    [32, 32]  // ローカルワークアイテム数
);

実行結果

cpu: 38 msecs, gpu: 3 msecs

おお〜なぜかわからないけど速くなりました。ようやくCPU比10倍に……。
(おそらく、列方向に分割したことでメモリフットプリントが小さくなったり、ワークアイテム数が増えたことでメモリアクセス中にも計算が捗るようになったのだと思います)

調子に乗って行列のサイズをさらに増やしてみます。

    enum {
        ROWS = 1000,
        COLS = 2000,
        RESULT_COLS = 3000
    }

cpu: 48344 msecs, gpu: 2100 msecs

CPU比でほぼ20倍ですね。

次回予告

なんとかGPUの優越を示せたものの、まだまだ明らかに改善できる部分があります。
次回以降は、今回のコードをベースにして、行列計算をどんどん最適化していきます。