CUDAの内部：パフォーマンスエンジニアリング

動機と振り返り

前回の記事 — Hello CUDA! — では、CUDAの概要とGPUアーキテクチャの全体像を探りました。その続きとして、この記事では、GPUのパワー（と複雑さ）を構成する中心的なトピックであるパフォーマンスエンジニアリングについて議論します。

GPUは、その優れた並列コンピューティング能力で長年知られています。数十億のパラメータを持つディープラーニングモデルのトレーニングから、分子動力学のシミュレーション、金融リスク分析に至るまで、すべてGPUを必要とします。

CUDAでは、最適化されていないコードと最適化されたコードのパフォーマンス差は、数十倍、あるいは数百倍になることもあります。このスピードアップを達成するには、単に「GPUコードを書く」だけでは不十分です。私たちは並列で考えることを学ばなければなりません。これには、並列アルゴリズムと、それがハードウェアとどのように直接相互作用するかについての深い理解が必要です。¹ 

この記事では、GPUコアをビジー状態に保ち、ハードウェアの全能力を引き出す方法を学びます。並列アルゴリズムを深く掘り下げる代わりに、各カーネルの実際のパフォーマンスを決定する基本的なテクニックに焦点を当てます。

ハードウェアと実験セットアップ#

パフォーマンスは常にハードウェアによって制限されます。ある最適化手法が特定のアーキテクチャで効果的に機能しても、別のものでは非効率的かもしれません。 

客観的に測定するため、2つの異なるハードウェアプラットフォームを使用します：

• NVIDIA GeForce MX330: 古く、低スペックのラップトップGPU。ローエンドハードウェア（Pascalアーキテクチャ）を代表します。
• NVIDIA RTX A4000: Ampereアーキテクチャに基づく、モダンで強力なワークステーションGPU。

両方で測定することで、どのテクニックが一般的な利点（両方で効果的）を提供し、どのテクニックが特定のアーキテクチャでのみ効果的なのかがわかります。

以下は、これら2つのGPUの主要な仕様の比較表です：

簡単な分析： 差は歴然です。A4000 (Ampere) は、SM数が16倍（48対3）であるだけでなく、メモリ帯域幅も約8倍（448対56.1 GB/s）です。これは、メモリバウンドなカーネルがA4000で大幅に恩恵を受け、コンピュートバウンドなカーネルが膨大な数のSM/CUDAコアから恩恵を受けることを示唆しています。

パフォーマンス最適化テクニック#

CUDAで最も重要な最適化テクニックのいくつかを、最も強力なツールから順に見ていきましょう。

1. 共有メモリ (SMEM)#

グローバルメモリ（GPUのVRAM）へのアクセスは、CUDAカーネルにおいて最も高価な操作の一つです。レイテンシが非常に高く、数百クロックサイクルに達する可能性があります。カーネルが常にグローバルメモリとの読み書きを行っていると、パフォーマンスは深刻な「ボトルネック」に見舞われます。

共有メモリ (SMEM) は、オンチップ（SMの内部に配置）のメモリ領域で、サイズは小さい（通常SMあたり48KB～128KB）ですが、アクセス速度と帯域幅が非常に高いです。SMEMのレイテンシはL2キャッシュよりも低く、レジスタよりわずかに高いだけで、ブロック内でのデータ再利用に効果的なツールとなります。¹

一般的な戦略は次のとおりです：

• データのブロックをグローバルメモリから共有メモリにロードします（一度だけ、そしてコアレッシングされた方法でロードしようと試みます - セクション3参照）。
• ブロック内のスレッドは、共有メモリ内のデータにアクセスして多くの計算を実行します。
• 最終結果をグローバルメモリに書き戻します。

例1：行列積

$C = A \times B$（サイズ $M \times N$ と $N \times K$）の2つの行列を乗算する問題を考えます。

ナイーブなカーネル（SMEMなし）

このカーネルは単純です。各スレッドがC行列の1つの要素を計算します。しかし、各スレッドは行列AからN回、行列BからN回の読み取りを実行する必要があり、これらはすべてグローバルメモリからのものです。

__global__ void matrix_multiplication_naive(
    const float *A, const float *B, float *C, 
    int M, int N, int K)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (row < M && col < K)
    {
        float sum = 0.0f;
        // N回ループし、毎回グローバルメモリにアクセス
        for (int k = 0; k < N; ++k)
        {
            sum += A[row * N + k] * B[k * K + col];
        }
        C[row * K + col] = sum;
    }
}

最適化されたカーネル（タイリングとSMEM）

「タイリングと呼ばれるテクニックを使用します。各ブロックは、C行列の1つの「タイル」、例えば $16 \times 16$ の計算を担当します。そのために、AとBから対応する $16 \times 16$ のタイルを共有メモリにループでロードします。

ブロック内の各スレッドの実行イメージ：

• ループ（タイル数）
• ロード： タイルAの1要素とタイルBの1要素を（グローバルメモリから）SMEMにロードします。
• 同期（sync）： ブロック内の全スレッドがロードを完了するのを待ちます（__syncthreads();）。
• 計算： （現在SMEMにある）サブタイルの内積を計算します。
• 同期（sync）： 次のタイルをロードする前に、全スレッドが計算を終えるのを待ちます（__syncthreads();）。
• 書き込み： 最終的な合計をグローバルメモリのC行列に書き込みます（最後に一度だけ）。

#define TILE_SIZE 16
 
__global__ void matrix_multiplication_smem(
    const float *A, const float *B, float *C, 
    int M, int N, int K)
{
    // AとBのタイルのための共有メモリを宣言
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];
 
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
 
    int tx = threadIdx.x; // タイル内のx座標
    int ty = threadIdx.y; // タイル内のy座標
 
    float sum = 0.0f;
 
    // 必要なすべてのタイルにわたってループ
    for (int t = 0; t < (N + TILE_SIZE - 1) / TILE_SIZE; ++t)
    {
        // 1. グローバルメモリから共有メモリへタイルをロード
        // グローバルメモリ内の 'A' の列インデックスと 'B' の行インデックスを計算
        int A_col = t * TILE_SIZE + tx;
        int B_row = t * TILE_SIZE + ty;
 
        // 各スレッドがAの1要素をAsにロード
        if (row < M && A_col < N) {
            As[ty][tx] = A[row * N + A_col]; // ty, tx は SMEM 内の座標
        } else {
            As[ty][tx] = 0.0f; // 範囲外の場合はパディング
        }
 
        // 各スレッドがBの1要素をBsにロード
        if (B_row < N && col < K) {
            Bs[ty][tx] = B[B_row * K + col];
        } else {
            Bs[ty][tx] = 0.0f; // 範囲外の場合はパディング
        }
 
        // 2. 全スレッドのロード完了を保証するために同期
        __syncthreads();
 
        // 3. 共有メモリから計算
        #pragma unroll // コンパイラにこのループを展開するようヒントを与える
        for (int i = 0; i < TILE_SIZE; ++i)
        {
            sum += As[ty][i] * Bs[i][tx];
        }
 
        // 4. 次のタイルに移動する前に同期
        __syncthreads();
    }
 
    // 5. 最終結果をグローバルメモリに書き込む
    if (row < M && col < K)
    {
        C[row * K + col] = sum;
    }
}

結果（例1：行列積）

分析： この測定データは、非常に重要な教訓を教えてくれます：最適化はアーキテクチャに強く依存するということです。

• MX330 (Pascal) の場合： 予測通り、タイリングとSMEMのテクニックは明確な利益をもたらし、ナイーブなカーネルと比較して約2.35倍高速化（37.240 / 15.846）しました。古いハードウェアでは、L2キャッシュが小さくメモリ帯域幅も狭いため、SMEMを使用してグローバルメモリアクセスを減らすことは、パフォーマンスを大幅に向上させます。

• RTX A4000 (Ampere) の場合： 驚くべき結果が出ました。私たちのタイリングとSMEMカーネル（2.077 ms）は、実際にはナイーブなカーネル（1.861 ms）よりも約1.1倍遅くなりました。

  • これはSMEMが役に立たないという意味ではありません。ナイーブなカーネルが、その単純さにもかかわらず、Ampereアーキテクチャのハードウェアメカニズムから大きな恩恵を受けていることを示しています。大きなL2キャッシュ（4MB）と非常に高いメモリ帯域幅（448 GB/s）により、A4000のハードウェアは、私たちが手動で介入しなくても、グローバルメモリアクセスのレイテンシのほとんどを自動的に隠蔽してくれた可能性があります。
  • 一方、タイリングとSMEMカーネルは、より複雑なインデックス計算ロジック、そして特に __syncthreads() 命令といった新しい「オーバーヘッド」コストを導入しました。Ampereのような強力なアーキテクチャでは、これらの同期バリアのコストが、SMEMを使用することで得られる利益よりも大きくなる可能性があります。特に、ナイーブなカーネルがすでに非常に高速に実行されている場合には。

• cuBLAS ライブラリ： どちらの場合も、cuBLAS が圧倒的です。特にA4000では、ナイーブなカーネルより約9倍、私たちのSMEMカーネルより約10倍高速です。これは、cuBLAS がアセンブリレベルで最適化されているだけでなく、Ampereアーキテクチャで利用可能なTensor Cores（テンソルコア）を最大限に活用するように設計されているためです。これは、私たちの基本的なfloatカーネルが使用していない、行列積演算に特化したハードウェアです。

例2：アトミック操作の最適化

atomicAdd() は、複数のスレッドがグローバルメモリ内の同じ変数を更新する際のレースコンディションを避けるために必要な関数です。¹

問題：配列要素のカウント

入力配列 input（サイズ N）が与えられたとき、配列内で特定の値 K に等しい要素の総数をカウントする必要があります。この総数は、グローバルメモリ内の単一の変数 output に蓄積され、格納されなければなりません。

ナイーブなカーネル（ナイーブなアトミック）

最も単純な解決策は、有効な要素を見つけた各スレッドが、グローバルメモリの output 変数に対して直接 atomicAdd() を呼び出すことです。

__global__ void count_equal_naive(
    const int *input, int *output, int N, int K)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N && input[idx] == K)
    {
        // 多数のアトミック操作が「ボトルネック」を引き起こす
        atomicAdd(output, 1);
    }
}

問題： もし input に K と等しい要素が多数含まれていると、数千のスレッドが同時に同じアドレスに対して atomicAdd を呼び出そうと「競合」します。これらのアクセスはシリアライズされ、並列性が完全に失われます。

最適化されたカーネル（SMEMと並列リダクション）

より良い戦略は、グローバルメモリへの atomicAdd 呼び出し回数を最小限にすることです。ブロックごとに1回だけ呼び出すようにします。

• グリッドストライドループ： 各スレッドは複数の要素を（strideだけ離れた）処理します。これは任意のサイズの配列を処理するための強力なパターンです。
• 各スレッドは、見つけた要素をローカルでカウントします。
• このローカルカウントの結果は、共有メモリ（I[]）内の配列に格納されます。
• 並列リダクション： ブロックのローカルカウントをすべて合計するために、共有メモリ内で直接並列リダクションサムを実行します。この「リダクション」はSMEMとブロック内のスレッドを使用し、グローバルメモリへの atomicAdd のような競合を完全に回避します。これは非常に効率的です。
• 最後に、ブロック内の1つのスレッドだけ（例：threadIdx.x == 0）が、ブロックの合計を output に加算するために atomicAdd を一度だけ呼び出します。

__global__ void count_equal_optimized(
    const int *input, int *output, int N, int K)
{
    // カーネル起動時にSMEMを動的に割り当てるために 'extern' を使用
    extern __shared__ int I[];
    
    int g_idx = blockIdx.x * blockDim.x + threadIdx.x;
    int l_idx = threadIdx.x;
    // ストライドはグリッド内の総スレッド数
    int stride = gridDim.x * blockDim.x;
 
    // 1. 各スレッドがローカルでカウント（グリッドストライドループ）
    int count = 0;
    while (g_idx < N)
    {
        if (input[g_idx] == K)
        {
            count++;
        }
        g_idx += stride; // 次の要素へジャンプ
    }
    
    // 2. ローカル結果をSMEMに保存
    I[l_idx] = count;
    __syncthreads();
 
    // 3. SMEM内で並列リダクションを実行
    for (int s = blockDim.x / 2; s > 0; s >>= 1)
    {
        if (l_idx < s)
        {
            I[l_idx] += I[l_idx + s];
        }
        __syncthreads();
    }
 
    // 4. スレッド0のみがブロックの合計をグローバルメモリに書き込む
    if (l_idx == 0)
    {
        atomicAdd(output, I[0]);
    }
}

結果（例2：アトミック操作の最適化）

分析： この結果は、ハードウェアアーキテクチャの重要性についての観察をさらに裏付けます。

• MX330 (Pascal) の場合： 最適化テクニック（グローバルメモリアトミックの負荷軽減）は、約1.34倍の利益（12.271 / 9.1230）をもたらしました。これは妥当な結果です。古いPascalアーキテクチャでは、グローバルメモリでの atomicAdd の競合は高価なボトルネックだからです。SMEMでのリダクションに切り替えることで、この競合が大幅に減少しました。

• RTX A4000 (Ampere) の場合： パフォーマンスの向上はほとんどなく、約1%の改善（1.000 / 0.990）にとどまりました。

  • Ampere (Compute Capability 8.6) のような最新のGPUアーキテクチャは、Pascalよりもはるかに効率的なアトミック処理システムを持っています。グローバルメモリの atomicAdd 操作は、L2キャッシュで効率的に処理・集約（coalesce）されることが多く、競合コストが大幅に削減されます。
  • 私たちの「最適化された」カーネルは、Ampereのハードウェアがすでにほぼ完全に解決していたボトルネック（グローバルアトミック）を、SMEMへの書き込み、複数の __syncthreads() 命令、複雑なロジックといった他の一連の操作に置き換えただけでした。
  • その結果、SMEMでの並列リダクションの総オーバーヘッド時間は、グローバルアトミックの負荷を軽減することで節約された時間とほぼ等しくなりました。言い換えれば、私たちは最新のハードウェアではもはや主要な問題ではない問題を「最適化」したのです。

一般的な結論： 古典的な最適化テクニック（SMEMタイリング、アトミックの削減など）は依然として非常に価値がありますが、その有効性はもはや絶対的なものではありません。あるテクニックが、あるハードウェア世代では「救世主」であっても、キャッシュ、帯域幅、専用処理ユニット（L2キャッシュアトミック、Tensor Coresなど）の改善により、次の世代では非効率（あるいは逆効果）になる可能性があります。

2. メモリコアレッシング#

これは、グローバルメモリアクセスを最適化する上で最も重要な概念の一つです。

ワープ（同期して実行される32スレッドのグループ）内の32スレッドがグローバルメモリにアクセスするとき、GPUハードウェアはこれらの32のリクエストをできるだけ少ないメモリトランザクションに**「コアレッシング（合体）」**しようとします。

• 理想的なケース（コアレッシング）： 32スレッドがメモリ内の32の連続した場所（例：A[idx], A[idx+1], ..., A[idx+31]）にアクセスします。GPUハードウェアは、メモリを32バイトまたは128バイトの「セグメント」単位で読み取ります。アクセスがコアレッシングされると、32スレッドのワープが32個のfloat値（32 * 4 = 128バイト）を読み取るのに、単一の128バイトトランザクション1回で済みます。
• 最悪のケース（非コアレッシング）： 32スレッドが32のランダムまたは不連続な場所（ストライドアクセス）にアクセスします。例えば、32スレッドが A[idx * 100] にアクセスすると、それらは32の異なるメモリセグメントに分散する可能性があります。GPUは32回の個別の128バイトトランザクションを実行する必要があり、帯域幅の31/32を無駄にします！

ナイーブな行列積の例では、A[row * N + k]（行単位）へのアクセスは通常コアレッシングされます（ワープ内の隣接スレッドは隣接したcol値を持つため）。逆に、B[k * K + col]（列単位）へのアクセスはしばしば非コアレッシングとなり、深刻な帯域幅の無駄を引き起こします。¹

3. バンクコンフリクト#

この問題は共有メモリで発生します。SMEMは一枚岩のブロックではなく、32のメモリ「バンク」に分割されています。ワープ内のスレッドは、異なるバンクにアクセスする限り、並列でSMEMにアクセスできます。

バンクコンフリクト： ワープ内の2つ以上のスレッドが、同じバンクに位置するアドレスにアクセスしようとすると発生します。これらのアクセスはシリアライズされ、SMEMの速度が失われます。

ルール： 具体的には、4バイトワードサイズ（floatやint用）のSMEMの場合、バンク $k$ は $(addr / 4) \pmod{32} = k$ となるアドレス $addr$ の4バイトワードを含みます。

典型的な例：

__shared__ float A[32][32];

• 行単位のアクセス： A[my_row][threadIdx.x]
  • 32スレッドが A[r][0], A[r][1], ..., A[r][31] にアクセス。
  • これらは連続したアドレスであり、32の異なるバンク（バンク0, 1, 2, ...）に分散します。
  • 結果：非常に高速、コンフリクトなし。
• 列単位のアクセス： A[threadIdx.x][my_col]
  • スレッド0が A[0][c]（例：バンク1）にアクセス
  • スレッド1が A[1][c]（A[0][c]から32要素離れている）にアクセス。そのアドレスも mod 32 の結果、バンク1になります。
  • スレッド2が A[2][c]（A[0][c]から64要素離れている）にアクセス。そのアドレスも mod 32 の結果、バンク1になります。
  • 結果：大惨事！ 32ウェイのバンクコンフリクト！ 32スレッドがバンク1にアクセスし、シリアライズされます。

32のバンクを32のレジ（checkout lanes）と想像してください：

• 行単位アクセス： 32人（スレッド）が32の異なるレジ（バンク0...31）に進みます。全員が並列に処理されます。
• 列単位アクセス： 32人全員（スレッド）が1つのレジ（例：バンク1）に並びます。31人が待たなければなりません。

回避方法： パディング。配列の次元を変更してアクセスパターンを崩します：

// 1列のパディングを追加
__shared__ float A[32][33]; 

これで、列単位で A[threadIdx.x][my_col] にアクセスすると：

• スレッド0が A[0][c]（バンク c % 32）にアクセス
• スレッド1が A[1][c]（33要素離れている）にアクセス。バンク (c + 33) % 32 = バンク (c+1) % 32
• スレッド2が A[2][c]（66要素離れている）にアクセス。バンク (c + 66) % 32 = バンク (c+2) % 32
• 結果：コンフリクト解消！ 並列アクセス速度と引き換えに、少量のSMEMを「犠牲」にしました。

4. オキュパンシー（占有率）#

オキュパンシー（Occupancyとは、SM上のアクティブなワープ数を、そのSMがサポートできる最大ワープ数で割った比率です（例：32アクティブワープ / 64最大ワープ = 50%オキュパンシー）。

オキュパンシーはレイテンシ隠蔽の核となる要素です。あるワープが停止（ストール）しなければならない場合（例：グローバルメモリからのデータを待つ）、SMのスケジューラは即座に準備完了状態にある別のワープの実行に切り替えることができます。

• 低オキュパンシー： SMに切り替えるための準備完了ワープが十分にありません。唯一のアクティブワープがストールすると、SM全体が「アイドル」状態になり、計算リソースが無駄になります。
• 高オキュパンシー： SMには選択肢となるワープが多数あり、メモリレイテンシを隠蔽し、コンピュートコアをビジー状態に保つのに役立ちます。

これが根本的な違いです：CPUはレイテンシ指向（1つのタスクを非常に速く完了させようとする）であるのに対し、GPUはスループット指向（単位時間あたりに最も多くのタスクを完了させようとする）です。高オキュパンシーは、このスループットモデルの鍵です。

オキュパンシーは、SM上のリソースのうち、最初に枯渇したものによって制限されます：

• レジスタ： 各スレッドが使用するレジスタが多すぎると、SMは多くのスレッドを保持するのに十分なレジスタを持てなくなり、オキュパンシーが低下します。
• 共有メモリ： 各ブロックが使用するSMEMが多すぎると、SMは多くのブロックを保持するのに十分なSMEMを持てなくなり、オキュパンシーが低下します。
• ブロックあたりのスレッド数： ブロックあたりのスレッド数を少なすぎる値（例：64）に設定すると、ワープは32スレッドであるため、決して高いオキュパンシーを達成できません。
• ハードウェアの制限： 各SMには、管理できるブロックとスレッドの最大数があります。例えば、A4000 (CC 8.6) は、最大 1536 スレッド / SM（48ワープに相当）および 32 ブロック / SM を処理できます。¹

5. カーネルフュージョン（カーネル融合）#

これは、（連続して実行される）複数のカーネルを単一のカーネルに統合するテクニックです。

例1： 以下を実行する代わりに：

• kernel_add(A, B, C); // C = A + B
• kernel_scale(C, D); // D = C * alpha

これらを次のように統合します：

• kernel_add_and_scale(A, B, D); // D = (A + B) * alpha

例2： 「SAXPY」操作（Y = a*X + Y）または Y = a*X + b。

• 最適化なし：
  • kernel_scale(X, a, Temp); // Temp = a*X (Temp をグローバルメモリに書き込む)
  • kernel_add(Temp, b, Y); // Y = Temp + b (Temp をグローバルメモリから読み込む)
• 最適化（フュージョン）：
  • kernel_fused(X, a, b, Y); // 単一のカーネルが Y[i] = a*X[i] + b[i] を計算。Temp はレジスタ内にのみ存在し、Temp のグローバルメモリへの読み書きを完全に排除します。

利点：

• カーネル起動オーバーヘッドの削減： すべての __global__ 呼び出しには小さなオーバーヘッドがあります。
• グローバルメモリアクセスの排除： これが最大の利点です。上記の例では、中間変数（C または Temp）をグローバルメモリに書き込み、次のカーネルがそれをグローバルメモリから読み戻す必要があります。融合されたカーネルは、この中間値をレジスタに一時的に保持でき、2回の高価なグローバルメモリアクセスを完全に排除します。

トレードオフ：

• 複雑さの増大： 融合されたカーネルは、記述やデバッグが難しくなります。
• リソース圧迫の増大： 新しいカーネルはより多くの作業を行うため、（中間変数を保持するために）より多くのレジスタを必要とする可能性があり、オキュパンシーを低下させる可能性があります。

カーネルフュージョンはトレードオフです。メモリ帯域幅の要件を減らす代わりに、SM上の計算リソース要件を増やす可能性があります。そのトレードオフが見合うかどうかは、測定して確認する必要があります。

結論#

GPUの最適化は、理論的な理解と実践的な経験の両方を必要とするプロセスです。それは、行列積や配列の総和といった、一見おなじみの問題を、まったく新しい並列の視点から見直すことを私たちに強います。

このプロセスの核心は、カーネルのパフォーマンスボトルネックを特定することです。それがメモリバウンドなのか、コンピュートバウンドなのか、あるいはレイテンシバウンドなのか。NVIDIA Nsight のようなツールを使用してカーネルを「プロファイリング」し、どこに労力を集中すべきかを正確に見つけることが不可欠です。

私たちは、基本的でありながら非常に重要なテクニックを見てきました。将来の記事では、ダイナミック・パラレリズム、ディープラーニングのためのTensor Coresのパワー（Tensor CoresはVoltaアーキテクチャ（CC ≥ 7.0）以降で利用可能）、そしてストリームを使った非同期操作など、より高度なトピックを探求するかもしれません。