これまでのあらすじ

ローカルメモリーの導入で計算速度を数倍に上げることはできましたが、まだまだ理論性能からは遠いです。

今回は頑張ってさらなる各種最適化を行い、なんとか理論性能に近づけられるよう頑張ります。

ソースコード

下準備

これからどんどん試行錯誤で最適化を行うので、以下の下準備をしました。

構成改善・表示改善

CPUとの比較テストに時間がかかって面倒なので、構成(dubのconfig)で有無を切り替えられるようにしました。

$ dub run --config=examples-product #CPUテストなし
$ dub run --config=examples-product-with-cpu-test #CPUテストあり

dub.sdl

configuration "examples-product" {
    targetType "executable"
    mainSourceFile "examples/product.d"
    stringImportPaths "examples/"
}

configuration "examples-product-with-cpu-test" {
    targetType "executable"
    mainSourceFile "examples/product.d"
    stringImportPaths "examples/"
    versions "DcltkWithCpuTest"
}

また、ベンチマークを本来の計算部分のみに絞りました。
結果表示についても、GFLOPSを出したりCPU比を出すように改善しました。

examples/product.d


    // 純粋な計算部分のみ計測するよう修正
    void productGpu() {
        cl_event event;
        cl.enqueueKernel(commandQueue, kernel, [32 * 4, 32 * 4], [32, 32], event);
        cl.waitAndReleaseEvents(event);
    }

    // FLOPSも算出する
    immutable gpuMsecs = benchmark!(() => productGpu())(4)[0].total!"msecs" / 4;
    immutable gpuFlops = (cast(real) ROWS) * RESULT_COLS * (COLS * 2.0) / ((cast(real) gpuMsecs) / 1000.0);
    
    // version指定によりCPUテスト有無を切り替える
    version(DcltkWithCpuTest) {
        cl_event event;
        cl.enqueueReadBuffer(commandQueue, resultBuffer, 0, gpuResult, event);
        cl.waitAndReleaseEvents(event);

        immutable cpuMsecs = benchmark!(() => productCpu(
                    lhs, rhs, cpuResult, ROWS, COLS, RESULT_COLS))(1)[0].total!"msecs";

        // CPUテスト時は速度のCPU比も出す
        writefln("cpu: %d msecs, gpu: %d msecs (%.1f GFLOPS, faster %.1f times)",
            cpuMsecs, gpuMsecs, gpuFlops / (10.0^^9), cast(real) cpuMsecs / cast(real) gpuMsecs);

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

        // check result values.
        foreach(i, e; cpuResult) {
            assert(approxEqual(e, gpuResult[i]));
        }
    } else {
        // GPUのみ
        writefln("gpu: %d msecs (%.1f GFLOPS)", gpuMsecs, gpuFlops / (10.0^^9));
    }

OpenCL Cのコード部分は別ファイル(product.cl)に切り出し、import("product.cl")で取り込むようにしました。

product.d

    auto program = cl.createProgramFromSource(context, import("product.cl"));
    scope(exit) cl.releaseProgram(program);
    cl.buildProgram(program, deviceIds);

product.cl

__kernel void product(
    __global const float *lhs,
    __global const float *rhs,
    __global float *result,
    uint rows,
    uint cols,
    uint resultCols,
    __local float *localRow,
    __local float *localCol) {
    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);

    const size_t localI = get_local_id(0);
    const size_t localRows = get_local_size(0);
    const size_t localJ = get_local_id(1);
    const size_t localCols = get_local_size(1);

    for(size_t i = 0; i < rows; i += groupRows) {
        for(size_t j = 0; j < resultCols; j += groupCols) {
            float value = 0.0f;

            for(size_t k = 0; k < cols; k += localCols) {

                barrier(CLK_LOCAL_MEM_FENCE);
                if((i + groupI) < rows && (k + localJ) < cols) {
                    localRow[localI * localCols + localJ] = lhs[(i + groupI) * cols + (k + localJ)];
                }
                if((j + groupJ) < resultCols && (k + localI) < cols) {
                    localCol[localI * localCols + localJ] = rhs[(k + localI) * resultCols + (j + groupJ)];
                }
                barrier(CLK_LOCAL_MEM_FENCE);

                if((i + groupI) < rows && (j + groupJ) < resultCols) {
                    for(size_t lk = 0; lk < localCols && (k + lk) < cols; ++lk) {
                        value += localRow[localI * localCols + lk] * localCol[lk * localCols + localJ];
                    }
                }
            }
            if((i + groupI) < rows && (j + groupJ) < resultCols) {
            	result[(i + groupI) * resultCols + (j + groupJ)] = value;
            }
        }
    }
}

CPUでの計算の高速化

CPUでの計算があまりにも遅すぎるので、std.paralellismを使って並列化しました。
これでCPUながらだいぶ速くなります。（複数のコアがあれば）

examples/product.d

import std.parallelism : parallel;
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 {
    foreach(i; parallel(iota(0, rows))) { // 変更箇所。行をCPUコアに振り分ける
        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;
        }
    }
}

修正前

cpu: 104929 msecs, gpu: 1367 msecs (8.8 GFLOPS, faster 76.8 times)

修正後

cpu: 24676 msecs, gpu: 1366 msecs (8.8 GFLOPS, faster 18.1 times)

こんなに手軽に5倍くらい速くなるのは嬉しいですね。

GPUをちゃんと動かす

さて、前回までなんとなくコーディングで話を進めてきましたが、そもそもGPUはちゃんと動いているのでしょうか……。

私はテスト用にGoogle Compute EngineのGPU(Tesla K80)を付けたインスタンスを使っているのですが、実はドキュメントにこんな記載があります。

GPU パフォーマンスの最適化

まだやってなかった。

実行前

cpu: 44077 msecs, gpu: 798 msecs

実行

$ sudo nvidia-smi -pm 1
Enabled persistence mode for GPU 00000000:00:04.0.
All done.
$ sudo nvidia-smi -ac 2505,875
Applications clocks set to "(MEM 2505, SM 875)" for GPU 00000000:00:04.0
All done.
$ sudo nvidia-smi --auto-boost-default=DISABLED
All done.

実行後

cpu: 50330 msecs, gpu: 648 msecs

なんと、これだけで20%も速度が改善されました。

実際にGPU借りて計算していると、これだけで20%も費用が削減できてしまいます。
ちゃんとドキュメントは読まないとダメですね。

行列のパディングと`if`の除去

歴史とGPGPUにifは禁物です。
GPGPUは同一処理を一斉に行うことで速度を上げていますが、ifが絡むと処理が全然同一でなくなって困ります。
なんでも、アーキテクチャによっては、if内を実行しないスレッドのみメモリを変更しないフラグを立てて処理を空回しして済ましているそうです。

現状のコードには、ワークアイテムの境界判定用のifがたくさん入っています。これを削減するにはどうするべきでしょうか。
境界判定が必要なのは、行列のサイズがワークグループのサイズと一致していないことが原因です。行列の端からはみ出るワークアイテムがあって、それがなんの疑問もなく動くと、メモリアクセス違反や意図しない計算の実行に繋がります。

今回、そこをどうにかして行列サイズをワークグループのサイズの倍数に合わせ、if文を排除します。

行列をパディング

さて、行列をワークグループのサイズに合わせるとなると、行数・列数は下記の通りになります。

現状
- 1000行・2000列 * 2000行・3000列 = 1000行・3000列
ワークグループサイズ(32*32)に合わせた場合
- 1024行・2016列 * 2016行・3008列 = 1024行・3008列

このワークグループサイズの倍数に合わせた行列サイズであれば、ワークアイテムが処理の中で余ることがなく、ifでの境界判定を行う必要が無くなります。
行列サイズによる計算時間の増加とトレードオフですが、そもそも無駄にifでスキップしている分が埋まるだけなので、たぶんきっとどうにかなると思います。

ワークグループサイズに合わせて増えた部分は、0を埋めておきます。

行列サイズを増やして計算結果は大丈夫か?

行列サイズを増やした場合、以下の手順を踏んで計算を行うことになります。

元の行列をワークグループサイズ倍の行列にコピー(範囲外は0埋め)
計算を行う。
ワークグループサイズ倍の結果行列から本来の結果の部分をコピー(範囲外は捨てる)

ここで、**行列のサイズを増やすと計算結果が変わってしまうのでは?**という疑問が生じます。(私は生じました)

しかし、以下の3パターンに分けて考えると大丈夫そうだという結論に至りました。

左辺の行列の行が増える場合 -> 結果の行列に0埋めの行が増えるのみ -> 結果の行列の切り取りで問題なし
右辺の行列の列が増える場合 -> 結果の行列に0埋めの列が増えるのみ -> 結果の行列の切り取りで問題なし
左辺の行列の列・右辺の行列の行が増える場合 -> 積和演算で0*0の項が増えるのみ -> 結果の行列の切り取りで問題なし

たぶん線形代数とかでちゃんと証明できると思うのですが、学がないので保留とさせてください……。

末尾に行または列を足すだけであれば問題ないので、パディングした行列でもちゃんと計算が行えそうです。

実装

以下の処理を行うよう実装していきます。

ワークグループサイズに合わせた行列のバッファーを確保
0埋めする
両辺の行列データをホストから矩形コピー
計算(あまり変わらない)
計算結果の行列をホストへ矩形コピー

0埋めにはclEnqueueFillBufferを、矩形コピーにはclEnqueueWriteBufferRectとclEnqueueReadBufferRect をそれぞれ使用しました。

product.d

    enum {
        // 行列サイズ
        ROWS = 1000,
        COLS = 2000,
        RESULT_COLS = 3000,

        // ワークグループの処理単位
        BATCH_ROWS = 32,
        BATCH_COLS = 32
    }

// 中略

    // 行列サイズをワークグループの処理単位(BATCH_COLS,BATCH_ROWS)に切り上げる。
    immutable bufferCols = cast(uint) roundUp(COLS, BATCH_COLS);
    immutable bufferRows = cast(uint) roundUp(ROWS, BATCH_ROWS);
    immutable bufferResultCols = cast(uint) roundUp(RESULT_COLS, BATCH_COLS);

    // 行列の要素数を計算
    immutable lhsSize = bufferCols * bufferRows;
    immutable rhsSize = bufferResultCols * bufferCols;
    immutable resultSize = bufferResultCols * bufferRows;

    // バッファの作成
    auto lhsBuffer = cl.createReadBuffer(context, lhsSize * float.sizeof);
    scope(exit) cl.releaseBuffer(lhsBuffer);
    auto rhsBuffer = cl.createReadBuffer(context, rhsSize * float.sizeof);
    scope(exit) cl.releaseBuffer(rhsBuffer);
    auto resultBuffer = cl.createWriteBuffer(context, resultSize * float.sizeof);
    scope(exit) cl.releaseBuffer(resultBuffer);

    // 0埋め
    cl.enqueueFillBuffer(
        commandQueue, lhsBuffer, [0.0f], 0, lhsSize);
    // ホストからバッファへの矩形コピー。なんか色々面倒だったのをラップして多少楽にした。
    cl.enqueueWriteBuffer(
        commandQueue,
        lhsBuffer,
        cl.Position(0, 0),
        cl.Region(COLS, ROWS),
        bufferCols,
        lhs);
    // 同上
    cl.enqueueFillBuffer(
        commandQueue, rhsBuffer, [0.0f], 0, rhsSize);
    cl.enqueueWriteBuffer(
        commandQueue,
        rhsBuffer,
        cl.Position(0, 0),
        cl.Region(RESULT_COLS, COLS),
        bufferResultCols,
        rhs);

// 中略

    // グローバルワークアイテム数を引数の行列サイズと一致させる。
    immutable(size_t)[] globalWorkSizes = [
        bufferResultCols,
        bufferRows
    ];

    // ローカルワークアイテム数は定数にした通り。
    immutable(size_t)[] localWorkSizes = [BATCH_COLS, BATCH_ROWS];

    // 計算実行
    void productGpu() {
        cl_event event;
        cl.enqueueKernel(commandQueue, kernel, globalWorkSizes, localWorkSizes, event);
        cl.waitAndReleaseEvents(event);
    }

カーネル関数はだいぶすっきりします。

product.cl

__kernel void product(
        __global const float *lhs,
        __global const float *rhs,
        __global float *result,
        uint rows,
        uint cols,
        uint resultCols,
        __local float *localLhs,
        __local float *localRhs) {
    const size_t localJ = get_local_id(0);
    const size_t localCols = get_local_size(0);
    const size_t localI = get_local_id(1);
    const size_t localRows = get_local_size(1);
    const size_t globalJ = get_global_id(0);
    const size_t globalI = get_global_id(1);

    float value = 0.0f;

    for(size_t k = 0; k < cols; k += localCols) {
        barrier(CLK_LOCAL_MEM_FENCE);
        localLhs[localI * localCols + localJ] = lhs[globalI * cols + (k + localJ)];
        localRhs[localI * localCols + localJ] = rhs[(k + localI) * resultCols + globalJ];
        barrier(CLK_LOCAL_MEM_FENCE);

        for(size_t lk = 0; lk < localCols; ++lk) {
            value += localLhs[localI * localCols + lk] * localRhs[lk * localCols + localJ];
        }
    }
    result[globalI * resultCols + globalJ] = value;
}

各ワークアイテムが自分の担当要素だけ気にすればよくなったので、外側のforや境界判定のifが無くなりました。これは期待できるかも……。

結果

$ cpu: 51624 msecs, gpu: 51 msecs (235.3 GFLOPS, faster 1012.2 times)

10倍速くなりました!
ようやくGPUらしい速度になってきましたね。

次回予告

次回は、さらにプライベートメモリー等を使用して速くしていきます。
できればTFLOPSの単位まで速くしてテラフロ野郎（漢字で書くと寺風呂野郎、英語で言えばテラフロッパー）を目指したいところです。