NVIDIA GPUとCUDAの完全解剖：なぜGPUがAIを支配するのか

CPUアーキテクチャ（レイテンシ最適化）:
┌────────────────────────────────────────────┐
│ Core 0  │ Core 1  │ Core 2  │ Core 3   │
│(高性能) │(高性能) │(高性能) │(高性能)  │
│ OOO実行 │ OOO実行 │ OOO実行 │ OOO実行  │
├─────────┴─────────┴─────────┴────────────┤
│           L3キャッシュ（32MB以上）          │
├────────────────────────────────────────────┤
│       メインメモリ（DDR5, 約50 GB/s）       │
└────────────────────────────────────────────┘

GPUアーキテクチャ（スループット最適化）— H100:
┌──┬──┬──┬──┬──┬──┬──┬──┬──┬──┬──┐
│SM│SM│SM│SM│SM│SM│SM│SM│SM│SM│...│  132個のSM
│  │  │  │  │  │  │  │  │  │  │   │  各SM = 128個のCUDAコア
└──┴──┴──┴──┴──┴──┴──┴──┴──┴──┴───┘
       L2キャッシュ（50MB）
       HBM3（80GB、3.35 TB/s）

GPU Device（H100 SXM5）
│
├── SM 0（Streaming Multiprocessor）
│   ├── 4個のWarpスケジューラ
│   ├── 4個のDispatch Unit
│   ├── 128個のCUDAコア（FP32）
│   ├── 64個のFP64コア
│   ├── 4個のTensor Core（第4世代：FP8/FP16/BF16/TF32対応）
│   ├── 1個の特殊関数ユニット（SFU：sin、cos、sqrtなど）
│   ├── LD/STユニット（Load/Store）
│   ├── L1キャッシュ + 共有メモリ = 228KB（比率設定可能）
│   └── レジスタファイル（65536個 × 32ビットレジスタ）
│
├── SM 1 ...
├── SM 2 ...
│   ...
├── SM 131 ...
│
├── L2キャッシュ（50MB）
└── HBM3メモリ（80GB、3.35 TB/s帯域幅）

Thread Block（例：256スレッド）
├── Warp 0: スレッド  0-31  → 同じ命令をロックステップで実行
├── Warp 1: スレッド 32-63  → 同じ命令をロックステップで実行
├── Warp 2: スレッド 64-95  → 同じ命令をロックステップで実行
├── ...
└── Warp 7: スレッド 224-255

// このコードはWarp Divergenceを引き起こす
__global__ void divergent_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx % 2 == 0) {
        // Warp内の偶数インデックスのスレッドだけがここを実行
        data[idx] = data[idx] * 2.0f;  // ステップ1
    } else {
        // Warp内の奇数インデックスのスレッドだけがここを実行
        data[idx] = data[idx] + 1.0f;  // ステップ2
    }
    // 結果：ステップ1実行中、奇数スレッドはアイドル（マスクオフ）
    //       ステップ2実行中、偶数スレッドはアイドル（マスクオフ）
    // 実効スループット：理論値の50%
}

// Divergenceなし：Warp内すべてのスレッドが同じパスを取る
__global__ void coherent_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // ブロック内のすべてのスレッドは同じblockIdx.xを共有するため
    // Warp内32スレッドは全員同じ分岐を取る
    if (blockIdx.x % 2 == 0) {
        data[idx] = data[idx] * 2.0f;
    } else {
        data[idx] = data[idx] + 1.0f;
    }
}

メモリ階層（高速 → 低速）:

┌─────────────────────────────────────────────────────┐
│  レジスタ（32ビット、SM当たり65536個）               │
│  レイテンシ：約1サイクル  帯域幅：膨大（ローカルアクセス）│
│  スコープ：単一スレッドのみ                          │
├─────────────────────────────────────────────────────┤
│  L1キャッシュ / 共有メモリ（H100ではSM当たり228KB）  │
│  レイテンシ：約32サイクル  帯域幅：約19 TB/s         │
│  スコープ：Thread Block内のすべてのスレッド          │
├─────────────────────────────────────────────────────┤
│  L2キャッシュ（H100では50MB）                       │
│  レイテンシ：約200サイクル  帯域幅：約12 TB/s        │
│  スコープ：全SMが共有                               │
├─────────────────────────────────────────────────────┤
│  グローバルメモリ HBM3（H100では80GB）              │
│  レイテンシ：約700サイクル  帯域幅：3.35 TB/s        │
│  スコープ：デバイス上のすべてのスレッド              │
└─────────────────────────────────────────────────────┘

// Naive行列乗算カーネル — なぜ遅いかを理解するためのバージョン
__global__ void matmul_naive(float* A, float* B, float* C, int N) {
    // 各スレッドが出力行列の1要素を担当
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            // 各イテレーションでグローバルメモリに2回アクセス → 約700サイクルのレイテンシ
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

// 起動設定の例
dim3 blockDim(16, 16);  // 256スレッド/ブロック
dim3 gridDim((N + 15) / 16, (N + 15) / 16);
matmul_naive<<<gridDim, blockDim>>>(A, B, C, N);

#define TILE_SIZE 16

// 共有メモリを活用したタイリング行列乗算
__global__ void matmul_tiled(float* A, float* B, float* C, int N) {
    // 共有メモリにタイルを宣言（HBMの約700サイクルに対し約32サイクル）
    __shared__ float tile_A[TILE_SIZE][TILE_SIZE];
    __shared__ float tile_B[TILE_SIZE][TILE_SIZE];

    int tx = threadIdx.x, ty = threadIdx.y;
    int row = blockIdx.y * TILE_SIZE + ty;
    int col = blockIdx.x * TILE_SIZE + tx;
    float sum = 0.0f;

    // K次元をタイルサイズ単位でイテレート
    for (int t = 0; t < N / TILE_SIZE; t++) {
        // ステップ1：各スレッドが1要素ずつ共有メモリにロード
        tile_A[ty][tx] = A[row * N + t * TILE_SIZE + tx];
        tile_B[ty][tx] = B[(t * TILE_SIZE + ty) * N + col];

        // ステップ2：すべてのスレッドのロード完了を待つ
        __syncthreads();

        // ステップ3：タイル内の内積を計算
        // すべてのアクセスが共有メモリ（約32サイクル）にヒット
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tile_A[ty][k] * tile_B[k][tx];
        }

        // ステップ4：次のタイルをロードする前に同期
        __syncthreads();
    }

    C[row * N + col] = sum;
}

タイリングの原理：

行列A:           行列B:
┌──┬──┬──┬──┐  ┌──┬──┬──┬──┐
│T0│T1│T2│T3│  │T0│T1│T2│T3│
├──┼──┼──┼──┤  ├──┼──┼──┼──┤
│  │  │  │  │  │  │  │  │  │
└──┴──┴──┴──┘  └──┴──┴──┴──┘

各Thread BlockはCの16×16出力タイルを1つ担当。
Aの行タイルとBの列タイルを順番に共有メモリにロードし、
計算後に次のタイルへ進む → グローバルメモリ転送量をTILE_SIZE倍削減。

// WMMAを使ったTensor Core直接アクセス
#include <mma.h>
using namespace nvcuda::wmma;

__global__ void tensor_core_matmul(half* a_ptr, half* b_ptr,
                                    float* c_ptr, int M, int N, int K) {
    // 各WarpがCの16×16出力タイル1つを担当
    // fragment = Tensor Coreが操作する行列の断片
    fragment<matrix_a, 16, 16, 16, half, row_major> a_frag;
    fragment<matrix_b, 16, 16, 16, half, col_major> b_frag;
    fragment<accumulator, 16, 16, 16, float> c_frag;

    // アキュムレータをゼロ初期化
    fill_fragment(c_frag, 0.0f);

    // K次元をイテレート
    for (int k = 0; k < K; k += 16) {
        // グローバルメモリからfragmentをロード
        load_matrix_sync(a_frag, a_ptr + /* offset */, K);
        load_matrix_sync(b_frag, b_ptr + /* offset */, N);

        // 単一命令で16×16×16の行列乗算累算
        // 内部でTensor Coreハードウェアを直接呼び出す
        mma_sync(c_frag, a_frag, b_frag, c_frag);
    }

    // 結果をグローバルメモリに書き戻す
    store_matrix_sync(c_ptr + /* offset */, c_frag, N, mem_row_major);
}

CUDA実行階層：

Grid（カーネル全体の実行）
├── Block (0,0): 256スレッド
│   ├── Warp 0: スレッド  0-31 ── SMで同時実行
│   ├── Warp 1: スレッド 32-63
│   └── ...
├── Block (1,0): 256スレッド
├── Block (2,0): 256スレッド
└── ...

各BlockはちょうどひとつのSMに割り当てられる。
ひとつのSMは複数のBlockを同時にホストできる（Occupancy）。

import torch

A = torch.randn(4096, 4096, device='cuda', dtype=torch.float16)
B = torch.randn(4096, 4096, device='cuda', dtype=torch.float16)

# この1行の内部で起きていること：
# 1. cuBLASまたはCUTLASSが最適なカーネルを選択
# 2. 実行構成: 例) Grid(256,256)、Block(16,16)
# 3. Tensor Coreを活用したTiled GEMM実行
# 4. 約4096個のThread Block、各256スレッド
result = torch.matmul(A, B)  # 内部でcublasSgemmExを呼び出す

SMリソース制限（H100）：
- 同時Thread Block最大数：32
- 同時スレッド最大数：2048
- レジスタ：65536個（常駐Warpの全スレッドで分割）
- 共有メモリ：228KB（常駐Block間で分割）

Occupancy = アクティブWarp数 / 最大可能Warp数
目標：50%以上（Warpスイッチングでメモリレイテンシを隠蔽するため）

import torch
import time

# 1. Tensor Coreを有効化するためFP16またはBF16を使用
A = torch.randn(8192, 8192, device='cuda', dtype=torch.bfloat16)
B = torch.randn(8192, 8192, device='cuda', dtype=torch.bfloat16)

# 2. ウォームアップ：最初の実行にはカーネルのコンパイル/選択に時間がかかる
for _ in range(5):
    _ = torch.matmul(A, B)
torch.cuda.synchronize()

# 3. 計測
start = time.perf_counter()
for _ in range(100):
    C = torch.matmul(A, B)
torch.cuda.synchronize()
end = time.perf_counter()

# 4. TFLOPSを計算
# 8192x8192の行列乗算：2 * N^3 = 2 * 8192^3 = 1.1e12 FLOP
flops = 2 * 8192**3 * 100
elapsed = end - start
tflops = flops / elapsed / 1e12
print(f"達成性能：{tflops:.1f} TFLOPS")
# 目標：700+ TFLOPS（理論値989の70%以上が達成可能）