Tiled Convolution 與 Halo Cells 處理

重點總覽 (Overview)

輸入塊與輸出塊 (Input Tile vs Output Tile)

• Output tile：每個 block 負責計算的 P 元素集合。
• Input tile：要算出該 output tile，需要載入的全部 N 元素——必須往四周各延伸 FILTER_RADIUS，才能涵蓋 output tile 邊緣元素所需的鄰居。
• 延伸出來的外圈即 halo cells;它們同時是相鄰 block 的內部元素(這正是 §7.5 cache 方案的依據)。

Input tile (左)  =  output tile  + 四周各 r 的 halo
            r=2 範例：output 4x4 → input 8x8

  H H H H H H H H        H = halo cell (鄰居 tile 的內部元素)
  H H H H H H H H        O = output tile 元素 (此 block 真正算的)
  H H O O O O H H
  H H O O O O H H        IN_TILE_DIM = OUT_TILE_DIM + 2*FILTER_RADIUS
  H H O O O O H H
  H H O O O O H H
  H H H H H H H H
  H H H H H H H H

兩種 Thread 組織 (Two Thread Organizations)

本章正文採用方案 A。

方案 A:Tiled Kernel 與外圈停用 (Deactivated Exterior Threads)

#define IN_TILE_DIM 32
#define OUT_TILE_DIM ((IN_TILE_DIM) - 2*(FILTER_RADIUS))
__constant__ float F_c[2*FILTER_RADIUS+1][2*FILTER_RADIUS+1];

__global__ void convolution_tiled_2D_const_mem_kernel(float *N, float *P,
                                                      int width, int height) {
    // 1) 每個 thread 對應的「載入 / 計算」座標 (-r 偏移到 input tile 左上角)
    int col = blockIdx.x*OUT_TILE_DIM + threadIdx.x - FILTER_RADIUS;
    int row = blockIdx.y*OUT_TILE_DIM + threadIdx.y - FILTER_RADIUS;

    // 2) 協同載入 input tile 進 shared memory;ghost cell 補 0
    __shared__ float N_s[IN_TILE_DIM][IN_TILE_DIM];
    if (row >= 0 && row < height && col >= 0 && col < width) {
        N_s[threadIdx.y][threadIdx.x] = N[row*width + col];
    } else {
        N_s[threadIdx.y][threadIdx.x] = 0.0f;     // ghost cell → 0
    }
    __syncthreads();                              // barrier:確保整塊載齊

    // 3) 只有「內圈」active threads 計算 output
    int tileCol = threadIdx.x - FILTER_RADIUS;
    int tileRow = threadIdx.y - FILTER_RADIUS;
    if (col >= 0 && col < width && row >= 0 && row < height) {        // 邊界
        if (tileCol >= 0 && tileCol < OUT_TILE_DIM &&
            tileRow >= 0 && tileRow < OUT_TILE_DIM) {                 // 停用外圈
            float Pvalue = 0.0f;
            for (int fRow = 0; fRow < 2*FILTER_RADIUS+1; fRow++)
                for (int fCol = 0; fCol < 2*FILTER_RADIUS+1; fCol++)
                    Pvalue += F_c[fRow][fCol]
                            * N_s[tileRow+fRow][tileCol+fCol];        // 全部讀 shared
            P[row*width + col] = Pvalue;
        }
    }
}

外圈停用 (deactivation) 圖解 — 3×3 filter (r=1)、8×8 block/input tile、6×6 output tile:

 block = input tile (8x8),  threadIdx 範圍 0..7
 +---------------------------+
 | .  .  .  .  .  .  .  .   |   '.' = 外圈被停用的 thread (只載入,不計算)
 | .  #  #  #  #  #  #  .   |   '#' = active thread, tx,ty ∈ [1..6]
 | .  #  #  #  #  #  #  .   |
 | .  #  #  #  #  #  #  .   |   active (tx,ty) → output (tx-r, ty-r)
 | .  #  #  #  #  #  #  .   |   patch 左上角 = N_s[ty-r][tx-r]
 | .  #  #  #  #  #  #  .   |
 | .  #  #  #  #  #  #  .   |   thread(1,1) → output(0,0), patch 左上角 N_s[0][0]
 | .  .  .  .  .  .  .  .   |   thread(6,6) → output(5,5), patch 左上角 N_s[5][5]
 +---------------------------+

算術強度分析 (Arithmetic-to-Global-Memory Ratio)

只計算「內部 block」(halo 都不是 ghost cell);large input + small filter 時 ghost cell 影響可忽略。

• 運算量:每個 output 元素做 (2r+1)² 次乘加 = 2(2r+1)² ops → 整塊 OUT_TILE_DIM² · (2r+1)² · 2
• 記憶體:所有 global 存取都移到「載入 input tile」階段,每元素載 4 bytes → 整塊 IN_TILE_DIM² · 4 = (OUT_TILE_DIM+2r)² · 4 bytes

漸近上界 (當 OUT_TILE_DIM ≫ r,分母 ≈ OUT_TILE_DIM):

直覺:原始 kernel 中每個 N 元素被約 (2r+1)² 個 thread 重複載入;tile 無限大時每元素只載一次,故上界即「重複次數 × 2 ops / 4 bytes」。

• 32×32 input tile 是目前 GPU 的上限 (1024 threads/block)。
• tile 越大 → 越接近上界(halo 佔比下降)。
• filter 越大 → 上界越高(每個 N 被更多 thread 用),但實際值與上界落差也越大(halo 變多 → output tile 被擠小)。

方案 B:用 Cache 處理 Halo (Caching for Halo Cells)

核心觀察:halo cells 是鄰居 tile 的內部元素,當本 block 需要它們時,很可能已被鄰居 block 的存取帶進 L2 cache,不必再吃 DRAM 流量。因此只把 tile 內部元素載進 shared memory,halo 直接從 global memory(實則 L2)讀。

→ input tile = output tile = TILE_DIM,block 也同尺寸 (可為 2 的次方)。

#define TILE_DIM 32
__constant__ float F_c[2*FILTER_RADIUS+1][2*FILTER_RADIUS+1];

__global__ void convolution_cached_tiled_2D_const_mem_kernel(float *N, float *P,
                                                             int width, int height) {
    int col = blockIdx.x*TILE_DIM + threadIdx.x;
    int row = blockIdx.y*TILE_DIM + threadIdx.y;

    // 載入:不含 halo → 條件只需檢查一般陣列邊界 (無 ghost cell 之虞)
    __shared__ float N_s[TILE_DIM][TILE_DIM];
    if (row < height && col < width)
        N_s[threadIdx.y][threadIdx.x] = N[row*width + col];
    else
        N_s[threadIdx.y][threadIdx.x] = 0.0f;
    __syncthreads();

    if (col < width && row < height) {
        float Pvalue = 0.0f;
        for (int fRow = 0; fRow < 2*FILTER_RADIUS+1; fRow++) {
            for (int fCol = 0; fCol < 2*FILTER_RADIUS+1; fCol++) {
                if (threadIdx.x-FILTER_RADIUS+fCol >= 0 &&
                    threadIdx.x-FILTER_RADIUS+fCol < TILE_DIM &&
                    threadIdx.y-FILTER_RADIUS+fRow >= 0 &&
                    threadIdx.y-FILTER_RADIUS+fRow < TILE_DIM) {
                    // 內部 → 讀 shared memory
                    Pvalue += F_c[fRow][fCol]
                            * N_s[threadIdx.y-FILTER_RADIUS+fRow]
                                 [threadIdx.x-FILTER_RADIUS+fCol];
                } else {
                    // halo → 檢查是否 ghost cell;非 ghost 才從 global memory 讀
                    if (row-FILTER_RADIUS+fRow >= 0 && row-FILTER_RADIUS+fRow < height &&
                        col-FILTER_RADIUS+fCol >= 0 && col-FILTER_RADIUS+fCol < width) {
                        Pvalue += F_c[fRow][fCol]
                                * N[(row-FILTER_RADIUS+fRow)*width
                                    + col-FILTER_RADIUS+fCol];
                    }
                }
            }
        }
        P[row*width + col] = Pvalue;
    }
}

本章總結 (Chapter Summary)

• Convolution 是核心平行 pattern,也是 stencil(08-Stencil/01-Stencil-Background-and-Basic-Kernel)、grid 力/位能計算(17-Iterative-MRI-Reconstruction/01-MRI-Background-and-Iterative-Reconstruction)、CNN(16-Deep-Learning/02-Convolutional-Neural-Networks-Layers) 的基礎。
• 三步演進:basic kernel (DRAM-bound, 0.25 OP/B) → constant memory 消除 filter 的 DRAM 存取 (0.5 OP/B) → tiling 用 shared memory 降低 N 的 DRAM 存取 (數倍提升,但增加 divergence 與程式複雜度) → L2 cache for halo 簡化程式。
• 算術強度分析是可遷移的重要技能;從中看出小 tile 的限制,在大 filter 與 3D convolution 尤其嚴重。
• 技巧可直接推廣到 1D / 3D,只是索引計算與迴圈巢狀更深。

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

Related Notes

• 07-Convolution/01-Convolution-Fundamentals-and-Basic-Kernel
• 07-Convolution/02-Constant-Memory-and-Caching
• 05-Memory-Architecture-And-Data-Locality/02-Tiling-and-Tiled-Matrix-Multiplication
• 08-Stencil/02-Shared-Memory-Tiling-for-Stencil
• 06-Performance-Considerations/04-Optimization-Checklist-and-Bottlenecks
• 04-Compute-Architecture-And-Scheduling/02-Warps-SIMD-and-Control-Divergence