FHD Kernel 的記憶體頻寬優化：Registers 與 Constant Memory

重點總覽 (Overview)

承接 17-Iterative-MRI-Reconstruction/02-FHD-Kernel-Parallelism-Structure 選定的 gather kernel（Fig. 17.10），本節是 Step 2：把這個 kernel 從記憶體頻寬綁定 (memory-bandwidth bound) 一路救回運算綁定邊緣。手法是逐步抬高 compute-to-global-memory-access ratio (OP/B，即 arithmetic intensity)。

瓶頸分析：Compute-to-Memory Ratio (Arithmetic Intensity)

逐項數 Fig. 17.10 內層迴圈每次迭代的 global memory 存取（共 14 次，每次 4 bytes = float）：

kx[m] ky[m] kz[m]        ← 3  (k-space 座標)
x[n]  y[n]  z[n]         ← 3  (voxel 座標)
rMu[m] (×2) iMu[m] (×2)  ← 4  (Mu，編譯器看成 2 個位置)
rFhD[n] read+write       ← 2
iFhD[n] read+write       ← 2
                         ───
                          14 accesses

約 13 個浮點運算（乘/加/三角）：

第一招：把陣列元素提升到 Registers (Register Promotion)

判斷哪些 global 存取可以變成 register：看 index 與 threadIdx 的關係。

__global__ void cmpFhD(/* ... */, int M) {
    int n = blockIdx.x * FHD_THREADS_PER_BLOCK + threadIdx.x;
    // 頻繁存取的座標與輸出 → automatic variables (registers)
    float xn_r = x[n]; float yn_r = y[n]; float zn_r = z[n];
    float rFhDn_r = rFhD[n]; float iFhDn_r = iFhD[n];
    for (int m = 0; m < M; m++) {
        float expFhD = 2*PI*(kx[m]*xn_r + ky[m]*yn_r + kz[m]*zn_r);
        float cArg = cos(expFhD);
        float sArg = sin(expFhD);
        rFhDn_r += rMu[m]*cArg - iMu[m]*sArg;   // 累加在 register
        iFhDn_r += iMu[m]*cArg + rMu[m]*sArg;
    }
    rFhD[n] = rFhDn_r; iFhD[n] = iFhDn_r;       // 圈後寫回一次
}

thread n 的 m-loop（register 重用）:

  x[n],y[n],z[n] ── 載入 1 次 ──► [xn_r][yn_r][zn_r]  (整圈唯讀重用)
                                        │
   m=0 ─► m=1 ─► m=2 ─► ... ─► m=M-1    │ 全程從 register 取
                                        ▼
  rFhDn_r / iFhDn_r  累加於 register ───► 圈結束才寫回 global 1 次

效果：每圈 global accesses 14 → 7，OP/B 0.23 → 0.46（13/(7×4)）。

第二招：k-space 放 Constant Memory + Chunking

kx[m] ky[m] kz[m] 為何適合 constant memory（而非 register）？

• 唯讀：kernel 不修改 k-space。
• index = m，與 threadIdx 無關 → warp 內 32 個 thread 同一時刻存取同一個元素。

constant cache 廣播 (broadcast)：
   kx_c[m] ──取一次進 cache──► broadcast 給 warp 內全部 32 threads
   ⇒ 每 32 次存取，≥31 次由 cache 服務 ⇒ 消除 96%+ 的 global accesses

容量限制 → Chunking 64KB

Constant memory 容量僅 64KB，但 k-space 樣本可達數百萬個 → 必須切塊 (chunk)，host 端多次呼叫 kernel，每次先把一塊塞進 constant memory：

__constant__ float kx_c[CHUNK_SIZE], ky_c[CHUNK_SIZE], kz_c[CHUNK_SIZE];

void main() {
    for (int i = 0; i < M/CHUNK_SIZE; i++) {
        cudaMemcpyToSymbol(kx_c, &kx[i*CHUNK_SIZE], 4*CHUNK_SIZE, cudaMemcpyHostToDevice);
        cudaMemcpyToSymbol(ky_c, &ky[i*CHUNK_SIZE], 4*CHUNK_SIZE, cudaMemcpyHostToDevice);
        cudaMemcpyToSymbol(kz_c, &kz[i*CHUNK_SIZE], 4*CHUNK_SIZE, cudaMemcpyHostToDevice);
        cmpFhD<<<N/FHD_THREADS_PER_BLOCK, FHD_THREADS_PER_BLOCK>>>(/* ... */, CHUNK_SIZE);
    }
    // 若 M 不是 CHUNK_SIZE 的整數倍，需再補一輪 memcpy + kernel
}

kernel 內 kx[m] 改成 global __constant__ 變數 kx_c[m]（不再放參數表）。剩下唯一的 global 存取是 rMu[m]、iMu[m]（編譯器把 4 次合併成 2 次 load）：

第三招：Array-of-Structures (AoS) 佈局調整 Cache

跑完 Fig. 17.13 在某些裝置上加速不如預期——問題出在 constant cache 設計 + k-space 記憶體佈局。

Cache line 會一次載入多個連續字。三個分離陣列 kx ky kz → 每次迭代需要 3 條 cache line（分屬三個陣列的不同區段）。

Working set 爆掉 → cache thrashing

1024 threads/block, 2 blocks/SM 同時常駐：
  warps/SM = (1024/32) × 2 = 64 warps
  每 warp 需 3 條 cache line（kx/ky/kz 各一）
  worst case working set = 64 × 3 = 192 cache lines
  (即使 3 個 warp 共享同一迭代) ≈ 64 cache lines

但有些裝置 constant cache 只有 ~32 條 line → 容不下 working set → 不同 warp 互相搶 line：下一圈要用的元素早被別的 warp 擠出去 → cache 無法消除 global accesses。

解法：SoA → AoS

把 x/y/z 三分量塞進一個 struct 的三個欄位，強制它們連續存放 → 一次迭代的 3 分量塞進同一條 cache line：

(A) Separate arrays (SoA)            (B) Array of structures (AoS)
  kx: [..x[i]..]  ─ line 1            k: [ x y z | x y z | x y z ... ]
  ky: [..y[i]..]  ─ line 2               └──┬──┘
  kz: [..z[i]..]  ─ line 3            一次迭代的 x,y,z ─► 同一條 cache line
  → 每 warp 需 3 lines                 → 每 warp 只需 1 line (working set ÷3)

struct kdata { float x, y, z; };
__constant__ struct kdata k[CHUNK_SIZE];
// kernel 內：kx[m] → k[m].x, ky[m] → k[m].y, kz[m] → k[m].z
float expFhD = 2*PI*(k[m].x*xn_r + k[m].y*yn_r + k[m].z*zn_r);

附帶好處：只剩一個陣列 → 只需 一次 cudaMemcpyToSymbol，搬移大小由 4*CHUNK_SIZE 改為 12*CHUNK_SIZE（一次搬三分量）。

考試/面試重點 (Exam / Test Patterns)

Related Notes

• 17-Iterative-MRI-Reconstruction/02-FHD-Kernel-Parallelism-Structure
• 17-Iterative-MRI-Reconstruction/04-Hardware-Trig-Accuracy-and-Performance-Tuning
• 17-Iterative-MRI-Reconstruction/01-MRI-Background-and-Iterative-Reconstruction
• 05-Memory-Architecture-And-Data-Locality/01-Memory-Access-Efficiency-and-CUDA-Memory-Types
• 07-Convolution/02-Constant-Memory-and-Caching
• 06-Performance-Considerations/01-Memory-Coalescing