30 天 GPU LeetCode 挑戰：從 CUDA 新手到 Kernel Leaderboard

平常使用 PyTorch 訓練模型時，只要把 Tensor 搬到 CUDA，就可以使用 GPU 加速。但 .cuda() 背後到底發生什麼事情？為什麼 GPU 可以比 CPU 更適合矩陣運算？當一個 CUDA Kernel 很慢時，又要從哪裡開始優化？

因此我決定用 30 天時間，從完全理解 CUDA 的基本執行模型開始，一路練習 GPU Kernel，最後挑戰 Kernel Leaderboard。

這系列不只是要讓程式「可以跑」，而是希望慢慢理解如何讓 Kernel 跑得更快。

一段原本跑在 CPU 上的程式，要怎麼搬到 GPU，並真正利用 GPU 的平行運算能力？

這篇會從一個最簡單的陣列相加開始，逐步理解：

• CPU 與 GPU 的差異
• CUDA Kernel 是什麼
• <<<grid, block>>> 到底代表什麼
• threadIdx.x、blockIdx.x、blockDim.x、gridDim.x
• 為什麼單純增加 thread 會產生 race condition
• index 與 stride 的概念
• Grid-Stride Loop
• Unified Memory 與 Page Fault
• Prefetching 為什麼可以讓 kernel 大幅變快
• 為什麼 Vector Add 最後會變成 Memory-Bandwidth Bound

1. 從 CPU 版本開始

先看最普通的 C++ 陣列相加：

#include <iostream>
#include <math.h>

// function to add the elements of two arrays
void add(int n, float *x, float *y)
{
    for (int i = 0; i < n; i++)
        y[i] = x[i] + y[i];
}

int main(void)
{
    int N = 1 << 20;  // 1,048,576 elements

    float *x = new float[N];
    float *y = new float[N];

    for (int i = 0; i < N; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    add(N, x, y);

    float maxError = 0.0f;
    for (int i = 0; i < N; i++)
        maxError = fmax(maxError, fabs(y[i] - 3.0f));

    std::cout << "Max error: " << maxError << std::endl;

    delete[] x;
    delete[] y;

    return 0;
}

這段程式做的事情很單純：

x = [1, 1, 1, 1, ...]
y = [2, 2, 2, 2, ...]

y[i] = x[i] + y[i]

最後：

y = [3, 3, 3, 3, ...]

其中：

int N = 1 << 20;

代表：

2^20 = 1,048,576

也就是大約一百萬個元素。

如果使用 Mac，可以直接用：

clang++ add.cpp -o add
./add

執行這個 CPU 版本。

2. CPU 與 GPU 思維的差異

CPU 版本的核心是：

for (int i = 0; i < n; i++)
    y[i] = x[i] + y[i];

概念上是一個 worker 依序處理：

i = 0
i = 1
i = 2
...

GPU 的思維則是：

Thread 0 → y[0]
Thread 1 → y[1]
Thread 2 → y[2]
...

也就是：

把同樣的工作拆成很多份，交給大量 threads 同時執行。

3. 第一個 CUDA Kernel

把原本的 add() 變成 CUDA kernel：

__global__
void add(int n, float *x, float *y)
{
    for (int i = 0; i < n; i++)
        y[i] = x[i] + y[i];
}

__global__ 代表：

這個函式由 CPU 啟動，但實際在 GPU 上執行。

接著使用：

add<<<1, 1>>>(N, x, y);

啟動 kernel。

CUDA 的 kernel launch syntax：

kernel<<<numBlocks, threadsPerBlock>>>(...);

所以：

add<<<1, 1>>>(N, x, y);

代表：

1 個 Block
×
1 個 Thread

這時其實還沒有真正平行化，只是把原本 CPU function 搬到 GPU。

4. 為什麼 <<<1, 256>>> 不能直接套上去？

如果直接改成：

add<<<1, 256>>>(N, x, y);

但 kernel 還是：

for (int i = 0; i < n; i++)
    y[i] = x[i] + y[i];

那會變成：

Thread 0 → 處理整個 array
Thread 1 → 處理整個 array
Thread 2 → 處理整個 array
...
Thread 255 → 處理整個 array

也就是每一個 thread 都會讀寫相同的 y[i]。

這會產生 Race Condition。

例如：

Thread 0: 讀 y[0] = 2
Thread 1: 讀 y[0] = 2

Thread 0: 寫 y[0] = 3
Thread 1: 可能又讀到 3，最後寫成 4

結果會依執行順序而不同。

5. cudaDeviceSynchronize() 為什麼不能解決 Race Condition？

cudaDeviceSynchronize() 的作用是：

讓 CPU 等待 GPU 完成工作。

例如：

add<<<...>>>(...);
cudaDeviceSynchronize();

代表：

CPU launch kernel
↓
GPU 開始執行
↓
CPU 等待
↓
GPU 完成
↓
CPU 繼續

它解決的是：

CPU ↔ GPU

之間的同步。

而 race condition 是：

GPU Thread 0 ↔ GPU Thread 1

同時存取相同 memory 的問題，所以 cudaDeviceSynchronize() 無法解決。

6. Thread、Block、Warp

CUDA 的結構可以先理解成：

Grid
  ↓
Block
  ↓
Thread

GPU 實際排程時，thread 又會以 Warp 為單位執行。

1 Warp = 32 Threads

因此 block size 通常會選 32 的倍數，例如：

128
256
512

如果 block size 是 256：

256 / 32 = 8 warps

也就是一個 block 共有 8 個 warps。

為什麼常看到 256？

因為 256 是一個很常見的起始值：

• 是 32 的倍數
• 有足夠多的 threads
• 不會讓單一 block 太大
• 在很多 workload 上通常有不錯的 occupancy

但要注意：

256 並不是永遠最佳，只是一個合理預設值。

7. 用 index 和 stride 分配工作

接著把 kernel 改成：

__global__
void add(int n, float *x, float *y)
{
    int index = threadIdx.x;
    int stride = blockDim.x;

    for (int i = index; i < n; i += stride)
        y[i] = x[i] + y[i];
}

如果：

add<<<1, 256>>>(N, x, y);

那：

blockDim.x = 256

所以：

stride = 256;

而不同 thread 的 threadIdx.x 會不同：

Thread 0 → index = 0
Thread 1 → index = 1
Thread 2 → index = 2
...
Thread 255 → index = 255

因此：

Thread 0 → 0,   256, 512, 768, ...
Thread 1 → 1,   257, 513, 769, ...
Thread 2 → 2,   258, 514, 770, ...
...

可以把它記成：

index  = 我從哪裡開始
stride = 我每次往後跳幾格

8. 從單一 Block 擴展到整個 Grid

單一 block 還是無法充分利用整張 GPU。

GPU 裡有很多 Streaming Multiprocessors（SM），而不同 blocks 可以被分配到不同 SM 上執行。

因此要增加 block 數量。

假設：

int blockSize = 256;
int numBlocks = (N + blockSize - 1) / blockSize;

再 launch：

add<<<numBlocks, blockSize>>>(N, x, y);

9. 為什麼 (N + blockSize - 1) / blockSize 要加 blockSize - 1？

這是整數版本的 向上取整除法（Ceiling Division）。

如果：

N = 1000
blockSize = 256

直接：

1000 / 256

C++ 整數除法會得到：

3

但：

3 × 256 = 768

不夠覆蓋 1000 個元素。

所以我們要：

ceil(1000 / 256) = 4

整數寫法：

(N + blockSize - 1) / blockSize

代入：

(1000 + 256 - 1) / 256
= 1255 / 256
= 4

如果剛好整除，例如：

N = 768

則：

(768 + 255) / 256
= 1023 / 256
= 3

不會多算一個 block。

10. Global Thread Index

有多個 blocks 後，只看：

threadIdx.x

已經不夠，因為每個 block 的 thread index 都會重新從 0 開始。

因此要算：

int index =
    blockIdx.x * blockDim.x + threadIdx.x;

附圖可以很直覺地理解：

假設：

blockIdx.x  = 2
blockDim.x  = 256
threadIdx.x = 3

則：

index
= 2 × 256 + 3
= 515

也就是：

Block 2 裡的 Thread 3，在整個 grid 中的 global index 是 515。

11. Grid-Stride Loop

完整 kernel：

__global__
void add(int n, float *x, float *y)
{
    int index =
        blockIdx.x * blockDim.x + threadIdx.x;

    int stride =
        blockDim.x * gridDim.x;

    for (int i = index; i < n; i += stride)
        y[i] = x[i] + y[i];
}

其中：

stride = blockDim.x * gridDim.x;

代表：

整個 Grid 中的總 thread 數。

如果：

blockDim.x = 256
gridDim.x  = 4096

則：

stride = 256 × 4096
       = 1,048,576

這種寫法稱為：

Grid-Stride Loop

它的好處是，即使資料數量比總 thread 數還多，每個 thread 仍然可以繼續處理：

index
index + stride
index + 2 × stride
...

12. 為什麼增加 Blocks 後沒有立刻變快？

這個 Vector Add：

y[i] = x[i] + y[i];

每個元素只有一次加法，但需要：

讀 x[i]
讀 y[i]
寫 y[i]

也就是：

計算量低、Memory Access 多。

因此它本質上很容易變成 Memory-Bandwidth Bound。

也就是：

GPU 算力不是瓶頸
Memory Bandwidth 才是瓶頸

13. Unified Memory

這份程式使用：

cudaMallocManaged(...)

配置 Unified Memory。

它讓 CPU 與 GPU 可以使用相同 pointer：

float *x;
cudaMallocManaged(&x, N * sizeof(float));

但「CPU / GPU 都能存取」不代表資料同時存在兩邊。

Unified Memory 本質上是 Virtual Memory。

某個 page 在某一時間可能 resident 在：

CPU Memory

或：

GPU Memory

14. Page Fault 為什麼會拖慢 GPU？

CPU 一開始初始化：

x[i] = 1.0f;
y[i] = 2.0f;

所以 memory page 一開始在 CPU。

GPU kernel 開始執行後，第一次讀這些 page 時：

GPU 要資料
↓
發現 page 在 CPU
↓
Page Fault
↓
搬移 CPU → GPU
↓
GPU Thread Stall
↓
繼續計算

如果反覆發生，就會變成：

算一下
↓
等資料
↓
算一下
↓
等資料

這會嚴重影響 kernel latency。

15. Prefetching

如果我們知道 kernel 等一下會使用 x 和 y，可以先：

cudaMemPrefetchAsync(
    x,
    N * sizeof(float),
    device
);

cudaMemPrefetchAsync(
    y,
    N * sizeof(float),
    device
);

概念變成：

CPU 初始化
↓
Prefetch
↓
CPU → GPU
↓
資料先準備好
↓
GPU Kernel

因此 kernel 不需要執行到一半才因 page fault 停下來搬資料。

要注意：

Prefetching 並不是讓 CPU → GPU 的資料傳輸時間消失。

真正改善的是：

把 data migration 移出 kernel 的 critical path，避免 GPU 在計算期間因 page fault stall。

16. 最終版本

#include <iostream>
#include <cmath>
#include <cuda_runtime.h>

__global__
void add(int n, float *x, float *y)
{
    int index =
        blockIdx.x * blockDim.x + threadIdx.x;

    int stride =
        blockDim.x * gridDim.x;

    for (int i = index; i < n; i += stride) {
        y[i] = x[i] + y[i];
    }
}

int main()
{
    int N = 1 << 20;

    float *x;
    float *y;

    cudaMallocManaged(&x, N * sizeof(float));
    cudaMallocManaged(&y, N * sizeof(float));

    for (int i = 0; i < N; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    int device = 0;
    cudaGetDevice(&device);

    cudaMemPrefetchAsync(
        x,
        N * sizeof(float),
        device
    );

    cudaMemPrefetchAsync(
        y,
        N * sizeof(float),
        device
    );

    int blockSize = 256;
    int numBlocks =
        (N + blockSize - 1) / blockSize;

    add<<<numBlocks, blockSize>>>(N, x, y);

    cudaDeviceSynchronize();

    cudaMemPrefetchAsync(
        y,
        N * sizeof(float),
        cudaCpuDeviceId
    );

    cudaDeviceSynchronize();

    float maxError = 0.0f;

    for (int i = 0; i < N; i++) {
        maxError =
            fmax(maxError, fabs(y[i] - 3.0f));
    }

    std::cout
        << "Max error: "
        << maxError
        << std::endl;

    cudaFree(x);
    cudaFree(y);

    return 0;
}

17. 效能結果

NVIDIA 教學中的測試結果：

從：

1 Thread

到：

256 Threads

效能大幅提升。

再從：

1 Block

到：

Multiple Blocks

讓更多 SM 同時有工作。

最後達到：

265 GB/s

已經超過 T4 峰值 bandwidth 320 GB/s 的 80%。

因此這個 kernel 最後的瓶頸已經不是 compute，而是：

Memory Bandwidth

18. Day 1 結論

今天最重要的不是記住 CUDA API，而是理解 GPU Programming 的思維轉換：

CPU：
一個 worker
用 loop 處理很多資料

GPU：
大量 workers
每個 thread 負責自己的資料

CUDA 中：

Grid
  ↓
Block
  ↓
Warp
  ↓
Thread

而一個 thread 的全域位置通常由：

blockIdx.x * blockDim.x + threadIdx.x

決定。

最後，GPU 的效能也不只是「thread 越多越快」。

當 Vector Add 已經有足夠的 parallelism 後，瓶頸會從：

Compute

轉成：

Memory Bandwidth

因此 GPU optimization 真正要做的是：

讓計算與資料搬移都能有效利用硬體，而不是單純增加 thread 數量。

Day 1 到這裡，已經可以回答幾個最基本的 CUDA 問題：

• Kernel 是什麼？
• <<<grid, block>>> 是什麼？
• 一個 Block 為什麼常使用 256 Threads？
• 為什麼 Warp 是 32 Threads？
• threadIdx / blockIdx 怎麼算 Global Index？
• Grid-Stride Loop 在解決什麼問題？
• 為什麼 Unified Memory 會有 Page Fault？
• Prefetching 為什麼能大幅降低 Kernel Time？
• 為什麼 Vector Add 最後是 Memory-Bandwidth Bound？

下一篇會進一步把 Thread / Block / Grid 的執行模型弄得更清楚，並開始實際解 GPU 題目。