Thread Coarsening 與 Register Tiling
重點總覽 (Overview)
針對 3D stencil，shared memory tiling 受限於 block size ≤ 1024 與 shared memory 容量 ∝ T³，使得 tile 邊長 T 只能取到 8，reuse 與 coalescing 都偏差。本節用 thread coarsening（z 方向） 解除這個限制，再用 register tiling 把只給單一 thread 使用的 z-鄰居搬進 register，紓解 shared memory 壓力。
| 版本 (Kernel) | Block size | On-chip 用量 | 可用 T | OP/B (七點 stencil) |
|---|---|---|---|---|
| Basic sweep (Fig 8.6) | T³ threads | 無 tiling | — | 0.46 |
| Shared-mem tiling (Fig 8.8) | T³ threads | shared T³ 元素 |
實際上限 T=8 | 1.37 (T=8) |
| Coarsening z (Fig 8.10) | T² threads | shared 3T² 元素 |
T=32 | 2.68 (T=32) |
| Coarsening + Register tiling (Fig 8.12) | T² threads | shared T² + 3 regs/thread |
T=32 | 2.68（同上，但更快） |
| 理論上限 (T→∞) | — | — | — | 3.25 |
Thread coarsening 提高 OP/B（讓 T 變大、降低 halo 比例）；register tiling 不改變 OP/B 也不改變 global memory 流量，它只是把 on-chip 的 reuse 從 shared memory 移到更快的 register，並省下 2/3 的 shared memory。
Z 方向 Thread Coarsening (Thread Coarsening in the Z Direction)
動機：3D seven-point stencil 的 tiled kernel 上限 3.25 OP/B，但硬體把 block 限制在 1024 threads，且 shared memory ∝ T³，逼得 T=8（512 threads）。此時：
- Halo 佔比高：8×8×8 input tile = 512，output tile 6³ = 216，halo = 296，約 58%（對比 2D 32×32 conv 僅約 12%）→ reuse 與 OP/B 被壓到 1.37。
- Coalescing 差：8 寬的列，一個 32-thread warp 要橫跨至少 4 個相距很遠的 row，無法合併存取，浪費 DRAM 頻寬。
核心想法：把每個 thread 的工作從「算一個 grid point」加重為「算 z 方向一整條 column」。Block 只覆蓋 input tile 的一個 x-y plane（T² threads），沿 z 逐層 iterate，每次算出一個 output plane。
peeled-away 視角 (T=6 → input plane 6×6, output plane 4×4)
z (into page)
│ input tile = T³ = 6³ = 216 grid points
│ output tile = (T-2)³ = 4³ = 64 grid points
▼
x ──►
┌───────────────┐ ← block = 1 個 input x-y plane = 6×6 threads
y │ . . . . . . . │ 其中內層 4×4 thread 才真正算 output
│ │ . ┌───────┐ . │
▼ │ . │ o o o │ . │ o = active output thread (4×4)
│ . │ o o o │ . │ . = halo-loading thread (周邊一圈)
│ . └───────┘ . │
└───────────────┘
block 沿 z 反覆 iterate：每個 iteration 算一層 output plane
每個 iteration 只需要 3 個 input plane 同時在片上：inPrev_s(z-1)、inCurr_s(z)、inNext_s(z+1)。算完一層後把它們**輪替（rotate）**再往 z 前進一格：
iteration i iteration i+1
inPrev_s ← plane z-1 inPrev_s ← (舊 inCurr_s) plane z
inCurr_s ← plane z inCurr_s ← (舊 inNext_s) plane z+1
inNext_s ← plane z+1 inNext_s ← 新載入 plane z+2
│ shift ──────────────────►│
// Fig 8.10: thread coarsening in z, 3D seven-point stencil
// IN_TILE_DIM = T, OUT_TILE_DIM = T-2
__global__ void stencil_kernel(float* in, float* out, unsigned int N) {
int iStart = blockIdx.z * OUT_TILE_DIM; // z 起點（output）
int j = blockIdx.y * OUT_TILE_DIM + threadIdx.y - 1; // y（-1 = 含 halo）
int k = blockIdx.x * OUT_TILE_DIM + threadIdx.x - 1; // x（-1 = 含 halo）
__shared__ float inPrev_s[IN_TILE_DIM][IN_TILE_DIM];
__shared__ float inCurr_s[IN_TILE_DIM][IN_TILE_DIM];
__shared__ float inNext_s[IN_TILE_DIM][IN_TILE_DIM];
// 預載前兩層 (z-1, z)
if (iStart-1 >= 0 && iStart-1 < N && j >= 0 && j < N && k >= 0 && k < N)
inPrev_s[threadIdx.y][threadIdx.x] = in[(iStart-1)*N*N + j*N + k];
if (iStart >= 0 && iStart < N && j >= 0 && j < N && k >= 0 && k < N)
inCurr_s[threadIdx.y][threadIdx.x] = in[ iStart *N*N + j*N + k];
for (int i = iStart; i < iStart + OUT_TILE_DIM; ++i) {
// 載入第三層 (z+1)
if (i+1 >= 0 && i+1 < N && j >= 0 && j < N && k >= 0 && k < N)
inNext_s[threadIdx.y][threadIdx.x] = in[(i+1)*N*N + j*N + k];
__syncthreads();
// 邊界 + halo thread 關閉
if (i >= 1 && i < N-1 && j >= 1 && j < N-1 && k >= 1 && k < N-1) {
if (threadIdx.x >= 1 && threadIdx.x < IN_TILE_DIM-1 &&
threadIdx.y >= 1 && threadIdx.y < IN_TILE_DIM-1) {
out[i*N*N + j*N + k] =
c0*inCurr_s[threadIdx.y][threadIdx.x]
+ c1*inCurr_s[threadIdx.y][threadIdx.x-1]
+ c2*inCurr_s[threadIdx.y][threadIdx.x+1]
+ c3*inCurr_s[threadIdx.y-1][threadIdx.x]
+ c4*inCurr_s[threadIdx.y+1][threadIdx.x]
+ c5*inPrev_s[threadIdx.y][threadIdx.x] // z-1
+ c6*inNext_s[threadIdx.y][threadIdx.x]; // z+1
}
}
__syncthreads();
// rotate：往 z 前進一格
inPrev_s[threadIdx.y][threadIdx.x] = inCurr_s[threadIdx.y][threadIdx.x];
inCurr_s[threadIdx.y][threadIdx.x] = inNext_s[threadIdx.y][threadIdx.x];
}
}
收益（兩件事）：
- Block size 從 T³ 降為 T² → 可取 T=32（1024 threads），OP/B 提升到 2.68。
- 任何時刻只需 3 層在片上 → shared memory 從
T³降為3T²元素。T=32 時3×32²×4B = 12KB/block（合理）。
Coarsening 的本質（見 06-Performance-Considerations/03-Thread-Coarsening）：把平行單位「部分序列化」進單一 thread，換取較少的重複代價——這裡換掉的是「每個 block 重複載入 halo」的低 reuse。
Register Tiling（把 z-鄰居放進 Register）
關鍵觀察：在 Fig 8.3 的所有 stencil（只取 x/y/z 軸鄰居、不含角落）中：
inPrev/inNext（z 鄰居）只被同 (x,y) 的單一 thread 使用 → 不需要放 shared memory。- 只有
inCurr的 x-y 鄰居會被多個 thread 共享 → 唯一真正需要 shared memory 的部分。
於是把 inPrev_s、inNext_s 改成 register 變數 inPrev、inNext，僅保留 inCurr_s：
Fig 8.10 (coarsening) Fig 8.12 (+ register tiling)
plane z-1 inPrev_s [shared] ──► inPrev [register]
plane z inCurr_s [shared] ──► inCurr_s [shared] (+ inCurr 一份 register)
plane z+1 inNext_s [shared] ──► inNext [register]
shared 用量 3·T² T² (=1/3)
// Fig 8.12: thread coarsening + register tiling
__global__ void stencil_kernel(float* in, float* out, unsigned int N) {
int iStart = blockIdx.z * OUT_TILE_DIM;
int j = blockIdx.y * OUT_TILE_DIM + threadIdx.y - 1;
int k = blockIdx.x * OUT_TILE_DIM + threadIdx.x - 1;
float inPrev; // register（z-1）
__shared__ float inCurr_s[IN_TILE_DIM][IN_TILE_DIM]; // 只剩這一層 shared
float inCurr; // register（z）
float inNext; // register（z+1）
if (iStart-1 >= 0 && iStart-1 < N && j >= 0 && j < N && k >= 0 && k < N)
inPrev = in[(iStart-1)*N*N + j*N + k];
if (iStart >= 0 && iStart < N && j >= 0 && j < N && k >= 0 && k < N) {
inCurr = in[iStart*N*N + j*N + k];
inCurr_s[threadIdx.y][threadIdx.x] = inCurr; // 同時保留一份在 shared 供共享
}
for (int i = iStart; i < iStart + OUT_TILE_DIM; ++i) {
if (i+1 >= 0 && i+1 < N && j >= 0 && j < N && k >= 0 && k < N)
inNext = in[(i+1)*N*N + j*N + k];
__syncthreads();
if (i >= 1 && i < N-1 && j >= 1 && j < N-1 && k >= 1 && k < N-1) {
if (threadIdx.x >= 1 && threadIdx.x < IN_TILE_DIM-1 &&
threadIdx.y >= 1 && threadIdx.y < IN_TILE_DIM-1) {
out[i*N*N + j*N + k] =
c0*inCurr // 自己（register）
+ c1*inCurr_s[threadIdx.y][threadIdx.x-1] // x-y 鄰居（shared）
+ c2*inCurr_s[threadIdx.y][threadIdx.x+1]
+ c3*inCurr_s[threadIdx.y-1][threadIdx.x]
+ c4*inCurr_s[threadIdx.y+1][threadIdx.x]
+ c5*inPrev // z-1（register）
+ c6*inNext; // z+1（register）
}
}
__syncthreads();
inPrev = inCurr; // register rotate
inCurr = inNext;
inCurr_s[threadIdx.y][threadIdx.x] = inNext; // 更新 shared 那一層
}
}
兩大好處 vs Fig 8.10：
- 大量 shared memory 讀寫改成 register 存取（latency 更低、bandwidth 更高）→ 更快。
- 每個 block 只用 1/3 的 shared memory（從
3T²→T²）。
代價：每個 thread 多用 3 個 register → 32×32 block 多 3072 registers/block。Stencil order 越高、z 方向鄰居越多，register 用量越大。
Register 是稀缺資源，過度使用會壓低 occupancy（見 04-Compute-Architecture-And-Scheduling/03-Resource-Partitioning-and-Occupancy）。若 register 變成瓶頸，可把部分 plane 改回 shared memory——這是典型的 shared memory ↔ register 取捨。
Register tiling 不是新技巧：matrix multiplication / convolution kernel 早就把每個 thread 算出的 output 值存在 register，整個 block 的 output tile 就是「集體存在 block 各 thread register 裡的 tile」。這裡只是首次拿 register 來存 input tile 的一部分——同一塊 tile 在運算過程中時而在 register、時而在 shared memory。
OP/B 與資源公式 (Quantitative Ratios)
| 量 | 公式 | 數值範例 |
|---|---|---|
| Tiled OP/B（七點 stencil） | 13·(T-2)³ / (4·T³) = (13/4)·(1 - 2/T)³ |
T=8→1.37；T=32→2.68 |
| OP/B 上限 (T→∞) | 13/4 |
3.25 OP/B |
| Shared mem（Fig 8.10 coarsening） | 3·T² 元素 |
T=32 → 3·32²·4B = 12KB |
| Shared mem（Fig 8.12 register tiling） | T² 元素 = 1/3 |
T=32 → 4KB |
| 額外 register（七點） | 3 regs/thread | 32×32 block → 3072 regs/block |
| Halo 佔比（order-1 3D, 邊長 T） | 1 - ((T-2)/T)³ |
T=8 → 58% |
七點 stencil 一次做 13 FLOP（7 乘 + 6 加）、讀 7 個 4-byte input；basic kernel 的 0.46 OP/B 即 13/(7·4)。Coarsening 讓 T 變大、halo 比例下降，是把 OP/B 推向 3.25 的手段；register tiling 只是讓既有 reuse 跑得更快。
考試/面試重點 (Exam / Test Patterns)
| 情境 / 關鍵字 | 答案 / 技巧 |
|---|---|
| 「為何 3D stencil 不能用大 tile？」 | block ≤ 1024 且 shared memory ∝ T³，逼 T=8，halo 佔 58%，OP/B 只 1.37 |
| Thread coarsening 用在哪個方向、解決什麼 | z 方向；block 降為 T² threads，可取 T=32，提高 reuse 與 coalescing |
| Coarsening 後 shared memory 需求 | 3T² 元素（只存 3 層）；T=32 → 12KB |
| 為何 z 鄰居可放 register | inPrev/inNext 只被單一 thread 用，不需跨 thread 共享 |
為何 inCurr 必須留在 shared memory |
它的 x-y 鄰居被多個 thread 共享 |
| Register tiling 對 OP/B / global 流量的影響 | 皆不變；只把 reuse 從 shared memory 移到 register，省 2/3 shared mem |
| Register tiling 的代價 | 多 3 regs/thread（七點）；高階 stencil 更多 → 可能壓 occupancy |
| Shared mem vs register 取捨 | register 不夠時把部分 plane 改回 shared memory |
| 計算 T=32 七點 stencil OP/B | (13/4)(1-2/32)³ = 3.25·0.824 ≈ 2.68 OP/B |
| 練習題 2（32×32 block、coarsen=16） | input tile 32·32·18=18432；output 30·30·16=14400；OP/B 13·14400/(4·18432)≈2.54；shared 無 reg-tiling 3·32²·4B=12KB，有 reg-tiling 32²·4B=4KB |
| 為何只能 register-tile「軸向」stencil | 五/七/十三點等不含角落，z 鄰居才會被單一 thread 獨用 |
Related Notes
- 08-Stencil/01-Stencil-Background-and-Basic-Kernel
- 08-Stencil/02-Shared-Memory-Tiling-for-Stencil
- 06-Performance-Considerations/03-Thread-Coarsening
- 05-Memory-Architecture-And-Data-Locality/02-Tiling-and-Tiled-Matrix-Multiplication
- 07-Convolution/03-Tiled-Convolution-and-Halo-Handling
- 17-Iterative-MRI-Reconstruction/03-FHD-Memory-Bandwidth-Optimization