第3章可扩展的并行执行 - Grid/Block/Thread (Scalable Parallel Execution - Grid/Block/Thread)

一、导读

3.1 本章学习目标

通过本章的学习，读者将能够：

深入理解CUDA编程模型中的三层线程层次结构：Grid、Block、Thread

掌握线程索引计算的方法，能够正确映射线程到数据元素

理解Block作为独立执行单元的特性及其在硬件上的调度机制

掌握Grid和Block维度的配置策略，针对不同问题规模选择合适的参数

理解线程同步原语__syncthreads()的工作原理和使用约束

学会分析kernel的资源需求（寄存器、共享内存）对并发执行的影响

掌握多维Grid和Block的使用方法，适用于多维数据处理

理解CUDA Occupancy的概念及其对性能的影响，能够使用Occupancy Calculator进行优化

3.2 为什么本章如此重要

CUDA的线程层次结构是GPU并行编程的基础抽象。与传统的多线程编程不同，CUDA引入了独特的Grid-Block-Thread三层模型，这一模型既反映了GPU硬件的架构特点，也为程序员提供了灵活的并行表达手段。

正确理解和使用线程层次结构对于编写高效、可扩展的CUDA代码至关重要。线程索引的计算错误是最常见的CUDA bug之一，可能导致数据访问越界或计算结果错误。Block大小的选择不当可能导致资源利用率低下，无法充分发挥GPU的计算能力。同步原语的误用可能引发死锁或竞态条件，导致程序崩溃或结果不可预测。

更重要的是，CUDA的线程层次结构设计体现了可扩展性的核心理念：代码无需修改即可在不同规模的GPU上运行。一个编写良好的CUDA kernel可以在拥有几个SM的移动GPU上正确执行，也可以在拥有上百个SM的数据中心GPU上获得极高的吞吐量。这种可移植性是CUDA成功的关键因素之一。

对于希望深入掌握GPU编程的开发者而言，透彻理解Grid、Block、Thread的概念及其与硬件的映射关系，是迈向高级优化的必经之路。

3.3 核心问题引导

CUDA的线程层次结构是如何组织的？Grid、Block、Thread之间的关系是什么？

如何根据数据规模计算合适的Grid和Block维度？

为什么Block的大小通常是32的倍数（warp大小的倍数）？

Block在SM上是如何调度和执行的？一个Block可以被拆分到多个SM上执行吗？

__syncthreads()的作用范围是什么？它能否用于跨Block同步？

什么是Occupancy？影响Occupancy的因素有哪些？

如何处理数据规模不是Block大小整数倍的情况？

多维Grid和Block适用于哪些应用场景？如何简化多维索引计算？

二、核心概念详解

3.1 CUDA线程层次结构

CUDA采用三层线程层次结构来组织并行执行：Thread（线程）、Block（线程块）、Grid（网格）。这一层次结构既是对并行任务的逻辑划分，也与GPU硬件的执行机制紧密对应。

Thread（线程）：线程是CUDA中最基本的执行单位。每个线程执行相同的kernel代码，但通过内置变量（如threadIdx）区分自己的身份，从而操作不同的数据。线程是最细粒度的并行单元，对应于SIMT模型中的单个执行流。

Block（线程块）：Block是一组线程的集合，这些线程可以相互协作，通过共享内存交换数据，并使用__syncthreads()进行同步。Block是CUDA中最重要的逻辑分组单位，具有以下关键特性：

同一Block内的线程可以同步：Block内的所有线程可以在执行过程中通过__syncthreads()进行屏障同步，确保所有线程都到达同步点后再继续执行。

同一Block内的线程共享共享内存：每个Block拥有自己的共享内存空间，Block内的所有线程都可以访问。这使得Block成为数据协作的基本单位。

Block是原子调度单位：一个Block一旦被调度到某个SM上执行，就不能被迁移到其他SM。Block内的所有线程必须在同一个SM上执行。

Block之间相互独立：不同Block的执行顺序是不确定的，可能并行执行，也可能串行执行。程序员不能假设Block之间的任何执行顺序或同步关系。

Grid（网格）：Grid是最高层次的线程组织单位，由多个Block组成。一次kernel启动会创建一个Grid，Grid中的所有Block执行相同的kernel函数。Grid的主要特点是：

Grid中的所有Block执行相同的kernel代码：这是SIMT模型的体现，所有线程执行相同的指令序列。

Grid的规模可以非常大：Grid可以包含数百万个Block，远超当前GPU可以同时执行的Block数量。GPU硬件会动态调度Block到可用的SM上执行。

Grid定义了问题的总体规模：Grid的维度通常对应于待处理数据的规模，确保每个数据元素都有一个对应的线程进行处理。

层次结构的映射：

Grid (Kernel Launch)
├── Block (0, 0)
│   ├── Thread (0, 0, 0)
│   ├── Thread (1, 0, 0)
│   ├── ...
│   └── Thread (blockDim.x - 1, blockDim.y - 1, blockDim.z - 1)
├── Block (1, 0)
│   ├── Thread (0, 0, 0)
│   └── ...
└── ...

这种层次结构的设计有多个好处：

可扩展性：代码无需关心实际有多少SM可用，GPU会自动调度Block到可用的SM上
局部性：Block内的线程可以共享数据和同步，有利于数据复用
灵活性：可以根据问题规模灵活选择Grid和Block的维度

3.2 线程索引计算

正确计算线程的全局索引是CUDA编程的基本技能。CUDA提供了多个内置变量来帮助线程确定自己的位置：

一维情况：

对于一维的Grid和Block，全局线程索引的计算公式为：

cpp

int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;

其中：

blockIdx.x：当前Block在Grid中的索引
blockDim.x：每个Block中的线程数（x维度）
threadIdx.x：当前线程在Block中的索引

二维情况：

对于二维的Grid和Block，需要将二维索引展平为一维：

cpp

// For 2D block within 2D grid
int globalIdxX = blockIdx.x * blockDim.x + threadIdx.x;
int globalIdxY = blockIdx.y * blockDim.y + threadIdx.y;
int globalIdx = globalIdxY * gridDim.x * blockDim.x + globalIdxX;

或者，如果数据本身是二维的（如图像），可以直接使用二维索引：

cpp

int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;

if (col < width && row < height) {
    // Process pixel at (row, col)
    output[row * width + col] = process(input[row * width + col]);
}

三维情况：

对于三维数据（如体素数据），需要处理三个维度：

cpp

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int z = blockIdx.z * blockDim.z + threadIdx.z;

if (x < dimX && y < dimY && z < dimZ) {
    int idx = z * dimY * dimX + y * dimX + x;
    // Process voxel at (x, y, z)
}

边界检查的重要性：

当数据规模不是Block大小的整数倍时，Grid中的某些线程可能没有对应的数据元素。因此，必须在kernel中加入边界检查：

cpp

int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < N) {
    // Only threads with valid data perform computation
    output[idx] = compute(input[idx]);
}

忽略边界检查可能导致数组越界访问，产生错误结果或程序崩溃。

线性索引与多维索引的选择：

对于一维数据（向量、数组），使用一维索引最简单
对于二维数据（图像、矩阵），使用二维索引更直观
对于三维数据（体素、3D网格），使用三维索引

选择与数据结构匹配的索引方式可以提高代码的可读性和可维护性。

3.3 Block的调度与执行

理解Block在硬件上的调度机制对于性能优化至关重要。

Block到SM的映射：

当kernel启动时，CUDA运行时将Grid中的Block分配给可用的SM执行。调度过程遵循以下原则：

Block是原子单位：一个Block的所有线程必须在同一个SM上执行，不能被拆分到多个SM。

Block可以分配到任意SM：调度器根据SM的资源可用性（寄存器、共享内存、线程槽位）决定将Block分配到哪个SM。

Block的执行顺序不确定：程序员不能依赖Block之间的执行顺序。不同次运行中，Block的执行顺序可能不同。

多个Block可以在同一SM上并发执行：只要SM有足够的资源，就可以同时驻留多个Block。这有助于提高资源利用率和隐藏延迟。

Warp调度：

在SM内部，Block被进一步划分为warp进行调度。每个warp包含32个线程，按照线程ID连续分配。例如，Block中的线程0-31构成第一个warp，32-63构成第二个warp，依此类推。

SM的warp调度器在每个时钟周期选择一个或多个就绪的warp发射指令。当一个warp等待内存访问或其他长延迟操作时，调度器切换到另一个就绪的warp，从而隐藏延迟。

Block的资源限制：

SM可以同时驻留的Block数量受到以下资源的限制：

线程数量限制：每个SM有最大线程数限制（如Ampere为1536或2048）。如果每个Block有256个线程，则SM最多可以驻留1536/256 = 6个Block。

Block数量限制：每个SM有最大Block数限制（如Ampere为16或32）。即使线程数量未达上限，Block数量也不能超过此限制。

寄存器限制：每个SM的寄存器文件总大小有限。如果每个线程使用较多寄存器，则会限制可以驻留的线程数量。

共享内存限制：每个SM的共享内存容量有限。如果每个Block使用较多共享内存，则会限制可以驻留的Block数量。

这些限制共同决定了SM的Occupancy，即活跃warp数量与最大支持warp数量的比值。

3.4 线程同步原语

__syncthreads()是CUDA中最基本的线程同步原语，用于Block内线程的屏障同步。

工作原理：

当线程执行到__syncthreads()时，它会暂停执行，直到Block内的所有其他线程也到达该同步点。一旦所有线程都到达同步点，所有线程继续执行后续代码。

cpp

__global__ void exampleKernel(float* data) {
    __shared__ float shared[256];
    
    int tid = threadIdx.x;
    
    // Phase 1: Load data into shared memory
    shared[tid] = data[tid];
    
    // Synchronize: Ensure all threads have loaded their data
    __syncthreads();
    
    // Phase 2: Process data from shared memory
    // All threads can safely read any element in shared[]
    float result = shared[tid] + shared[(tid + 1) % 256];
    data[tid] = result;
}

使用约束：

仅限Block内同步：__syncthreads()只能用于同一Block内的线程同步，不能跨Block同步。尝试跨Block同步会导致死锁，因为不同Block可能在不同时间执行，甚至在不同SM上执行。

条件同步的危险：如果__syncthreads()放在条件分支内，必须确保Block内所有线程都进入或不进入该分支。否则，部分线程等待同步，而其他线程永远不到达同步点，导致死锁。

cpp

// WRONG: May cause deadlock if not all threads take the same path
if (threadIdx.x < 128) {
    // Do something
    __syncthreads();  // Deadlock! Threads >= 128 never reach here
}

// CORRECT: All threads reach __syncthreads()
if (threadIdx.x < 128) {
    // Do something
}
__syncthreads();  // All threads synchronize here

Volta及以后的独立线程调度：在Volta架构之前，warp内的线程严格锁步执行，__syncthreads()足以保证warp内和warp间的同步。Volta引入独立线程调度（ITS）后，线程可以更独立地推进，对于warp级别的协作，可能需要使用__syncwarp()进行更细粒度的同步。

其他同步原语：

__syncwarp(mask)：warp级别的同步，仅同步掩码指定的线程。比__syncthreads()轻量，适用于warp内协作。

atomicAdd()、atomicMax()等：原子操作，用于线程间的安全累加或比较交换。原子操作可以用于跨Block的间接同步（通过全局内存标志），但效率较低，应谨慎使用。

cudaDeviceSynchronize()：主机端函数，等待设备上所有kernel完成。用于主机和设备之间的同步。

3.5 Grid和Block维度的配置策略

选择合适的Grid和Block维度是CUDA编程中的关键决策，直接影响性能和可扩展性。

Block大小的选择：

常见的Block大小是128、256、512，选择时应考虑以下因素：

Warp对齐：Block大小应该是32（warp大小）的倍数，以避免warp内的线程浪费。例如，Block大小为256意味着每个Block有8个完整的warp。

Occupancy目标：为了最大化资源利用率，Block大小应该使得SM可以驻留足够多的Block。例如，如果SM最多支持1536个线程，Block大小为256时可以有6个Block，Block大小为512时只有3个Block。较高的Block数量有助于更好的负载平衡。

共享内存需求：如果每个Block需要大量共享内存，较小的Block大小可以允许更多Block并发执行。

寄存器使用：如果每个线程使用较多寄存器，可能需要调整Block大小以平衡Occupancy。

问题特性：对于需要大量线程间协作的问题，较大的Block可以提供更多的共享内存和同步机会。对于独立的数据并行任务，较小的Block可能更灵活。

Grid大小的计算：

Grid大小应该覆盖整个数据范围，通常按以下方式计算：

cpp

int blockSize = 256;
int gridSize = (N + blockSize - 1) / blockSize;  // Ceiling division

这里使用向上取整除法，确保即使N不是blockSize的整数倍，也有足够的线程处理所有数据。

多维Grid和Block的应用：

多维配置特别适用于多维数据：

图像处理：使用二维Block和Grid，自然映射到图像的行列
3D体渲染：使用三维Block和Grid，映射到体素的xyz坐标
矩阵运算：使用二维配置，行和列分别对应Block的两个维度

cpp

// 2D configuration for image processing
dim3 blockSize(16, 16);  // 256 threads per block
dim3 gridSize((width + 15) / 16, (height + 15) / 16);

processImage<<<gridSize, blockSize>>>(d_input, d_output, width, height);

动态配置：

对于需要在不同GPU上运行的代码，可以根据设备属性动态配置Block大小：

cpp

cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, deviceId);

int maxThreadsPerBlock = prop.maxThreadsPerBlock;
int maxThreadsPerSM = prop.maxThreadsPerMultiProcessor;

// Choose block size based on device capabilities
int blockSize = 256;
if (maxThreadsPerSM >= 2048) {
    blockSize = 256;  // More blocks for higher occupancy
} else {
    blockSize = 128;  // Smaller blocks for older devices
}

3.6 Occupancy理论与优化

Occupancy是衡量GPU资源利用率的重要指标，定义为活跃warp数量与SM最大支持warp数量的比值。

Occupancy的计算：

Occupancy = Active Warps per SM / Max Warps per SM

例如，如果SM最多支持48个warp，而当前有24个warp活跃，则Occupancy为50%。

影响Occupancy的因素：

每Block线程数：更多的线程每Block意味着更少的Block可以驻留，可能降低Occupancy。

每线程寄存器使用量：更多的寄存器使用限制了可以驻留的线程数量。

每Block共享内存使用量：更多的共享内存使用限制了可以驻留的Block数量。

Occupancy与性能的关系：

高Occupancy通常有利于性能，因为更多的活跃warp可以更好地隐藏内存延迟。然而，高Occupancy并不是性能好的充分条件：

如果kernel是计算受限的（compute-bound），增加Occupancy可能不会带来显著提升
如果每个线程需要做大量工作，较低的Occupancy可能允许每个线程使用更多寄存器和共享内存，从而简化代码逻辑
过高的Occupancy可能导致缓存压力增加，降低缓存命中率

Occupancy优化工具：

CUDA提供了Occupancy Calculator（Excel表格）和Nsight Compute工具，可以帮助分析kernel的资源使用和理论Occupancy。通过这些工具，可以尝试不同的Block大小和资源配置，找到最佳平衡点。

`__launch_bounds__`编译指示：

cpp

__global__ void __launch_bounds__(256, 4) myKernel(float* data) {
    // Kernel code
}

这个编译指示告诉编译器：预期每个Block最多256个线程，每个SM最多4个Block。编译器会根据这些信息优化寄存器分配，避免因寄存器使用过多而导致Occupancy过低。

三、程序员的视角

3.1 在实际项目中的应用

大规模数据处理：在处理GB级别的数据时，合理配置Grid和Block维度至关重要。例如，在日志分析系统中，可以将每条日志记录映射到一个线程，使用较大的Grid覆盖所有记录，使用适中的Block大小（如256）平衡资源利用。

图像处理流水线：对于实时视频处理，通常使用二维Block和Grid，每个线程处理一个像素。Block大小选择16x16或32x8，以适应图像的宽高比。通过共享内存暂存相邻像素，可以实现高效的滤波和卷积操作。

深度学习推理：在自定义的神经网络层实现中，Batch维度通常映射到Grid的一个维度，特征维度映射到Block。通过合理配置，可以确保每个SM上有足够的Block并发执行，最大化吞吐量。

科学模拟：在粒子模拟或流体动力学仿真中，每个粒子或网格点映射到一个线程。由于模拟通常需要多次迭代，Block内的线程可以通过共享内存高效交换边界信息，减少全局内存访问。

3.2 常见陷阱与误区

忘记边界检查：当数据规模不是Block大小的整数倍时，多余的线程可能访问越界内存。始终在kernel开头添加if (idx < N)检查。

错误的索引计算：混淆blockIdx、threadIdx、blockDim的使用，导致索引计算错误。建议在复杂的多维情况下，先写出数学公式，再转换为代码。

条件同步导致死锁：将__syncthreads()放在条件分支内，且不同线程走不同分支，导致部分线程永远等待。确保所有线程要么都执行同步，要么都不执行。

假设Block执行顺序：不同Block的执行顺序是不确定的，不能假设Block 0一定在Block 1之前完成。如果需要Block间的依赖，必须使用多个kernel调用，通过主机端同步。

Block过大导致低Occupancy：选择过大的Block大小（如1024）可能导致每个SM只能驻留很少的Block，降低延迟隐藏能力。通常256是一个较好的默认选择。

忽视共享内存限制：每个Block请求过多的共享内存可能导致SM只能驻留很少的Block。检查设备的共享内存总量，合理分配每Block的使用量。

3.3 与其他层级的联系

与内存层次的联系：Block是共享内存的作用域，理解Block的边界有助于设计有效的数据复用策略。同一Block内的线程可以通过共享内存高效交换数据，而不同Block之间的数据交换需要通过全局内存。

与同步机制的联系：__syncthreads()的作用范围是Block，理解这一点对于设计正确的并行算法至关重要。跨Block的同步需要使用更复杂的机制，如多kernel启动或原子操作。

与性能分析的联系：Nsight Compute等工具可以显示每个kernel的Occupancy、寄存器使用、共享内存使用等指标。理解这些指标与Grid/Block配置的关系，有助于进行有针对性的优化。

四、代码实战

4.1 通用向量操作模板

展示正确处理边界和多维索引的通用模板。

cpp

#include <cuda_runtime.h>
#include <stdio.h>

// Generic 1D vector operation template
template<typename T>
__global__ void vectorOp1D(const T* input, T* output, int N, 
                            void (*op)(const T&, T&)) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Boundary check is critical
    if (idx < N) {
        op(input[idx], output[idx]);
    }
}

// Generic 2D matrix operation template
template<typename T>
__global__ void matrixOp2D(const T* input, T* output, 
                            int width, int height,
                            void (*op)(const T&, T&)) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col < width && row < height) {
        int idx = row * width + col;
        op(input[idx], output[idx]);
    }
}

// Example operation: square each element
__device__ void squareOp(const float& in, float& out) {
    out = in * in;
}

// Host wrapper for 1D vector operation
void vectorSquareGPU(const float* h_input, float* h_output, int N) {
    size_t size = N * sizeof(float);
    
    float *d_input, *d_output;
    cudaMalloc((void**)&d_input, size);
    cudaMalloc((void**)&d_output, size);
    
    cudaMemcpy(d_input, h_input, size, cudaMemcpyHostToDevice);
    
    int blockSize = 256;
    int gridSize = (N + blockSize - 1) / blockSize;
    
    // Use a simple inline kernel instead of function pointer
    // (CUDA doesn't support device function pointers well in all versions)
    auto kernel = []<<<gridSize, blockSize>>>(const float* input, float* output, int N) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx < N) {
            output[idx] = input[idx] * input[idx];
        }
    };
    
    // Actual implementation:
    squareKernel<<<gridSize, blockSize>>>(d_input, d_output, N);
    
    cudaMemcpy(h_output, d_output, size, cudaMemcpyDeviceToHost);
    
    cudaFree(d_input);
    cudaFree(d_output);
}

__global__ void squareKernel(const float* input, float* output, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        output[idx] = input[idx] * input[idx];
    }
}

// Test function
void testVectorOps() {
    int N = 1000000;
    float *h_input = (float*)malloc(N * sizeof(float));
    float *h_output = (float*)malloc(N * sizeof(float));
    
    for (int i = 0; i < N; i++) {
        h_input[i] = (float)i;
    }
    
    vectorSquareGPU(h_input, h_output, N);
    
    // Verify
    bool correct = true;
    for (int i = 0; i < N; i++) {
        float expected = h_input[i] * h_input[i];
        if (fabs(h_output[i] - expected) > 1e-3) {
            correct = false;
            break;
        }
    }
    printf("Vector square test: %s\n", correct ? "PASSED" : "FAILED");
    
    free(h_input);
    free(h_output);
}

int main() {
    testVectorOps();
    return 0;
}

4.2 使用共享内存的Block级归约

展示Block内线程协作和同步的典型模式。

cpp

#include <cuda_runtime.h>
#include <stdio.h>

// Optimized parallel reduction using shared memory
__global__ void reductionSharedMem(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Load input into shared memory
    float val = (idx < N) ? input[idx] : 0.0f;
    shared[tid] = val;
    __syncthreads();
    
    // Perform reduction in shared memory
    // Tree-based reduction pattern
    for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();  // Critical: ensure all updates are complete
    }
    
    // First thread in block writes result
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

// Further optimized: unroll last warp to avoid synchronization
__global__ void reductionUnrolled(const float* input, float* output, int N) {
    extern __shared__ float shared[];
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x * 2 + threadIdx.x;
    
    // Load two elements per thread for better utilization
    float val = 0.0f;
    if (idx < N) val += input[idx];
    if (idx + blockDim.x < N) val += input[idx + blockDim.x];
    
    shared[tid] = val;
    __syncthreads();
    
    // Tree reduction
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            shared[tid] += shared[tid + stride];
        }
        __syncthreads();
    }
    
    // Unroll last warp (no sync needed within warp)
    if (tid < 32) {
        volatile float* vshared = shared;
        vshared[tid] += vshared[tid + 32];
        vshared[tid] += vshared[tid + 16];
        vshared[tid] += vshared[tid + 8];
        vshared[tid] += vshared[tid + 4];
        vshared[tid] += vshared[tid + 2];
        vshared[tid] += vshared[tid + 1];
    }
    
    if (tid == 0) {
        atomicAdd(output, shared[0]);
    }
}

// Host function
float gpuReduction(const float* h_input, int N) {
    size_t size = N * sizeof(float);
    
    float *d_input, *d_output;
    cudaMalloc((void**)&d_input, size);
    cudaMalloc((void**)&d_output, sizeof(float));
    
    cudaMemcpy(d_input, h_input, size, cudaMemcpyHostToDevice);
    cudaMemset(d_output, 0, sizeof(float));
    
    int blockSize = 256;
    int gridSize = (N + 2 * blockSize - 1) / (2 * blockSize);
    int sharedMemSize = blockSize * sizeof(float);
    
    reductionUnrolled<<<gridSize, blockSize, sharedMemSize>>>(
        d_input, d_output, N);
    
    float h_output;
    cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
    
    cudaFree(d_input);
    cudaFree(d_output);
    
    return h_output;
}

int main() {
    int N = 10000000;
    float *h_input = (float*)malloc(N * sizeof(float));
    
    float expected = 0.0f;
    for (int i = 0; i < N; i++) {
        h_input[i] = 1.0f;
        expected += 1.0f;
    }
    
    float result = gpuReduction(h_input, N);
    printf("Reduction result: %.0f (expected: %.0f)\n", result, expected);
    printf("Error: %.6f\n", fabs(result - expected));
    
    free(h_input);
    return 0;
}

4.3 多维数据处理示例

展示二维图像处理的完整示例。

cpp

#include <cuda_runtime.h>
#include <stdio.h>

// 2D Gaussian blur kernel
__global__ void gaussianBlur2D(const unsigned char* input, 
                                unsigned char* output,
                                int width, int height,
                                const float* filter,
                                int filterRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col >= width || row >= height) return;
    
    float sum = 0.0f;
    float weightSum = 0.0f;
    
    // Apply 2D convolution
    for (int fr = -filterRadius; fr <= filterRadius; fr++) {
        for (int fc = -filterRadius; fc <= filterRadius; fc++) {
            int r = row + fr;
            int c = col + fc;
            
            // Boundary handling: clamp
            r = max(0, min(height - 1, r));
            c = max(0, min(width - 1, c));
            
            float weight = filter[(fr + filterRadius) * (2 * filterRadius + 1) + 
                                  (fc + filterRadius)];
            sum += input[r * width + c] * weight;
            weightSum += weight;
        }
    }
    
    output[row * width + col] = (unsigned char)(sum / weightSum);
}

// Separable Gaussian blur (optimized)
__global__ void gaussianBlurHorizontal(const unsigned char* input,
                                        float* tempOutput,
                                        int width, int height,
                                        const float* filter,
                                        int filterRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col >= width || row >= height) return;
    
    float sum = 0.0f;
    
    for (int f = -filterRadius; f <= filterRadius; f++) {
        int c = col + f;
        c = max(0, min(width - 1, c));
        
        sum += input[row * width + c] * filter[f + filterRadius];
    }
    
    tempOutput[row * width + col] = sum;
}

__global__ void gaussianBlurVertical(const float* tempInput,
                                      unsigned char* output,
                                      int width, int height,
                                      const float* filter,
                                      int filterRadius) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (col >= width || row >= height) return;
    
    float sum = 0.0f;
    
    for (int f = -filterRadius; f <= filterRadius; f++) {
        int r = row + f;
        r = max(0, min(height - 1, r));
        
        sum += tempInput[r * width + col] * filter[f + filterRadius];
    }
    
    output[row * width + col] = (unsigned char)sum;
}

// Host function for separable blur
void separableGaussianBlur(const unsigned char* h_input,
                           unsigned char* h_output,
                           int width, int height,
                           int filterRadius) {
    size_t imageSize = width * height;
    
    unsigned char *d_input, *d_output;
    float *d_temp;
    
    cudaMalloc((void**)&d_input, imageSize);
    cudaMalloc((void**)&d_output, imageSize);
    cudaMalloc((void**)&d_temp, imageSize * sizeof(float));
    
    cudaMemcpy(d_input, h_input, imageSize, cudaMemcpyHostToDevice);
    
    // Create Gaussian filter
    int filterSize = 2 * filterRadius + 1;
    float *h_filter = (float*)malloc(filterSize * sizeof(float));
    float sigma = filterRadius / 2.0f;
    float sum = 0.0f;
    
    for (int i = 0; i < filterSize; i++) {
        float x = i - filterRadius;
        h_filter[i] = expf(-(x * x) / (2 * sigma * sigma));
        sum += h_filter[i];
    }
    
    // Normalize
    for (int i = 0; i < filterSize; i++) {
        h_filter[i] /= sum;
    }
    
    float *d_filter;
    cudaMalloc((void**)&d_filter, filterSize * sizeof(float));
    cudaMemcpy(d_filter, h_filter, filterSize * sizeof(float), 
               cudaMemcpyHostToDevice);
    
    // Configure 2D grid and block
    dim3 blockSize(16, 16);
    dim3 gridSize((width + 15) / 16, (height + 15) / 16);
    
    // Horizontal pass
    gaussianBlurHorizontal<<<gridSize, blockSize>>>(
        d_input, d_temp, width, height, d_filter, filterRadius);
    
    // Vertical pass
    gaussianBlurVertical<<<gridSize, blockSize>>>(
        d_temp, d_output, width, height, d_filter, filterRadius);
    
    cudaMemcpy(h_output, d_output, imageSize, cudaMemcpyDeviceToHost);
    
    // Cleanup
    cudaFree(d_input);
    cudaFree(d_output);
    cudaFree(d_temp);
    cudaFree(d_filter);
    free(h_filter);
}

int main() {
    int width = 1920;
    int height = 1080;
    
    unsigned char *h_input = (unsigned char*)malloc(width * height);
    unsigned char *h_output = (unsigned char*)malloc(width * height);
    
    // Initialize with random data
    for (int i = 0; i < width * height; i++) {
        h_input[i] = (unsigned char)(i % 256);
    }
    
    separableGaussianBlur(h_input, h_output, width, height, 3);
    
    printf("Gaussian blur completed for %dx%d image\n", width, height);
    
    free(h_input);
    free(h_output);
    
    return 0;
}

五、深入思考

5.1 历史背景

CUDA的线程层次结构设计受到了早期并行编程模型的影响。OpenMP使用团队（team）和工作共享（work-sharing）的概念，MPI使用进程组和通信子。CUDA的创新在于将硬件执行模型（warp、SM）与编程抽象（Thread、Block、Grid）紧密结合，提供了既高效又易用的并行编程框架。

最初的CUDA架构（Tesla）只支持一维的Block和Grid。随着应用需求的多样化，Fermi架构引入了真正的多维支持，使得处理图像、体积数据等多维数据更加直观。Kepler架构增加了动态并行，允许kernel启动子kernel，进一步丰富了线程层次的表现力。

5.2 前沿方向

细粒度同步：Volta引入的独立线程调度打破了严格的warp锁步执行，使得线程同步更加灵活。未来可能会有更细粒度的同步原语，支持更复杂的协作模式。

异步编程模型：CUDA Graphs允许预先录制kernel启动序列，然后一次性提交，减少了启动开销。这对于小kernel的批量执行特别有用。

多GPU协同：随着多GPU系统的普及，如何在多个GPU之间划分Grid和协调Block执行成为一个重要研究方向。NVLink和GPUDirect等技术正在改善多GPU通信的效率。

5.3 延伸阅读

CUDA C Programming Guide - Thread Hierarchy

NVIDIA CUDA Best Practices Guide - Occupancy

"Programming Massively Parallel Processors" Chapter 3

NVIDIA Developer Blog: "How to Optimize Data Transfers and Kernel Execution"

六、实验笔记

实验1：Block大小对性能的影响

目的：研究不同Block大小对kernel执行时间和Occupancy的影响。

步骤：

编写一个简单的向量加法kernel

分别使用Block大小64、128、256、512、1024进行测试

使用Nsight Compute测量每个配置的Execution Time、Occupancy、Registers per Thread

分析结果，找出最佳Block大小

预期：对于简单的内存带宽受限kernel，中等Block大小（256）通常表现最好。过小（64）可能导致延迟隐藏不足，过大（1024）可能降低Occupancy。

实验2：边界条件处理正确性验证

目的：验证不同边界处理策略的正确性和性能。

步骤：

实现三种边界处理策略：截断（clamp）、环绕（wrap）、补零（zero-pad）

测试数据规模不是Block大小整数倍的情况

验证每种策略的输出正确性

比较性能差异

预期：截断策略最简单，性能最好；环绕和补零需要额外的条件判断，可能稍慢。正确性验证需要仔细设计测试用例。

七、本章小结

本章深入探讨了CUDA的线程层次结构，从Thread、Block到Grid的三层模型构成了GPU并行编程的基础。我们详细分析了每个层次的语义和特性：Thread是基本执行单位，Block是协作和同步的单位，Grid是问题规模的抽象。

线程索引计算是CUDA编程的基本技能，我们介绍了一维、二维、三维情况下的索引计算方法，强调了边界检查的重要性。Block的调度机制决定了它在SM上的执行行为，理解这一机制有助于做出合理的配置决策。

同步原语__syncthreads()是Block内协作的关键，但使用时需要注意避免条件同步导致的死锁。Grid和Block维度的配置需要综合考虑warp对齐、Occupancy、资源限制等因素，没有放之四海而皆准的最优值，需要根据具体应用调优。

通过三个代码示例，我们实践了通用向量操作模板、Block级归约和二维图像处理，展示了线程层次结构在实际应用中的使用方法。最后，我们回顾了线程层次结构的历史演进，展望了细粒度同步和异步编程等前沿方向。

掌握线程层次结构是编写高效CUDA代码的前提。只有在正确理解Thread、Block、Grid的语义和关系的基础上，才能设计出可扩展的并行算法，充分发挥GPU的计算潜力。

八、习题

索引计算：对于一个大小为N=1000000的数组，如果选择Block大小为256，Grid大小应该是多少？最后一个Block中有多少个活跃线程？多少个空闲线程？

多维映射：一个图像尺寸为1920x1080，使用Block大小16x16。计算Grid的维度。线程(100, 200)对应的Block索引和线程索引分别是多少？

同步分析：以下代码是否有问题？如果有，说明问题所在并修正：

cpp

__global__ void kernel(float* data, int N) {
    __shared__ float shared[256];
    int tid = threadIdx.x;
    if (tid < N) {
        shared[tid] = data[tid];
        __syncthreads();
        if (tid > 0) {
            shared[tid] += shared[tid - 1];
        }
    }
}

Occupancy计算：假设SM最多支持1536个线程和16个Block。如果每个Block有256个线程，每个线程使用32个寄存器，SM的寄存器文件有65536个寄存器。计算理论最大Occupancy（以warp为单位，SM最多支持48个warp）。

资源限制分析：一个kernel每个Block使用4KB共享内存，每个线程使用24个寄存器。SM有64KB共享内存和65536个寄存器。计算SM最多可以驻留多少个Block？限制因素是什么？

边界处理：编写一个kernel，对数组进行循环移位（circular shift）。处理当shift量大于0和小于0的情况，确保边界访问正确。

Block独立性：解释为什么CUDA不允许Block之间直接同步。如果需要实现跨Block的依赖，有哪些替代方案？各自的优缺点是什么？

性能调优：一个kernel的Occupancy只有25%，分析可能的原因。列出至少三种提高Occupancy的方法，并说明每种方法可能带来的副作用。

动态配置：编写一个函数，根据设备属性（查询cudaDeviceProp）自动选择最优的Block大小。考虑maxThreadsPerBlock、maxThreadsPerMultiProcessor、warpSize等参数。

综合设计：设计一个并行前缀和（scan）算法，使用多层次Block归约。说明如何划分数据到Block，如何在Block内进行scan，以及如何合并各Block的结果。

第3章 可扩展的并行执行 - Grid/Block/Thread (Scalable Parallel Execution - Grid/Block/Thread)