Hardware 三角函數、精度驗證與效能調校
重點總覽 (Overview)
本筆記接續 17-Iterative-MRI-Reconstruction/03-FHD-Memory-Bandwidth-Optimization，涵蓋 cmpFhD 最佳化的最後兩步 (Step 3、Step 4) 與全章總結。
| 項目 | 手段 | 關鍵數字 / 結論 |
|---|---|---|
| Step 3 硬體三角函數 | sin()/cos() → __sin()/__cos() (SFU intrinsics) |
在熱迴圈內，throughput 大增 |
| 精度驗證 | 與 phantom 比對 PSNR | HW trig：27.6 → 27.5 dB (可接受) |
| 精度對照 | double / single / HW trig | single 與 double 同為 27.6 dB；gridding/iFFT 僅 16.8 dB |
| Step 4 效能調校 | threads/block + #pragma unroll |
全組合窮舉比啟發式再快 ~20% |
| 調校張力 | unroll ↑ → 指令額外負擔 ↓ 但 register ↑ → occupancy ↓ | 參數須聯合評估 |
| 全章加速 | scatter→gather + register + constant + HW trig | 基礎版 → 最佳版 ~10× |
| Amdahl 後果 | FHD 從 ~100% 降到 ~50% | 另 50% 變成 CG solver 主導 |
SFU 硬體三角函數 (**sin /**cos Intrinsics)
- CUDA 提供數學函式的硬體實作，throughput 遠高於軟體版；最初動機是加速繪圖中的視角轉換。
- 這些函式由 SFU (Special Function Units) 以硬體指令執行。
- 使用方式極簡：把
sin()/cos()改成__sin()/__cos()(前綴兩個底線)。它們是 compiler 認得的 intrinsic function，會被翻成 SFU 指令。 - 因為呼叫位於高度執行的迴圈本體，改動帶來顯著加速。
#define FHD_THREADS_PER_BLOCK 1024
__global__ void cmpFhD(/* ... */ int M) {
int n = blockIdx.x * FHD_THREADS_PER_BLOCK + threadIdx.x;
float xn_r = x[n], yn_r = y[n], zn_r = z[n]; // 入 register
float rFhDn_r = rFhD[n], iFhDn_r = iFhD[n];
for (int m = 0; m < M; m++) {
float expFhD = 2*PI*(k[m].x*xn_r + k[m].y*yn_r + k[m].z*zn_r);
float cArg = __cos(expFhD); // SFU intrinsic
float sArg = __sin(expFhD); // SFU intrinsic
rFhDn_r += rMu[m]*cArg - iMu[m]*sArg;
iFhDn_r += iMu[m]*cArg + rMu[m]*sArg;
}
rFhD[n] = rFhDn_r; iFhD[n] = iFhDn_r;
}
為何 SFU 是關鍵瓶頸對策：在原本的迴圈裡，浮點算術對三角函數的比例只有 13:2，若不降低 trig 成本，效能會被 sin/cos 的長延遲、低吞吐主宰。
每個 thread 的 inner m-loop:
expFhD = 2*PI*(k.x*xn + k.y*yn + k.z*zn) ──▶ FP32 cores (FMA pipeline)
│
├─ __cos(expFhD) ─┐
├─ __sin(expFhD) ─┼──▶ SFU (獨立硬體單元，與 FP pipeline 並行)
▼ ▼
rFhDn += rMu*cArg - iMu*sArg ──▶ FP32 cores (可與 SFU 重疊)
iFhDn += iMu*cArg + rMu*sArg
| 面向 | 軟體版 sin()/cos() |
硬體版 __sin()/__cos() |
|---|---|---|
| 實作 | math library，用約 13/12 個 FP op 的 5 項 Taylor 級數 | 單一 SFU 指令 |
| 執行單元 | 一般 FP32 cores | 專屬 SFU |
| Throughput | 低 | 高 |
| 精度 | 高 | 略低 (見 CUDA C Programming Guide) |
計算密度 (compute-to-memory ratio) 在此區間是 best 0.75 OP/B / worst 0.25 OP/B：best case 假設 trig 用 13/12 個 FP op 的 Taylor 級數，worst case 假設每個 trig 是單一硬體 op。改用 SFU 反而把 FP-op 計數推向 worst-case 0.25 OP/B，但整體更快——因為 SFU 把 trig 從 FP pipeline 卸載、且本身吞吐高。OP/B 數字不能脫離「哪個單元在算」單獨解讀。
精度驗證 (PSNR-Based Accuracy Validation)
從軟體函式換成硬體函式必須驗證精度損失。HW 實作目前精度低於 software library，醫學影像不能盲改。
驗證流程 (phantom test)：用一個虛構物體的「完美」影像 I0 (稱為 phantom)，反向合成出對應的「掃描」k-space 資料，再丟進重建系統得到重建影像 I，最後把 I0 與 I 的 voxel 值代入 PSNR 公式比較。
phantom (I0) ──反向合成──▶ 合成 k-space data
│ │
│ 重建 (cmpFhD + CG)
│ ▼
└──────── PSNR(I0, I) ◀── reconstructed image (I)
- 通過門檻依應用而定；本案與臨床 MRI 專家確認 HW 函式造成的 PSNR 變化在可接受範圍內。
- 供醫師判讀的影像還需目視檢查影像品質 (artifacts)。
| 重建/實作方式 | PSNR | 說明 |
|---|---|---|
| Iterative recon, CPU double | 27.6 dB | 基準 |
| Iterative recon, CPU single | 27.6 dB | 與 double 無可測量差異 |
Iterative recon, HW trig (__sin/__cos) |
27.5 dB | 僅輕微下降，仍可接受 |
| 簡單 gridding / iFFT (bilinear interp) | 16.8 dB | 明顯較差、嚴重 artifacts |
此表同時佐證兩件事：(1) single-precision 對本應用足夠；(2) iterative reconstruction (27.6 dB) 比傳統 gridding/iFFT (16.8 dB) 大幅減少 artifacts，這正是 17-Iterative-MRI-Reconstruction/01-MRI-Background-and-Iterative-Reconstruction 提出該方法的動機。
實驗性效能調校 (Experimental Performance Tuning)
兩個 kernel 設定參數需要決定 (Step 4)：
| 參數 | 作用 | 增大太多的代價 |
|---|---|---|
threads / block (FHD_THREADS_PER_BLOCK) |
要足夠填滿每個 SM 的 thread 容量 | 與其他資源競爭 |
loop unroll 次數 (#pragma unroll N) |
減少 overhead 指令、降低每筆 k-space 的 clock cycle | register 使用上升 → 每 SM 可容納 block 數下降 (occupancy ↓) |
for (int m = 0; m < M; m++) {
#pragma unroll 4 // 在 line 07 的 for-loop 前指定展開次數
// ... loop body ...
}
增大一個參數會吃掉本可給另一參數的資源 (register / occupancy)。因此要以實驗方式聯合評估，組合數可能很多。
- 對 FHD：窮舉所有組合選最佳實測 runtime，比只探索部分趨勢的啟發式搜尋再快 ~20%。
- Ryoo et al. (2008) 提出以 Pareto-optimal curve 為基礎的方法，先篩掉多數劣等組合 (program optimization carving)。
整體加速與 CG/Amdahl 總結 (Overall ~10x Speedup & CG/Amdahl Summary)
最佳化進程 (compute-to-global-memory-access ratio)：
| 版本 | 關鍵技巧 | OP/B |
|---|---|---|
| Gather 基礎版 (Fig 17.10) | loop fission + interchange，避開 atomics | 0.23 |
| + Registers (Fig 17.11) | 把 x/y/z[n]、rFhD/iFhD[n] 放進 register |
0.46 |
| + Constant memory (Fig 17.13/16) | kx/ky/kz 入 constant cache + array-of-structs 佈局 |
1.63 |
| + HW trig (Fig 17.17) | __sin/__cos (SFU) |
解除 trig 吞吐瓶頸 |
從基礎版到最終最佳化版約有 10× 加速。
- 最佳化之前：FHD 幾乎佔 100% 執行時間。
- 最佳化之後：FHD 只剩約 50%，另約 50% 落在原本 <1% 的 CG solver ("find ρ")。
- 這是平行化真實應用的常見現象 (Amdahl's Law)：某階段被大幅加速後，執行時間轉由原本不顯著的其他階段主宰；要再快就得去加速 CG solver。
整章方法論摘要：scatter → gather (避開 global atomics) 是核心，靠 loop fission + loop interchange 達成；再疊加 register 提升、constant memory/cache、HW 函式三類最佳化。
考試/面試重點 (Exam / Test Patterns)
| 情境 / 關鍵字 | 答案 / 技巧 |
|---|---|
__sin vs sin 差別 |
__sin 是 intrinsic，編成 SFU 硬體指令，吞吐高但精度略低 |
| 為何 HW trig 仍是必要最佳化 | FP 算術:trig = 13:2；不降 trig 成本，效能被其長延遲/低吞吐主宰 |
| HW trig 對 OP/B 的影響 | 把 trig 變單一 op，OP/B 反而趨近 worst-case 0.25，但因 SFU 卸載+高吞吐而更快 |
| 如何驗證精度沒壞 | 用 phantom (I0) 反向合成 k-space → 重建 I → 算 PSNR(I0,I)；必要時目視檢查 |
| PSNR 公式 | PSNR = 20·log10( max(I0) / sqrt(MSE) )，MSE 為均方誤差 |
| 數字記憶 | double=single=27.6 dB；HW trig=27.5 dB；gridding/iFFT=16.8 dB |
#pragma unroll 的取捨 |
減少 overhead 指令，但 register ↑ → occupancy ↓；需與 threads/block 聯合調 |
| 為何參數要聯合窮舉 | 參數共享 register/occupancy 資源，互相耦合；FHD 窮舉比啟發式快 ~20% |
| 全章總加速 / Amdahl | ~10×；FHD 由 ~100% 降到 ~50%，CG solver 成新瓶頸 |
Related Notes
- 17-Iterative-MRI-Reconstruction/03-FHD-Memory-Bandwidth-Optimization
- 17-Iterative-MRI-Reconstruction/02-FHD-Kernel-Parallelism-Structure
- 17-Iterative-MRI-Reconstruction/01-MRI-Background-and-Iterative-Reconstruction
- 19-Parallel-Programming-And-Computational-Thinking/01-Goals-of-Parallel-Computing-and-Amdahls-Law
- 04-Compute-Architecture-And-Scheduling/03-Resource-Partitioning-and-Occupancy
- 24-Numerical-Considerations/02-Representable-Numbers-Precision-and-Accuracy