第3章：CUDA编程模型基础——线程、块、网格与核函数

CUDA编程模型说白了就是教你怎么把手里的并行任务切成小份，然后扔到GPU上几千个核心去跑。我第一次接触这个概念时也有点懵——线程、块、网格，三层结构到底怎么用？别急，我们一步步拆开看。

3.1 线程、块、网格——GPU的“三级组织架构”

我们先从硬件角度想：GPU里动辄几千个CUDA核心，你要怎么调度它们？线程是最小的执行单元，每个线程跑同一段代码（核函数），但处理的数据不同。就像工厂里每个工人各干各的活。光有线程还不够——几千个线程如果散乱无章，管理起来太麻烦。所以CUDA引入了线程块的概念：你把线程分成小组，比如256个线程一组，组内可以共享数据、同步执行。这叫块。最后，这些块再组成一个网格，整个网格就是一次核函数调用覆盖的全部线程。

我个人习惯把这三级结构想象成部队：网格是整个战役，块是连队，线程就是每个士兵。士兵只能在同一连队里互相支援（共享内存、同步），连队之间则独立作战。

你可能会问：为什么要有块这个中间层？说白了，这是为了在硬件上高效映射。GPU的SM（流多处理器）可以同时执行多个块，但每个SM内的共享内存和寄存器是有限的。块太大可能塞不进一个SM，块太小又浪费资源。我经常在调优时去调整块大小，一个常用的起步值是256个线程每块。

3.2 核函数——在GPU上运行的代码

核函数（kernel）就是你要让GPU执行的代码。它跟普通C函数长得很像，但有几个硬性要求：必须用 __global__ 修饰，返回值必须是 void，不能是类的成员函数。每个线程都执行同样的核函数代码——但通过内置变量 threadIdx、blockIdx 和 blockDim 来区分谁处理哪部分数据。

举个例子，我要把两个数组X和Y逐元素相加：

// 核函数定义
__global__ void vecAdd(float *X, float *Y, float *Z, int n) {
    // 计算全局线程索引
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    // 边界检查
    if (idx < n) {
        Z[idx] = X[idx] + Y[idx];
    }
}

注意这里每个线程只处理一个元素。这叫做“线程到数据的映射”。我在项目中见过有人忘了做边界检查，结果数组大小不是块大小的整数倍时，线程跑越界，读到了垃圾数据，整加结果偏掉——嗯，排查了好久才找到这个坑。

3.3 <<<>>> 语法——启动核函数的“开关”

在主机端（CPU）调用核函数时，要用一个特殊的语法：

vecAdd<<>>(d_X, d_Y, d_Z, N);

<<>> 就是告诉CUDA运行时：你开一个网格，网格里有 gridDim 个块，每个块里有 blockDim 个线程。这两个可以是整数（一维），也可以是 dim3 类型（二维、三维）。

比如：

// 一维：256个线程一维，共100个块
vecAdd<<<100, 256>>>(...);

// 二维：16x16的线程块，8x8个块
dim3 blockDim(16, 16);
dim3 gridDim(8, 8);
kernel<<>>(...);

3.4 内置变量——线程如何“认路”

每个线程在核函数内部可以访问一些只读变量，用来计算自己的唯一ID：

变量	含义	示例值（一维）
threadIdx.x	线程在本块内的x索引	0 ~ (blockDim.x - 1)
blockIdx.x	本块在网格内的x索引	0 ~ (gridDim.x - 1)
blockDim.x	每个块在x维度的线程数	由 <<<>>> 指定
gridDim.x	网格在x维度的块数	由 <<<>>> 指定

计算全局线程ID的经典公式：

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

如果是二维块和二维网格：

int globalIdx_x = threadIdx.x + blockIdx.x * blockDim.x;
int globalIdx_y = threadIdx.y + blockIdx.y * blockDim.y;
int globalIdx = globalIdx_y * (gridDim.x * blockDim.x) + globalIdx_x;

我最开始总弄混 blockIdx 和 threadIdx 的维度顺序。有个简单记忆法：先想块内的线程（threadIdx），再想网格里的块（blockIdx）。公式里就是 threadIdx + blockIdx * blockDim，跟矩阵索引很像。

3.5 完整示例——向量加法一步到位

下面是一个完整的向量加法代码框架，包括主机端内存分配、数据拷贝、核函数启动和结果收回：

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

// 核函数：Y = X + Y（就地加法）
__global__ void vecAdd(float *X, float *Y, int n) {
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < n) {
        Y[idx] = X[idx] + Y[idx];
    }
}

int main() {
    int n = 1 << 20;        // 约100万个元素
    size_t bytes = n * sizeof(float);

    // 分配主机内存
    float *h_X = (float*)malloc(bytes);
    float *h_Y = (float*)malloc(bytes);
    // 初始化数据（略）

    // 分配设备内存
    float *d_X, *d_Y;
    cudaMalloc(&d_X, bytes);
    cudaMalloc(&d_Y, bytes);

    // 拷贝到设备
    cudaMemcpy(d_X, h_X, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_Y, h_Y, bytes, cudaMemcpyHostToDevice);

    // 定义网格和块
    int threads_per_block = 256;
    int blocks_per_grid = (n + threads_per_block - 1) / threads_per_block;

    // 启动核函数
    vecAdd<<<blocks_per_grid, threads_per_block>>>(d_X, d_Y, n);

    // 同步并检查错误
    cudaDeviceSynchronize();

    // 拷回结果
    cudaMemcpy(h_Y, d_Y, bytes, cudaMemcpyDeviceToHost);

    // 清理
    cudaFree(d_X); cudaFree(d_Y);
    free(h_X); free(h_Y);
    return 0;
}

这个 blocks_per_grid 的计算方式是经典写法：向上取整 + 边界检查。为什么不用 n / threads_per_block？万一 n 不是256的整数倍，少开的块就会漏掉尾部元素，所以一定要 (n + threads - 1) / threads。而且核函数里要有 if (idx < n) 兜底。

3.6 两个容易忽视的细节

1. 网格和块的最大尺寸：每个维度的 blockDim 最大1024，网格 gridDim 最大 2^31-1（实际上受显卡限制）。块内线程总数不能超过1024。如果你设了 blockDim(32, 32)，那就是1024个线程，刚刚好。
2. 同步问题：不同块之间的线程无法直接同步，只能通过核函数结束来隐含同步。块内的线程可以用 __syncthreads() 同步，但必须保证所有线程都执行到这条语句，否则会死锁。我有一次在条件分支里用 __syncthreads() 没处理好，程序直接卡住，从此养成了“分支内慎用同步”的习惯。