遞迴範例：Quadtree
重點總覽 (Overview)
Quadtree 用「遞迴的 device-side kernel launch」實作：一個 quadrant 對應一個 block，block 若發現點數仍 > 門檻就 build_quadtree_kernel<<<4, ...>>> 啟動 4 個 child block,把空間遞迴細分成四象限。
| 項目 | 內容 | 對應程式碼 |
|---|---|---|
| 映射 | 1 node (quadrant) ↔ 1 block;1 次 launch = 4 blocks (四象限) | nodes[blockIdx.x] |
| 終止條件 | num_points ≤ min 或 depth ≥ max_depth 即 return |
check_num_points_and_depth() |
| 計數 | block 內 threads 用 shared-memory atomicAdd 數四象限的點數 | count_points_in_children() |
| 位移 | 對 4 個計數器做 4-element exclusive scan + 父象限 global offset | scan_for_offsets() |
| 重排 | 依 offset 把點搬到下一個 buffer,同象限的點連續排在一起 | reorder_points() |
| 遞迴 | 由 block 的最後一個 thread 準備 child node 並 launch | prepare_children() + <<<4,...>>> |
| Double buffer | points[0]/points[1] ping-pong;最終結果必須在 buffer 0 |
point_selector |
Important
這是本章唯一的「遞迴 (recursion)」案例:child kernel 與 parent kernel 是同一個 function(build_quadtree_kernel 自我呼叫)。對照兄弟筆記的 Bezier 例子,那裡 parent 與 child 是兩個不同 kernel。
演算法與空間細分 (Quadtree Subdivision: Block-per-Quadrant)
- Quadtree:把 2D 空間遞迴切成 4 個等大象限 (quadrant);每個 quadrant 是一個 tree node,存一段點。
- 規則:若某 quadrant 點數 > 固定 minimum(範例設
min = 2)且 depth 未達max_depth,就再切成 4 個 child node。 - 平行映射:depth d 的每個 active node = 一個 block;一次 child launch 一律是
<<<4, blockDim>>>(分支度固定為 4)。
Depth=0 [ Block 0 : 整個空間, 全部點 ]
│ 切 4 象限, launch <<<4>>>
▼
Depth=1 [00] [01] [02] [03] ← 00,02 只有 2 點 → 不再切
│ │
Depth=2 [010..013] [030..033] ← 010,011,012,013,030,031,032 ≤2 → 停
│
Depth=3 [0330][0331][0332][0333] ← 全部 ≤2, 結束
(陰影 = active block;非 active 的象限不啟動 child grid)
象限判定(以 bounding box 中心 center 為界):
x<cx x>=cx
┌─────────┬─────────┐
y>=cy │ TL=[0] │ TR=[1] │ center = ( (min.x+max.x)/2 ,
├─────────┼─────────┤ (min.y+max.y)/2 )
y<cy │ BL=[2] │ BR=[3] │
└─────────┴─────────┘
Tip
「一個 quadrant = 一個 block」讓 SPMD 程式天然支援不規則的工作量:點多的象限自己再 launch 更多 block,點少的象限直接 return,不必由 host 統一安排。這正是 21-CUDA-Dynamic-Parallelism/01-Dynamic-Parallelism-Fundamentals 所說的「device 動態發掘新工作」。
遞迴 Kernel 主結構 (Recursive Kernel)
__global__ void build_quadtree_kernel(Quadtree_node *nodes,
Points *points, Parameters params) {
__shared__ int smem[8]; // [0..3] 計數器, [4..7] scan 後的 offset
Quadtree_node &node = nodes[blockIdx.x]; // 本 block 負責的 node
node.set_id(node.id() + blockIdx.x);
int num_points = node.num_points();
// 1) 終止判斷 (點數 / depth);停止時把點搬回 buffer 0
bool exit = check_num_points_and_depth(node, points, num_points, params);
if (exit) return;
// 2) 算 bounding box 中心 (四象限分界)
const Bounding_box &bbox = node.bounding_box();
float2 center; bbox.compute_center(center);
int range_begin = node.points_begin();
int range_end = node.points_end();
const Points &in_points = points[ params.point_selector ]; // 輸入 buffer
Points &out_points = points[(params.point_selector+1)%2 ]; // 輸出 buffer
// 3) 計數 → 4) scan offset → 5) 重排
count_points_in_children(in_points, smem, range_begin, range_end, center);
scan_for_offsets(node.points_begin(), smem);
reorder_points(out_points, in_points, smem, range_begin, range_end, center);
// 6) 由最後一個 thread 準備 child node 並遞迴 launch 4 個 block
if (threadIdx.x == blockDim.x - 1) {
Quadtree_node *children = &nodes[params.num_nodes_at_this_level];
prepare_children(children, node, bbox, smem);
build_quadtree_kernel<<<4, blockDim.x, 8*sizeof(int)>>>
(children, points, Parameters(params, true)); // depth+1, 切換 selector
}
}
流程圖 (Fig. 21.8 右側):
Assign block→node
│
num_points>min ── No ──► Exit (必要時把點寫回 buffer 0)
&& depth<max?
│ Yes
Compute center of bounding box
│
Count points in children (count_points_in_children)
│
Scan for offsets (scan_for_offsets)
│
Reorder points (reorder_points)
│
Launch 4 children (<<<4, blockDim>>>)
Warning
只有 block 的最後一個 thread (threadIdx.x == blockDim.x-1) 真正做 launch。若每個 thread 都 launch 會產生 blockDim.x × 4 個多餘 grid,並可能耗盡 pending launch pool(見 21-CUDA-Dynamic-Parallelism/04-Execution-Considerations-and-Summary)。
計數:Shared-Memory Atomics (Counting Points)
block 內所有 thread 協作 grid-stride 掃過 [range_begin, range_end),用 atomicAdd 到 shared memory 累加四象限計數器。
__device__ void count_points_in_children(const Points &in_points, int* smem,
int range_begin, int range_end, float2 center) {
if (threadIdx.x < 4) smem[threadIdx.x] = 0; // 4 個計數器歸零
__syncthreads();
for (int it = range_begin + threadIdx.x; it < range_end; it += blockDim.x) {
float2 p = in_points.get_point(it);
if (p.x < center.x && p.y >= center.y) atomicAdd(&smem[0], 1); // TL
if (p.x >= center.x && p.y >= center.y) atomicAdd(&smem[1], 1); // TR
if (p.x < center.x && p.y < center.y) atomicAdd(&smem[2], 1); // BL
if (p.x >= center.x && p.y < center.y) atomicAdd(&smem[3], 1); // BR
}
__syncthreads();
}
- 計數器只有 4 個且在 shared memory,atomic 競爭遠比寫 global memory 便宜 → 屬於 09-Parallel-Histogram/02-Histogram-Optimizations-Privatization-Coarsening-Aggregation 的 privatization 思路(每 block 私有 histogram)。
- grid-stride 迴圈
it += blockDim.x讓相鄰 thread 讀相鄰點 → memory coalescing。
Scan 求位移 (Four-Element Scan + Global Offset)
把 4 個象限計數轉成「每象限在 buffer 中的起始位置」:對 4 個計數器做 exclusive scan,再加上父 quadrant 的 global offset (node.points_begin())。
__device__ void scan_for_offsets(int node_points_begin, int* smem) {
int* smem2 = &smem[4]; // 輸出 offset 寫到 smem[4..7]
if (threadIdx.x == 0) {
for (int i = 0; i < 4; i++) // 4-element exclusive scan (序列)
smem2[i] = (i == 0) ? 0 : smem2[i-1] + smem[i-1];
for (int i = 0; i < 4; i++) // 加上父象限的 global offset
smem2[i] += node_points_begin;
}
__syncthreads();
}
counts (smem[0..3]) : [ nTL | nTR | nBL | nBR ]
exclusive scan : [ 0 | nTL |nTL+nTR| nTL+nTR+nBL ]
+ node_points_begin : ───────► smem2[0..3] = 各象限在 buffer 的起始 index
Tip
元素只有 4 個,所以單一 thread 序列 scan 就夠了,不必動用 11-Prefix-Sum-Scan/01-Scan-Foundations-and-Kogge-Stone 的平行 scan。這是「演算法選擇要看資料規模」的好例子(見 19-Parallel-Programming-And-Computational-Thinking/02-Algorithm-Selection-and-Tradeoffs)。
跨 Buffer 重排與 Double Buffering (Reorder & Ping-Pong)
用 scan 算出的 offset,把點 scatter 到輸出 buffer:對 smem2[q] 做 atomicAdd 取得目的 index,確保同象限的點連續、不互相覆蓋。
__device__ void reorder_points(Points& out_points, const Points &in_points,
int* smem, int range_begin, int range_end, float2 center) {
int* smem2 = &smem[4]; // 各象限「下一個可寫位置」(由 scan 初始化)
for (int it = range_begin + threadIdx.x; it < range_end; it += blockDim.x) {
int dest; float2 p = in_points.get_point(it);
if (p.x < center.x && p.y >= center.y) dest = atomicAdd(&smem2[0], 1); // TL
if (p.x >= center.x && p.y >= center.y) dest = atomicAdd(&smem2[1], 1); // TR
if (p.x < center.x && p.y < center.y) dest = atomicAdd(&smem2[2], 1); // BL
if (p.x >= center.x && p.y < center.y) dest = atomicAdd(&smem2[3], 1); // BR
out_points.set_point(dest, p); // scatter 到輸出 buffer
}
__syncthreads();
}
Double buffer ping-pong:演算法用 points[0]、points[1] 兩個 buffer 輪流當 in/out;每深一層 point_selector 切換一次。
Depth0: in=Buf0 ──reorder──► out=Buf1 (Fig 21.9 A→B)
Depth1: in=Buf1 ──reorder──► out=Buf0
Depth2: in=Buf0 ──reorder──► out=Buf1 ... 依此交替
最終:結果可能停在 Buf1 → 終止時必須複製回 Buf0
__device__ bool check_num_points_and_depth(Quadtree_node &node,
Points *points, int num_points, Parameters params) {
if (params.depth >= params.max_depth ||
num_points <= params.min_points_per_node) {
// 停止遞迴:若目前點在 buffer 1,必須搬回 buffer 0 (對外結果固定在 Buf0)
if (params.point_selector == 1) {
int it = node.points_begin(), end = node.points_end();
for (it += threadIdx.x; it < end; it += blockDim.x)
points[0].set_point(it, points[1].get_point(it));
}
return true;
}
return false;
}
Warning
為什麼要 swap? 對外 API 約定「結果在 buffer 0」。若某 leaf 停在奇數 depth(資料在 buffer 1),return 前必須把自己負責的那段點複製回 buffer 0,否則結果不完整(Fig 21.9E)。
| 比較 | Scatter (本例) | Gather |
|---|---|---|
| 寫入位置 | 由 atomicAdd(offset) 動態決定 |
每 thread 固定寫自己位置 |
| 需要 atomic | 是(搶象限內的下一格) | 否 |
| 適用 | 輸出長度需動態分群(reorder/bucket) | 一對一映射 |
更廣的取捨見 18-Electrostatic-Potential-Map/01-DCS-Scatter-vs-Gather 與 13-Sorting/01-Sorting-Foundations-and-Parallel-Radix-Sort(radix sort 也用「count→scan→scatter」三步)。
遞迴 Launch 與 Node 配置 (Child Launch & Node Array)
prepare_children() 由最後一個 thread 設定 4 個 child node 的 id、bounding box、點範圍,再做 build_quadtree_kernel<<<4, blockDim>>>。
__device__ void prepare_children(Quadtree_node *children, Quadtree_node &node,
const Bounding_box &bbox, int *smem) {
int child_offset = 4 * node.id(); // 在下一層 node 區段中的位移
children[child_offset+0].set_id(4*node.id()+0); // 4 個 child 的 id (= 陣列索引)
children[child_offset+1].set_id(4*node.id()+1);
children[child_offset+2].set_id(4*node.id()+2);
children[child_offset+3].set_id(4*node.id()+3);
const float2 &p_min = bbox.get_min(), &p_max = bbox.get_max();
float2 center; bbox.compute_center(center);
// 各 child 的 bounding box = 四象限 (TL / TR / BL / BR)
children[child_offset+0].set_bounding_box(p_min.x, center.y, center.x, p_max.y); // TL
children[child_offset+1].set_bounding_box(center.x, center.y, p_max.x, p_max.y); // TR
children[child_offset+2].set_bounding_box(p_min.x, p_min.y, center.x, center.y);// BL
children[child_offset+3].set_bounding_box(center.x, p_min.y, p_max.x, center.y);// BR
// 各 child 的點範圍 = scan 出來的 [offset[q], offset[q+1])
children[child_offset+0].set_range(node.points_begin(), smem[4+0]);
children[child_offset+1].set_range(smem[4+0], smem[4+1]);
children[child_offset+2].set_range(smem[4+1], smem[4+2]);
children[child_offset+3].set_range(smem[4+2], smem[4+3]);
}
- Quadtree 存成一個 node array:不用指標串接,而是用
4*node.id()算出 child 在陣列中的位置(類似把完全四元樹壓平存放)。 params透過Parameters(params, true)傳給 child:depth+1、切換point_selector、更新num_nodes_at_this_level。
遞迴規模量化(以 N 個等距點、min=2、完美細分為例):
| N (點數) | depth 層數(含 root) | child kernel launch 數 |
|---|---|---|
| 公式 | ≈ log₄N 直到象限點數 ≤ min |
Σ 4^i = (4^D − 1)/3 |
| 64 (8×8) | 4 (depth 0,1,2,3) | 1+4+16 = 21 |
Warning
child grid 太小會嚴重浪費 GPU。書中建議「child grid 要多 block,或至少 block 內多 thread」;而硬體 nesting depth 上限 24 層,所以只有相對淺的樹能高效實作(細節見 21-CUDA-Dynamic-Parallelism/04-Execution-Considerations-and-Summary)。
考試/面試重點 (Exam / Test Patterns)
| 場景 / 關鍵字 | 答案 / 技巧 |
|---|---|
| 「quadtree 一個 block 負責什麼?」 | 一個 quadrant (node);一次 launch 固定 4 個 block(四象限) |
| 「block 何時 launch children?」 | num_points > min 且 depth < max_depth 兩者皆成立 |
| 64 等距點的 quadtree 最大深度(含 root) | 4(depth 0/1/2/3;64→16→4→1,4≤?... 4>2 仍切,得 1 點停) |
| 同上,total child kernel launches | 21 = 1+4+16 =(4³−1)/3 |
| 「四象限計數器放哪?用什麼?」 | shared memory + atomicAdd(privatization,避免 global atomic) |
| 「為何需要 scan?」 | 把 4 個計數轉成 buffer 內各象限起始 offset(exclusive scan + parent global offset) |
| 「為何只用 1 個 thread 做 scan?」 | 只有 4 個元素,序列 scan 即可,平行 scan overhead 不划算 |
| 「為何用兩個 buffer?結束要做什麼?」 | ping-pong 重排;結果須在 buffer 0,停在 buffer 1 時要複製回 buffer 0 |
| 「reorder 用 scatter 還 gather?」 | scatter:atomicAdd(offset) 取目的 index 把點搬到輸出 buffer |
| 「哪個 thread 做遞迴 launch?」 | 最後一個 thread threadIdx.x==blockDim.x-1(避免重複 launch) |
| count→scan→reorder 三步像哪個演算法 | 與 radix sort 的 count→scan→scatter 同構 |
Related Notes
- 21-CUDA-Dynamic-Parallelism/01-Dynamic-Parallelism-Fundamentals
- 21-CUDA-Dynamic-Parallelism/02-Bezier-Curves-Example
- 21-CUDA-Dynamic-Parallelism/04-Execution-Considerations-and-Summary
- 11-Prefix-Sum-Scan/01-Scan-Foundations-and-Kogge-Stone
- 13-Sorting/01-Sorting-Foundations-and-Parallel-Radix-Sort
- 09-Parallel-Histogram/02-Histogram-Optimizations-Privatization-Coarsening-Aggregation