Skip to content
Published on

GPUカーネルを自分の手で直すということ — 転置カーネルを5倍速くするまで

🇰🇷KO🇺🇸EN🇯🇵JA🇨🇳ZH
Split View필사 모드
シェア
Authors

はじめに — プロファイラが示したカーネルがライブラリの中に見つからないとき

学習ステップひとつが420msで、そのうち170msが名前も知らないカーネルひとつに入っている状況に出くわすと、選択肢が一気に狭まります。cuBLAS呼び出しなら手の出しようがなく、PyTorchの演算子なら別の演算子に置き換えてみればよいでしょう。しかしそのカーネルが、自分のモデル特有の風変わりなマスキングロジックのために誰かが急いで書いたものだとしたら、直せるのは自分しかいません。

この記事はその地点から始まります。GPUカーネルを「自分の手で直す」とは実際どのような作業なのか、何を知っておくべきで何を測るべきかを扱います。結論を先に言えば、カーネル最適化の8割は演算を減らすことではなく、メモリを動かす順序を変えることです。この記事の実習で直すカーネルも、算術演算がひとつも入っていないカーネルです。それでも5倍以上速くなります。

基準環境はCUDA Toolkit 13.3 Update 1、Nsight Compute 2026.2系です。概念は世代を問いませんが、ツールのフラグ名は変わるので、コマンドラインは自分のバージョンのNsight Compute CLIドキュメントと照らし合わせてください。

カーネルとは何か、そしてなぜ直すことになるのか

カーネルはGPU上で実行される関数ひとつです。CPUの関数と違う点は、一度の呼び出しで数万個のインスタンスが同時に開始されることです。このインスタンスひとつひとつがスレッドであり、スレッドは自分が何番目かを組み込み変数から把握し、それに対応するデータの断片だけを処理します。

カーネルを直接直す価値がある状況は、実務では思ったより狭いものです。順に点検するとこうなります。

状況先にすべきことカーネルを直す理由になるか
標準的なGEMM、畳み込みが遅いcuBLAS、cuDNNのバージョンとデータ型、テンソルコア経路を確認ほぼならない
小さな演算子が何個も連なって実行されるtorch.compileでフュージョンを誘導大抵ならない
アテンションの変形が必要FlashAttention系ライブラリにその変形があるか確認なければなる
カーネルひとつが理論帯域幅の20パーセントしか使わないNsight Computeで原因を確認なる
自分のドメイン特有のインデキシング、マスキング、疎パターン代替できるライブラリがないなる

肝心なのは最後の二行です。ライブラリで代替できるなら、代替する方が常に有利です。ライブラリのカーネルには自分が使える時間よりはるかに多くのチューニングが入っており、新しいアーキテクチャが出れば自分の代わりに更新されます。自分で書いたカーネルは、自分が一生面倒を見なければなりません。

実行モデル — スレッド、ワープ、ブロック、グリッド

カーネルを直すには、ハードウェアがスレッドをどうまとめるかを知る必要があります。階層は四段階です。

グリッド(grid)         カーネル呼び出し一回分の全スレッド集合
  └ ブロック(block)    同じSMに配置され、シェアードメモリと __syncthreads() を共有
      └ ワープ(warp)   32スレッド。スケジューリングと命令発行の実際の単位
          └ スレッド    自分のレジスタを持つ実行フローひとつ

ここで実務上もっとも重要な事実は、ワープこそが本当の単位であるということです。プログラマはスレッド単位でコードを書きますが、ハードウェアは32個をひとかたまりにして発行します。ここから二つの帰結が生まれます。

第一に、分岐発散です。ひとつのワープ内のスレッドが互いに異なる分岐を取ると、ハードウェアは両方の経路を順番に実行し、該当しないスレッドを無効化します。ワープの中で分かれる条件文は実行時間を足します。ワープの境界に合わせて分かれる条件文は無料です。

第二に、メモリアクセスもワープ単位でまとめられます。これがコアレッシングであり、次節のテーマです。

ブロックサイズを決めるときは32の倍数にするのが基本です。ブロックサイズが33だと、ハードウェアはワープを二つ発行し、そのうち31レーンが遊ぶことになります。

メモリ階層 — コアレッシングとバンクコンフリクト

GPUメモリは層になっており、層ごとにレイテンシと帯域幅が桁違いです。

階層範囲おおまかな性質
レジスタスレッドひとつ最速。スレッドあたりの個数が占有率を左右する
シェアードメモリブロックひとつSM内部のSRAM。プログラマが直接管理するキャッシュ
L1 / テクスチャキャッシュSMひとつシェアードメモリと物理的に同じ記憶装置を分け合う
L2キャッシュGPU全体すべてのSMが共有。HBM前の最後の防衛線
HBM(グローバルメモリ)GPU全体容量は大きいがレイテンシが数百サイクル

コアレッシング

グローバルメモリアクセスは32バイト単位のトランザクションとして処理されます。ひとつのワープの32スレッドが連続したfloat32個、つまり128バイトを読めば、トランザクション4回で終わります。逆に同じ32スレッドが4096間隔で離れたfloat32個を読むと、トランザクションが32回必要になり、そのたびに32バイトを取ってきて4バイトだけ使い、残りを捨てます。帯域幅の8分の1しか使っていないことになります。

これがカーネル性能差のもっとも大きな単一要因です。実習で正確にこの8倍の無駄を目にすることになります。

バンクコンフリクト

シェアードメモリは32個のバンクにインターリーブされています。4バイトワード基準で、アドレスを4で割った商を32で割った余りがバンク番号です。ひとつのワープのスレッドが互いに異なるバンクに触れれば一サイクルで処理され、同じバンクの異なるアドレスに触れればその数だけ直列化されます。

典型的な事故は正方形のシェアード配列への列方向アクセスです。tile[32][32]でtile[i][0]をiについて走査すると、全要素がバンク0に落ちます。32-wayコンフリクトです。解法は配列をひと列広げることです。tile[32][33]と宣言すると、行ごとの開始バンクが一つずつずれ、列アクセスが32個のバンクに均等に分散します。シェアードメモリを32 * 4バイト、つまりブロックあたり128バイト余分に使い、直列化をなくす取引です。

占有率とルーフライン — 何を見て何を無視するか

占有率は目標ではなく症状である

占有率は、SMが常駐させられる最大ワープ数に対する実際の常駐ワープ数です。初心者にありがちな誤解は、これを最大化すべき目標とみなすことです。そうではありません。

占有率がやることはひとつだけです。レイテンシを隠すことです。あるワープがHBMの応答を待っている間に別のワープを実行すればSMは遊びません。つまり占有率は「レイテンシを隠すだけのワープが十分にあるか」という問いの代理指標であって、それ自体が性能ではありません。

低い占有率の方がむしろ速い場合が二つあります。

第一に、スレッドあたりのレジスタを多く使って命令レベル並列性を確保している場合です。スレッドひとつが独立したロードを4つ同時に飛ばしておけば、ワープ数が4分の1でも進行中のメモリ要求の総量は変わりません。Vasily VolkovのBetter Performance at Lower Occupancyがこの論旨を正面から扱っています。2010年の資料ですが、論旨は今も有効です。

第二に、カーネルがすでに帯域幅の上限に張り付いている場合です。メモリパイプが飽和状態なら、ワープを増やしても入る余地がありません。

そこで実務のルールはこうなります。占有率は低いときだけ見ます。25パーセントを下回っていてカーネルがレイテンシに縛られているなら、レジスタ使用量やシェアードメモリ割り当てを疑います。60パーセント出ているのにカーネルが遅いなら、占有率は犯人ではないので他を見るべきです。

ルーフライン — 大半のカーネルは帯域幅に縛られている

どこを見るかを決める道具がルーフラインです。軸は算術強度、すなわち動かしたバイトひとつあたりの演算回数です。

達成性能(FLOP/s)
   ^
   |            ______________  演算上限
   |           /
   |          /   傾き = メモリ帯域幅
   |         /
   +--------+-------------------> 算術強度 (FLOP/Byte)
          転換点

転換点 = (演算上限 FLOP/s) / (メモリ帯域幅 Byte/s)

転換点はハードウェアの特性です。最新のデータセンターGPUではこの値は数十から数百FLOP/Byteの範囲です。ところが実際に使うカーネルの算術強度を計算してみると、大半は一桁です。

演算おおまかな算術強度位置
要素ごとの加算1 FLOP / 12 Byte極端なメモリバウンド
活性化関数数FLOP / 8 Byteメモリバウンド
行列転置0 FLOP / 8 Byte純粋なメモリ
LayerNorm一桁台FLOP / Byteメモリバウンド
GEMM(大きな行列)タイルサイズに比例して数百まで演算バウンド
LLMのデコード段階バッチが小さいと2未満メモリバウンド

読み方は単純です。算術強度が転換点よりはるかに小さければ、そのカーネルの性能上限はすでに決まっています。演算命令をいくら減らしても無駄で、唯一意味のある改善は動かすバイトを減らすか、動かし方を直すことです。

だからカーネル最適化の最初の問いは常にこれです。「このカーネルの理論最小トラフィックは何バイトで、今実際に何バイト動かしているか」。

実習 — 転置カーネルを四段階で直す

では実際に直してみます。対象は4096 x 4096のfloat行列転置です。演算が0回なのでメモリの話だけが残り、教材として理想的です。

理論最小トラフィックは明確です。一回読んで一回書くので2 * 4096 * 4096 * 4バイト、約134MBです。これより少なく動かす方法はありません。したがって性能の上限は「転置せずにただコピーするだけのカーネル」であり、これを先に測って基準線とします。

全体コード

// transpose.cu
// ビルド: nvcc -O3 -arch=sm_80 transpose.cu -o transpose
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

static const int TILE = 32;
static const int BLOCK_ROWS = 8;   // ブロックあたり 32x8 = 256 スレッド
static const int N = 4096;

#define CHECK(x) do { cudaError_t e_ = (x); if (e_ != cudaSuccess) { \
    printf("CUDA error: %s (line %d)\n", cudaGetErrorString(e_), __LINE__); \
    exit(1); } } while (0)

// 0段階。上限線: 転置せずコピーだけする。
__global__ void copyKernel(float *out, const float *in) {
  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = in[(y + j) * N + x];
}

// 1段階。naive: 読みはコアレッシングされるが、書きがN間隔のストライドになる。
__global__ void transposeNaive(float *out, const float *in) {
  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[x * N + (y + j)] = in[(y + j) * N + x];
}

// 2段階。シェアードメモリタイル: 転置をSRAM内で終わらせ、
//        グローバルの読みと書きの両方をコアレッシングさせる。
__global__ void transposeShared(float *out, const float *in) {
  __shared__ float tile[TILE][TILE];

  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * N + x];

  __syncthreads();

  // ブロック座標を入れ替えて、書きも連続アドレスになるようにする。
  x = blockIdx.y * TILE + threadIdx.x;
  y = blockIdx.x * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = tile[threadIdx.x][threadIdx.y + j];
}

// 3段階。パディング一列でシェアードメモリのバンクコンフリクトをなくす。
__global__ void transposePadded(float *out, const float *in) {
  __shared__ float tile[TILE][TILE + 1];   // 唯一の違い

  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * N + x];

  __syncthreads();

  x = blockIdx.y * TILE + threadIdx.x;
  y = blockIdx.x * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = tile[threadIdx.x][threadIdx.y + j];
}

typedef void (*Kern)(float *, const float *);

static void bench(const char *name, Kern k, float *d_out, const float *d_in,
                  const float *h_ref, float *h_out, bool checkTranspose) {
  dim3 grid(N / TILE, N / TILE), block(TILE, BLOCK_ROWS);
  const int WARMUP = 5, ITERS = 50;
  const double bytes = 2.0 * N * N * sizeof(float);

  for (int i = 0; i < WARMUP; i++) k<<<grid, block>>>(d_out, d_in);
  CHECK(cudaDeviceSynchronize());

  cudaEvent_t t0, t1;
  CHECK(cudaEventCreate(&t0));
  CHECK(cudaEventCreate(&t1));
  CHECK(cudaEventRecord(t0));
  for (int i = 0; i < ITERS; i++) k<<<grid, block>>>(d_out, d_in);
  CHECK(cudaEventRecord(t1));
  CHECK(cudaEventSynchronize(t1));

  float ms = 0.f;
  CHECK(cudaEventElapsedTime(&ms, t0, t1));
  double perIter = ms / ITERS;
  double gbs = bytes / (perIter * 1.0e-3) / 1.0e9;

  // 正確性検証のない性能数値には意味がない。
  CHECK(cudaMemcpy(h_out, d_out, (size_t)N * N * sizeof(float),
                   cudaMemcpyDeviceToHost));
  long bad = 0;
  for (long r = 0; r < N && bad == 0; r++)
    for (long c = 0; c < N; c++) {
      float want = checkTranspose ? h_ref[c * N + r] : h_ref[r * N + c];
      if (h_out[r * N + c] != want) { bad++; break; }
    }

  printf("%-18s %8.3f ms   %8.1f GB/s   %s\n", name, perIter, gbs,
         bad ? "FAIL" : "ok");
  CHECK(cudaEventDestroy(t0));
  CHECK(cudaEventDestroy(t1));
}

int main() {
  size_t bytes = (size_t)N * N * sizeof(float);
  float *h_in = (float *)malloc(bytes), *h_out = (float *)malloc(bytes);
  for (long i = 0; i < (long)N * N; i++) h_in[i] = (float)(i % 1000);

  float *d_in, *d_out;
  CHECK(cudaMalloc(&d_in, bytes));
  CHECK(cudaMalloc(&d_out, bytes));
  CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

  cudaDeviceProp p;
  CHECK(cudaGetDeviceProperties(&p, 0));
  printf("%s  peak HBM = %.1f GB/s\n\n", p.name,
         2.0 * p.memoryClockRate * (p.memoryBusWidth / 8) / 1.0e6);

  bench("copy (upper bound)", copyKernel,      d_out, d_in, h_in, h_out, false);
  bench("naive",              transposeNaive,  d_out, d_in, h_in, h_out, true);
  bench("shared tile",        transposeShared, d_out, d_in, h_in, h_out, true);
  bench("shared + padding",   transposePadded, d_out, d_in, h_in, h_out, true);

  cudaFree(d_in); cudaFree(d_out); free(h_in); free(h_out);
  return 0;
}

測定方法で重要なこと

数値を信じるには、まずハーネスが正直でなければなりません。上のコードが守っている規則は五つです。

  • ウォームアップを捨てます。 最初の呼び出しにはコンテキスト生成とモジュールロードが混ざり込みます。
  • 繰り返し測定して平均を取ります。 カーネルひとつが1msくらいだと、クロック変動だけで10パーセント揺れます。
  • CPUタイマーではなくcudaEventを使います。 カーネル実行は非同期なので、CPU時間は発行時間しか測れません。
  • 正確性を検証します。 誤ったインデキシングは大抵の場合速くなります。検証のないGB/sは数字遊びです。
  • 時間ではなく有効帯域幅に換算します。 絶対時間はサイズと機材によって変わりますが、理論帯域幅に対して何パーセントかはどこでも比較可能です。

有効帯域幅の式は簡単です。動かすべき最小バイトを実際にかかった時間で割ります。ここで「動かすべき」という点が重要です。無駄にして実際に行き来したバイトではなく、アルゴリズム上必要なバイトを使います。そうしてこそ無駄が数値に現れます。

結果の形

以下はA100 80GB(sm_80)級の機材でこのハーネスを回したときに出る典型的な形です。絶対値は機材、ドライバ、クロック状態によって大きく変わるので、そのまま引用せず、上のコードを自分のGPUで実際に動かして自分の基準線を作ってください。 意味があるのは段階間の相対比率です。

段階有効帯域幅コピー比ボトルネック
copy(上限線)基準値100100パーセントなし。HBM飽和
naive約18約18パーセント書きが非コアレッシング。トランザクションあたり4バイトしか使わない
shared tile約63約63パーセントシェアードメモリの32-wayバンクコンフリクト
shared + padding約93約93パーセント実質なし。タイル境界効果だけが残る

読み取るべきことは三つです。

第一に、naiveが上限の5分の1水準だという点です。演算は0回で動かすデータ量も同じなのに5倍遅い。違いは順序だけです。書きがN間隔で散らばっているため、32バイトのトランザクションごとに4バイトだけ書いて28バイトを捨てます。8倍の無駄が他の効果と混ざって5倍の差として現れます。

第二に、シェアードタイルが大きな幅を回復するものの最後まで行かない点です。グローバルアクセスは両方直したのに、ボトルネックがSRAMの中に移りました。tile[threadIdx.x][threadIdx.y + j]は列方向アクセスであり、32x32の正方配列では列アクセスが全部同じバンクに落ちます。

第三に、最後の段階のコードの違いが配列宣言一文字だという点です。[TILE]を[TILE + 1]に変えたのがすべてです。カーネル最適化はしばしばこういう形をしています。アルゴリズムではなく、データ配置をひとマスずらすことです。

よくある失敗パターン

この実習で実際によく踏む地雷です。

  • __syncthreads()の欠落。 シェアードタイルを埋めて読む前に同期がないと、結果が非決定的に間違います。小さな入力では偶然合うことが多く、より危険です。
  • __syncthreads()を分岐の中に入れること。 ブロック内の一部のスレッドしか到達しない位置に置くと未定義動作です。
  • インデックス計算でブロック座標を入れ替えないこと。 シェアードタイルだけ入れて出力インデックスをそのままにすると、書きが再びストライドになり2段階目の利得が消えます。結果は合っているのに速くならない、もっとも気づきにくい形です。
  • -O3なしで測ること。 ホストコードの最適化が抜けると、検証ループが測定時間を支配して結論が逆転します。
  • Nが小さすぎる。 カーネル実行オーバーヘッドが数マイクロ秒なので、総時間が数十マイクロ秒だとオーバーヘッドを測ることになります。

Nsight Computeで何を見るか

数値が悪かったとき原因を教えてくれるのがプロファイラです。コマンドラインから見ていきます。

# 全セクション収集。カーネルひとつに数百msかかるので対象を絞る必要がある。
ncu --set full \
    --kernel-name regex:transpose \
    --launch-skip 5 --launch-count 1 \
    -o transpose_report \
    ./transpose

# 原因特定用の指標だけをピンポイントで抜く(はるかに速い)
ncu --metrics \
  sm__throughput.avg.pct_of_peak_sustained_elapsed,\
gpu__dram_throughput.avg.pct_of_peak_sustained_elapsed,\
l1tex__data_bank_conflicts_pipe_lsu_mem_shared.sum,\
l1tex__average_t_sectors_per_request_pipe_lsu_mem_global_op_ld.ratio \
  --kernel-name regex:transpose --launch-count 1 ./transpose

# GUIで開く
ncu-ui transpose_report.ncu-rep

--launch-skipでウォームアップ実行を飛ばすことが重要です。プロファイラは最初の実行のキャッシュコールド状態をそのまま見せるので、飛ばさないと定常状態でないものを分析することになります。

レポートを開くとセクションがいくつもありますが、実務では見る順序が決まっています。

1. Speed of Light。 演算処理量とメモリ処理量をハードウェアの理論値に対するパーセントで見せます。ここで方向が決まります。メモリが80パーセント以上なら帯域幅に張り付いており、両方とも30パーセント未満ならレイテンシか占有率の問題です。私たちの転置カーネルのnaive版はこの画面で両方とも低く出ます。パイプが飽和しているのではなく、無駄にしているからです。

2. Memory Workload Analysis。 ここが核心です。リクエストひとつあたり何個のセクタを取ってきたかを見せてくれます。完璧にコアレッシングされた32スレッドのfloatロードはリクエストあたり4セクタです。naiveカーネルの書きはリクエストあたり32セクタ出ます。8倍の無駄がこの一行にそのまま刻まれます。この指標ひとつが「コアレッシングの問題か」に即答します。

3. Shared Memory関連指標。 バンクコンフリクトの回数が出ます。2段階目のカーネルはここが大きく跳ね上がり、3段階目は0に近く落ちます。パディングが実際に効いたかを確認する場所です。

4. Warp State Statistics。 ワープが止まっていた理由を種類別に見せます。Stall Long Scoreboardが圧倒的ならグローバルメモリの応答待ちであり、Stall MIO Throttleならシェアードメモリや特殊関数ユニット側の混雑です。原因がメモリか命令かを切り分けてくれます。

5. Occupancy。 最後に見ます。前の四つがきれいなのに依然として遅いとき、そこで常駐ワープが不足していないか確認します。この順序を守らずに占有率から先に見ると、大抵見当違いの方向に最適化してしまいます。

Nsight Computeにはルーフラインセクションもあり、自分のカーネルが坂の上にいるのか平地にいるのかを図で見せてくれます。転置カーネルは算術強度が0なのでいちばん左端に打たれ、それだけで「演算の最適化はやることがない」という結論が出ます。

この作業をいつやめるか

転置カーネルを上限の93パーセントまで上げたなら、残り7パーセントを追う理由はほとんどありません。判断基準をあらかじめ決めておく方がよいでしょう。

  • 理論最小トラフィック比90パーセントを超えたらやめます。 メモリバウンドなカーネルではそこから上はタイル境界とTLB効果なので、努力に対する回収が急激に悪化します。
  • 全体実行時間に占めるこのカーネルの比重を測り直します。 170msだったものを40msに縮めたなら、今度は別の場所がボトルネックです。アムダールの法則はカーネル最適化でもそのまま働きます。
  • 維持コストを計算します。 自分で書いたカーネルは新しいアーキテクチャが出るたびに再検証の対象です。ライブラリより20パーセント速いカーネルが2年後には30パーセント遅くなっているかもしれません。
  • ひとつ上の層で解決するかを先に見ます。 次の記事で扱うTritonで同じカーネルを20行で書いて似た性能が出るなら、CUDA C++版を維持する理由が減ります。

おわりに — カーネルを直すとは、データが動く順序を直すことである

この記事の実習には浮動小数点演算が一度も登場しませんでした。それでも最初のバージョンと最後のバージョンの間に5倍の差が出ました。変わったのは、同じデータをどんな順序で読み、どこに一時的に置いておき、どんな順序で書くか、それだけです。

この事実がGPUカーネル作業の性格を規定します。演算量を減らすアルゴリズム的改善は、大抵ライブラリがすでにやっているか、自分たちの問題では変えられません。私たちに残された梃子はメモリ階層内でのデータ配置と移動順序であり、幸いそちらの方がはるかに大きな梃子です。

作業順序を一文にまとめるとこうなります。理論最小トラフィックを計算し、ハーネスで今何パーセント使っているかを測り、Nsight Computeで無駄の位置を特定し、配置を直し、また測る。勘で直して速くなったと言うのは、この順序から測定を二回抜いたものであり、そうやって得た結論は次の機材で覆ります。

参考資料

Discuss on Twitter • View on GitHub