Stencil 基礎與基本平行核心 (Background & Basic Kernel)

重點總覽 (Overview)

離散化與結構化網格 (Discretization & Structured Grids)

把連續、可微的函數轉成電腦可處理的離散表示 (discrete representation):

(A) 連續函數              (B) 結構化網格 (h = π/6)       (C) 離散陣列 F
 y = sin(x)               x = i·h, 等距格點             F[0..6]
      .-''-.              o   o   o   o   o   o   o     F[2] = sin(2·π/6) = 0.87
    /        \            |   |   |   |   |   |   |     x 值「隱含」= i·h
   0 ........ π           0  π/6 2π/6 ...        π      (不另存 x)

• 格距 h 愈小 → 表示愈精確,但需要更多儲存與更多運算(解 PDE 時)。
• 非格點處的值用 interpolation(linear / splines)近似。
• 精度也取決於浮點位寬:

有限差分與 Stencil 的階數 (Finite Difference & Stencil Order)

Stencil 定義:在 structured grid 每個點上套用的一組幾何權重模式 (geometric pattern of weights),規定如何由中心點與鄰點的值,用數值近似求出該點的輸出。

一階導數的經典有限差分近似:

• 誤差 O(h²) → 與 h 的平方成正比,h 愈小近似愈好(圖 8.1 中 h = π/6 ≈ 0.52)。
• 因為 x = i·h,可寫成對陣列的運算:

→ 係數 [-1/2h, 0, 1/2h] 定義了一個 1D three-point stencil。

Order(階數)= 中心點「每一側」的格點數 = 近似導數的階數:

1D stencils (● = 中心, ○ = 鄰點)
 order 1 (3-point):        ○──●──○            (含至一階導數)
 order 2 (5-point):     ○──○──●──○──○         (含至二階導數)
 order 3 (7-point):  ○──○──○──●──○──○──○      (含至三階導數)

通則: 1D 點數 = 2·order + 1;含至 n 階導數 → 每側 n 點

基本平行 Stencil Sweep 核心 (Basic Parallel Kernel)

Stencil sweep:對所有相關輸入格點套用 stencil,一次產生所有輸出格點。基本核心的簡化假設:

1. 一次 sweep 內,輸出格點之間無相依(可全平行)。
2. 邊界格點存 boundary condition,輸出 = 輸入,不更新(只算內部陰影區)。

Output grid (2D 示意,每個 block 算 4×4 output tile)
┌─────────────────────────┐
│ □ □ □ □ □ □ □ □ □ □ □ □ │  □ = boundary cell:輸出 = 輸入(thread 關閉)
│ □ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ □ │  ▓ = 內部格點:由 stencil 計算
│ □ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ □ │
│ □ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ ▓ □ │  邊界形成 grid 表面的一層
│ □ □ □ □ □ □ □ □ □ □ □ □ │
└─────────────────────────┘

3D seven-point stencil(中心 + x/y/z 各 2 鄰點):

            in[i+1][j][k]   (z+)
                 │
   in[i][j-1][k] │ in[i][j+1][k]
         (y-)  \ │ /  (y+)
in[i][j][k-1]──[center]──in[i][j][k+1]    center = in[i][j][k]
       (x-)   / │ \   (x+)
                 │
            in[i-1][j][k]   (z-)
       7 loads → c0·center + c1..c6·neighbors

Thread → data 映射(3D,row-major 線性化,見 03-Multidimensional-Grids-And-Data/02-Mapping-Threads-to-Multidimensional-Data):

__global__ void stencil_kernel(float* in, float* out, unsigned int N) {
    // i = z(最外/慢), j = y, k = x(最內/快)
    unsigned int i = blockIdx.z*blockDim.z + threadIdx.z;
    unsigned int j = blockIdx.y*blockDim.y + threadIdx.y;
    unsigned int k = blockIdx.x*blockDim.x + threadIdx.x;

    // 只算內部格點;邊界 (index 0 或 N-1) 保持不變
    if (i >= 1 && i < N-1 && j >= 1 && j < N-1 && k >= 1 && k < N-1) {
        out[i*N*N + j*N + k] =  c0 * in[ i   *N*N +  j   *N +  k  ]   // center
                              + c1 * in[ i   *N*N +  j   *N + (k-1)]   // x-
                              + c2 * in[ i   *N*N +  j   *N + (k+1)]   // x+
                              + c3 * in[ i   *N*N + (j-1)*N +  k  ]    // y-
                              + c4 * in[ i   *N*N + (j+1)*N +  k  ]    // y+
                              + c5 * in[(i-1)*N*N +  j   *N +  k  ]    // z-
                              + c6 * in[(i+1)*N*N +  j   *N +  k  ];   // z+
    }
}
// c0..c6 由所解的微分方程決定;線性索引 = i*N*N + j*N + k

算術強度:0.46 OP/B (Arithmetic Intensity)

對基本核心計算 floating-point 對 global memory access 比值:

Roofline 直覺:
 compute-bound  ┤            ╱──────  峰值算力
                │          ╱
                │        ╱
                │      ╱   ← kernel 落在這條斜線上(memory-bound)
 memory-bound   │    ╱ ▲ 0.46 OP/B (遠在左側)
                └──────┴──────────────► OP/B
                      0.46

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

Related Notes

• 08-Stencil/02-Shared-Memory-Tiling-for-Stencil
• 08-Stencil/03-Thread-Coarsening-and-Register-Tiling
• 07-Convolution/01-Convolution-Fundamentals-and-Basic-Kernel
• 05-Memory-Architecture-And-Data-Locality/01-Memory-Access-Efficiency-and-CUDA-Memory-Types
• 03-Multidimensional-Grids-And-Data/02-Mapping-Threads-to-Multidimensional-Data
• 24-Numerical-Considerations/01-Floating-Point-Data-Representation