範例：Bezier Curves 與動態工作量

重點總覽 (Overview)

本節以 adaptive subdivision of spline curves（Bezier 曲線自適應細分）示範「每個工作單位的工作量會動態變化」時，Dynamic Parallelism 如何同時 (1) 抽取更多平行度、(2) 消除 control divergence、(3) 用 device-side cudaMalloc 精準配置記憶體。

Bezier 曲線數學背景 (Bezier Curve Math)

Bezier 曲線由一組 control points P0…Pn 定義，n 為階數（order）。首尾控制點一定在曲線上，中間控制點通常不在曲線上。

• 本範例使用 quadratic（3 個控制點）。
• curvature（曲率）越大 → 需要越多 tessellation points 才能畫得平滑，這就是每組控制點工作量不同的來源。

tessellation 取樣：t ∈ [0,1] 等分成 nTessPoints 個點
u = idx / (nTessPoints − 1)        // idx = 0 .. nTessPoints-1
B3u[0]=(1-u)²   B3u[1]=2(1-u)u   B3u[2]=u²
position = B3u[0]*P0 + B3u[1]*P1 + B3u[2]*P2

無 Dynamic Parallelism 的版本 (Without Dynamic Parallelism)

設計：一個 block 處理一組控制點，block 先算 curvature 推得點數，再用 for-loop（步進 blockDim.x）讓 block 內所有 thread 協作填滿 vertex 陣列。

#define MAX_TESS_POINTS 32
struct BezierLine {
    float2 CP[3];                       // 3 control points (quadratic)
    float2 vertexPos[MAX_TESS_POINTS];  // 固定配最壞情況大小！
    int    nVertices;
};

__global__ void computeBezierLines(BezierLine *bLines, int nLines) {
    int bidx = blockIdx.x;                       // <-- block 為單位
    if (bidx < nLines) {
        float curvature = computeCurvature(bLines);
        int nTessPoints = min(max((int)(curvature*CURVE_TESS_FACTOR),4),
                              MAX_TESS_POINTS);   // 工作量 ∝ curvature
        bLines[bidx].nVertices = nTessPoints;
        // for-loop：迭代數因 block 而異 -> 工作量不均
        for (int inc = 0; inc < nTessPoints; inc += blockDim.x) {
            int idx = inc + threadIdx.x;
            if (idx < nTessPoints) {
                float u = (float)idx/(float)(nTessPoints-1);
                float omu = 1.0f - u, B3u[3] = {omu*omu, 2.0f*u*omu, u*u};
                float2 position = {0,0};
                for (int i=0;i<3;i++) position = position + B3u[i]*bLines[bidx].CP[i];
                bLines[bidx].vertexPos[idx] = position;
            }
        }
    }
}

兩個問題：

block0 (curvature 小, 4 pts)   : [#-------]  迭代 1 次後閒置
block1 (curvature 大, 32 pts)  : [########]  迭代 32 次
block2 (curvature 小, 6 pts)   : [##------]  迭代 2 次後閒置
   ^ block 間工作量差異 -> SM 利用率低
   ^ 同一 warp 內 thread 迭代數不同 -> control divergence
   ^ vertexPos 每條都吃 MAX_TESS_POINTS 空間 -> 記憶體浪費

有 Dynamic Parallelism 的版本 (With Dynamic Parallelism)

把原 kernel 拆成 parent（發現工作量） 與 child（執行計算）。

關鍵差異 1：thread 為單位 vs block 為單位

parent 每組控制點的工作（算 curvature + launch）很小，所以用一條 thread做就好（原版用一個 block）。因此 host 端 launch 時把 grid 大小除以 BLOCK_DIM：

// host 端
computeBezierLines_parent<<<ceil((float)N_LINES/BLOCK_DIM), BLOCK_DIM>>>(bLines_d, N_LINES);

關鍵差異 2：device-side cudaMalloc 精準配置

struct BezierLine {
    float2  CP[3];
    float2 *vertexPos;   // <-- 改成指標，動態配置
    int     nVertices;
};

__global__ void computeBezierLines_parent(BezierLine *bLines, int nLines) {
    int lidx = threadIdx.x + blockDim.x*blockIdx.x;   // <-- thread 為單位 (vs Fig.21.6 的 blockIdx.x)
    if (lidx < nLines) {
        float curvature = computeCurvature(bLines);
        bLines[lidx].nVertices = min(max((int)(curvature*CURVE_TESS_FACTOR),4),
                                     MAX_TESS_POINTS);
        // 在 device 端依實際點數精準配置記憶體
        cudaMalloc((void**)&bLines[lidx].vertexPos,
                   bLines[lidx].nVertices*sizeof(float2));
        // launch 一個大小可變的 child grid（依 nVertices）
        computeBezierLines_child<<<ceil((float)bLines[lidx].nVertices/32.0f), 32>>>
            (lidx, bLines, bLines[lidx].nVertices);
    }
}

__global__ void computeBezierLines_child(int lidx, BezierLine *bLines, int nTessPoints) {
    int idx = threadIdx.x + blockDim.x*blockIdx.x;   // 每 thread 只算 1 個 vertex
    if (idx < nTessPoints) {
        float u = (float)idx/(float)(nTessPoints-1);
        float omu = 1.0f - u, B3u[3] = {omu*omu, 2.0f*u*omu, u*u};
        float2 position = {0,0};
        for (int i=0;i<3;i++) position = position + B3u[i]*bLines[lidx].CP[i];
        bLines[lidx].vertexPos[idx] = position;       // 無迴圈、無 divergence
    }
}

thread → child grid 映射

parent grid (N_LINES 條 thread, 每 thread = 一組控制點)
 t0 ─launch→ child grid A: [#-#-]     (4 pts  -> 1 block)
 t1 ─launch→ child grid B: [########################] (32 pts -> 1 block,32thr)
 t2 ─launch→ child grid C: [#-#-#-]   (6 pts  -> 1 block)
 ...
 每條 child grid 大小不同 = variable child grid launches
 child grid 內「1 thread = 1 vertex」-> 完美 load balance, 0 divergence

Device 端記憶體釋放 (Device-Side cudaFree)

device kernel 用 cudaMalloc 配置的記憶體，必須由 device kernel cudaFree 釋放（host 不能直接 free 它）。因此另寫一個 kernel 由 host 呼叫：

__global__ void freeVertexMem(BezierLine *bLines, int nLines) {
    int lidx = threadIdx.x + blockDim.x*blockIdx.x;
    if (lidx < nLines)
        cudaFree(bLines[lidx].vertexPos);   // 平行釋放每條線的 vertex 記憶體
}
// host 端
freeVertexMem<<<ceil((float)N_LINES/BLOCK_DIM), BLOCK_DIM>>>(bLines_d, N_LINES);

工作量、Pool 與 Stream 的關聯 (Workload, Pool & Streams)

此範例是後續「重要執行考量」的 running example，這裡只點出與工作量直接相關的兩個陷阱（細節見 21-CUDA-Dynamic-Parallelism/04-Execution-Considerations-and-Summary）：

• Pending launch pool：parent 每條 thread 各 launch 一個 child grid → 共 launch N_LINES 個 child grid。預設 pool 容量 2048；若 N_LINES 超過會落入 virtualized pool，效能可能慢一個數量級。對策：

  cudaDeviceSetLimit(cudaLimitDevRuntimePendingLaunchCount, N_LINES);
  

• Named streams：stream scope 為 block 私有；不指定 stream 時，同一 block 內所有 launch 都用 NULL stream → 被序列化（即使來自不同 thread）。要讓同 block 各 thread 的 child grid 並行，須每 thread 建立 non-blocking named stream：

  cudaStream_t stream;
  cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);
  computeBezierLines_child<<<ceil((float)bLines[lidx].nVertices/32.0f), 32, 0, stream>>>
      (lidx, bLines, bLines[lidx].nVertices);
  cudaStreamDestroy(stream);
  

關鍵公式與計數 (Key Formulas & Counts)

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

Related Notes

• 21-CUDA-Dynamic-Parallelism/01-Dynamic-Parallelism-Fundamentals
• 21-CUDA-Dynamic-Parallelism/03-Recursive-Quadtree-Example
• 21-CUDA-Dynamic-Parallelism/04-Execution-Considerations-and-Summary
• 04-Compute-Architecture-And-Scheduling/02-Warps-SIMD-and-Control-Divergence
• 15-Graph-Traversal/02-Vertex-Centric-and-Edge-Centric-BFS