CUDA三维线程层级：Grid/Block/Thread索引与Numba实现

最早我学 CUDA 时，被一套 Grid、Block、Thread 的三维坐标绕了很久。尤其在 Python 里用 Numba 写 kernel，明明 CPU 上只是两层 for 循环的事，一到 GPU 就变成 blockIdx、blockDim、threadIdx 的排列组合，稍不留神，图像输出直接变成一堆错乱的条纹。这篇是同一个系列里的第三篇，之前讲过基础 kernel 写法、内存访问，这次专门把最核心也最容易翻车的 Grid/Block/Thread 三维层次结构拆开讲清楚。

看完这篇，你能搞明白三件事：第一，blockIdx、blockDim、threadIdx 这三个内置变量到底怎么配合使用；第二，1D、2D、3D 核函数的索引换算公式怎么推；第三，用 Python 写三维体数据内核时，块大小、网格大小怎么定，遇到数据尺寸超过 GPU 限制时又该怎么处理。

1. 先分清 Grid、Block、Thread 各自管的是哪一块

1.1 三层盒子：线程从内到外的归属关系

CUDA 的并行模型，用一句话概括就是：你启动一个 kernel，本质上是在一堆“虚拟线程”上同时执行同一段代码。这堆虚拟线程不是一马平川的大平层，而是按照三层结构组织起来：

• Thread 是最小执行单元，真正跑你写的函数体。
• Block 是一组 Thread 的集合，同一个 Block 内的线程可以通过共享内存协作。
• Grid 是所有 Block 的集合，一次 kernel 启动对应的就是整个 Grid。

启动语法里，kernel[(gridDim,), (blockDim,)]前面的元组定义外层盒子装多少个 Block，后面定义每个 Block 里有多少个 Thread。Python 里最常写的vector_add[(blocks,), (threads,)](a, b, c)，就是这一层关系的直接体现。

对应地，kernel 内部有三个只读变量：

• blockIdx：当前 Block 在整个 Grid 里的编号，有 x、y、z 三个分量。
• blockDim：每个 Block 的尺寸，同样有 x、y、z 三个分量。
• threadIdx：当前 Thread 在当前 Block 里的编号，也是 x、y、z 三轴。

很多人第一次看到这三个变量时，会误以为它们是一组“全局坐标”。其实它们只覆盖了各自的层级：blockIdx 是“小区里的哪栋楼”，threadIdx 是“楼里的哪个房间”，而 blockDim 是“这栋楼每层有多少个房间”。

1.2 为什么 CUDA 要把层次结构设计成三维

一个重要原因是，现实数据不全是线性数组。图像天然是二维的，医学影像、气象网格、物理场模拟天然是三维的。如果一律把所有数据拉成一维再计算，光做边界判断就能把人搞疯。

三维 Grid/Block/Thread 的意义，是让程序员能把“数据的空间结构”直接映射到“线程的组织结构”上：

• 一维数据，比如波形、时间序列，用 x 轴就够；
• 二维图像，宽和高分别对应 x、y；
• 三维体数据，x、y、z 三个轴各司其职，代码结构几乎和 CPU 上的三层循环一一对应。

这也解释了为什么 CUDA 的线程层次结构看起来比普通并行框架“啰嗦”。每个维度都预留了 x、y、z，其实是在给不同场景留表达空间，而不是强迫你所有维度都塞满。

1.3 维度上限要心里有数

三维设计再好，也是有边界的。下面这些上限不是“规范里的数学建议”，而是实际运行时硬性的检查点：

• Block 内线程总数不能超过 1024。也就是说，如果 blockDim 是 (256, 4, 1)，总线程数 1024，可以；如果是 (1024, 2, 1)，总线程数 2048，直接启动失败。
• blockDim 每一维也有上限，常规是 x 最多 1024，y 最多 1024，z 最多 64。但就算 z 是 64，总线程乘积也不能超过 1024，所以实际使用中很少有人真的把三维 block 都塞满。
• Grid 每一维的上限，x 很大，通常超过 20 亿；y 和 z 在老一些的架构上通常是 65535。对绝大多数二维图像处理来说，这个量级够用；但一旦图像高度或体数据深度超过 65535，就不能直接拿 gridDim.y/z 硬顶。

超出上限的解决办法，我后面会在 grid-stride loop 那一节专门讲。这里先记住：三维层次是给表达自由度用的，不是让你无限塞数据。

2. 索引换算：从一维到三维的坐标公式

2.1 最常用的一维换算公式

在 kernel 里，最常见的需求是把当前的线程映射到数组下标上。一维情况下，公式极其简单：

idx = blockIdx.x * blockDim.x + threadIdx.x

含义是：当前 Block 前面有几组 Block，每组 Block 里有 blockDim.x 个线程，当前线程在自己 Block 里的偏移量是 threadIdx.x，两者相加就是全局一维下标。

Numba 提供了一个更省事的写法：cuda.grid(1)，它内部做的就是这件事。所以绝大多数一维 kernel 写成：

@cuda.jit def some_kernel(arr): i = cuda.grid(1) if i < arr.size: arr[i] = arr[i] * 2

cuda.grid(1)返回的就是blockIdx.x * blockDim.x + threadIdx.x。还有一个配套函数cuda.gridsize(1)，它返回整个 Grid 里的线程总数，计算公式是gridDim.x * blockDim.x。后面讲循环遍历时会用到它。

2.2 二维换算和内存布局的对应关系

二维情况下，如果直接用x, y = cuda.grid(2)，那么：

• x等于blockIdx.x * blockDim.x + threadIdx.x
• y等于blockIdx.y * blockDim.y + threadIdx.y

这里最容易踩的坑是：到底谁是“行”、谁是“列”。建议把它固定成一套习惯：x永远对应最内侧的内存轴，也就是“列”，对应 NumPy 数组的最后一个维度；y对应“行”，对应数组的第一个维度。

如果是行主序的二维数组，某一行上相邻元素在内存中是连续的，那么 kernel 里用x作为列下标访问arr[y, x]，就能让相邻线程访问相邻显存地址，这正是 GPU 喜欢的合并访问模式。反过来，如果把x当成行下标拿去做转置一类操作，虽然功能上也能算，但访问效率会一落千丈。

在边界判断上，二维数组的形状必须写清楚：

if x < arr.shape[1] and y < arr.shape[0]: out[y, x] = arr[y, x] * scale

shape[0]是行数，对应y的范围；shape[1]是列数，对应x的范围。这个顺序如果写反，典型症状是输出图像发生 90 度旋转外加索引越界。

2.3 三维换算与体素访问

三维情况进一步扩展：

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

Numba 里用x, y, z = cuda.grid(3)一步到位。边界判断依然要对齐数组形状：

if x < vol.shape[0] and y < vol.shape[1] and z < vol.shape[2]: out[x, y, z] = vol[x, y, z] * factor

假设vol的 shape 是(W, H, D)，那么x是宽方向，y是高方向，z是深度方向。这个约定在 kernel 代码里是一致的，问题不大。

真正麻烦的是当你需要把三维坐标映射回一维连续内存时。NumPy 默认的 C 内存布局是行主序，最后一个维度变化最快，所以三维体数据的一维偏移量是：

idx = (z * H + y) * W + x

也就是先转 z，再转 y，最后转 x。如果 kernel 接收的是一个一维数组，又想按体素坐标访问，你会经常用到这个换算。反过来，如果只用vol[x, y, z]这种高级索引，Numba 会在底层帮你算好，但代价是每次都要多做一次坐标换算，性能敏感的内核里还是建议自己算清楚。

2.4 大数组绕行方案：grid-stride loop

前面提到 Grid 的 y、z 维度有 65535 上限，即便 x 维度很大，很多场景也未必够用。更常见的问题是，数据量本身就比gridDim.x * blockDim.x大，那怎么办？

一种思路是往 Grid 里塞更多 Block，直到覆盖所有数据；另一种更稳的思路是 grid-stride loop，让每个线程处理多个元素：

@cuda.jit def saxpy(x, y, alpha): start = cuda.grid(1) stride = cuda.gridsize(1) for i in range(start, x.size, stride): y[i] = alpha * x[i] + y[i]

这里start是每个线程的起点，stride是整个网格的线程总数。比如网格总共有 1024 个线程，数据长度是 10000，那么线程 0 处理 0、1024、2048…… 直到末尾，线程 1 处理 1、1025、2049……这样每个线程都均匀地“跳过”处理。

这种模式有两个好处：一是网格尺寸可以设成固定大小，不依赖数据长度；二是天然支持动态数据，哪怕数组长度在多次启动之间变化，也只需要调整循环条件。缺点是同一个线程处理多个不相邻元素，如果算法里有邻居依赖，比如图像滤波，就不能轻易套用。

3. Python 实操：把 1D、2D、3D 的 Grid/Block/Thread 跑通

3.1 环境准备：用 Numba 还是 PyCUDA

Python 里写 CUDA 内核，最常见的两条路线是 Numba 和 PyCUDA。

PyCUDA 的直接玩法是把 CUDA C 代码写成字符串，再传给 Python 编译器，灵活是真灵活，但调试起来更像在 Python 里套了一层 C 开发，心智负担偏重。Numba 则是用装饰器@cuda.jit直接标记一个 Python 函数，通常更贴近普通 Python 用户的习惯。我自己日常写算法原型用 Numba，只有需要精细控制某些底层特性时才上手 PyCUDA。

安装环境时，如果你用 conda，建议直接装了编译器再装 Numba：

conda install numba

之后检查 GPU 是否被识别：

from numba import cuda cuda.detect()

能列出可用 GPU 就算成功。如果 Numba 报驱动相关错误，多数是 CUDA 驱动与 Numba 版本不匹配，升级或降级其中一个通常能解决。

3.2 一维向量加法：最快上手的模板

一维向量加法是 CUDA 的入门标准动作。我用这个模板验证环境和索引习惯，测试成本最低。

from numba import cuda import numpy as np @cuda.jit def vector_add(a, b, c): tid = cuda.grid(1) if tid < a.size: c[tid] = a[tid] + b[tid] size = 10000 a = np.arange(size, dtype=np.float32) b = np.arange(size, dtype=np.float32) * 2 c = np.empty_like(a) threads = 256 # 向上取整，保证覆盖所有数据 blocks = (size + threads - 1) // threads vector_add[(blocks,), (threads,)](a, b, c) cuda.synchronize() np.testing.assert_allclose(c, a + b, rtol=1e-5) print("OK")

这里有个小细节：(size + threads - 1) // threads是整数除法的向上取整写法。因为 GPU 线程数通常不是size的整数倍，最后一个 Block 可能只有部分线程在工作，所以 kernel 内部一定要加if tid < a.size做边界检查。

返回值方面，我不急着看一眼结果就完事，而是用np.testing.assert_allclose做自动校验，这样以后改块大小、改网格配置，只要结果比对报错，马上就知道哪一步出错。

3.3 二维矩阵缩放：把索引映射到图像坐标

二维情况，我用矩阵缩放来演示。假设有一个(H, W)的矩阵，每个元素乘一个系数。

@cuda.jit def matrix_scale(mat, out, scale): x, y = cuda.grid(2) if x < mat.shape[1] and y < mat.shape[0]: out[y, x] = mat[y, x] * scale H, W = 512, 768 mat = np.random.randn(H, W).astype(np.float32) out = np.empty_like(mat) threads = (16, 16) # 每个 Block 256 个线程 bx = (W + 15) // 16 by = (H + 15) // 16 matrix_scale[(bx, by), threads](mat, out, 2.5) cuda.synchronize()

这段代码的核心是先确定 Block 的二维尺寸，再推算出 Grid 尺寸。网格的设计逻辑是：让每个线程负责矩阵中的一个元素，Block 像积木一样平铺满整个矩阵。

这里的threads = (16, 16)是一个很常用的二维块配置，总线程数 256，同时宽高方向都能覆盖数据。x 方向 16 个线程对应一行里的 16 个连续列，读取mat[y, x]时，这 16 个线程访问的就是相邻内存地址，合并访问效率很高。

3.4 三维体素处理：真正的三维 Block 怎么用

三维体素处理是标题里的“重头戏”。我以一个简单的体素阈值化为例，输入一个三维数组，输出一个 0/1 掩膜。

这里我约定体数据 shape 是(W, H, D)，也就是宽、高、深度。

@cuda.jit def volume_threshold(vol, out, threshold): x, y, z = cuda.grid(3) if x < vol.shape[0] and y < vol.shape[1] and z < vol.shape[2]: out[x, y, z] = 1.0 if vol[x, y, z] > threshold else 0.0 W, H, D = 128, 128, 128 vol = np.random.randn(W, H, D).astype(np.float32) out = np.zeros((W, H, D), dtype=np.float32) threads = (8, 8, 4) # 8*8*4 = 256 个线程 bx = (W + 7) // 8 by = (H + 7) // 8 bz = (D + 3) // 4 volume_threshold[(bx, by, bz), threads](vol, out, 0.0) cuda.synchronize()

三维 Block 的配置 (8, 8, 4)，乘积正好也是 256 个线程，符合“总线程数不超过 1024”的硬约束。x、y、z 三个方向各分配 8、8、4 个线程，是因为三维数据在深度方向往往没有宽高那么大，把 4 放在 z 轴，能让整个 Block 的形状更匹配数据分布。

如果你的体数据非常大，深度超过 65535，那 gridDim.z 就不够用了。此时建议先把体数据用z * (H * W) + y * W + x压平成一维数组，然后改用一维 grid 和前面说过的 grid-stride loop 来处理，这样既绕开维度上限，又能在一定程度上保持代码简洁。

3.5 网格与块大小的换算速查

这只是一组默认起点。更好的做法是写一个工具函数，根据数据 shape 自动算网格，省得每次手算还容易出错。

4. 线程组织方式对性能的影响

4.1 为什么块大小要取 32 的倍数

GPU 真正的硬件执行单位是 warp，一个 warp 固定包含 32 个线程。调度的时候，32 个线程是一起取指令、一起执行的。所以如果你的 Block 大小不是 32 的倍数，最后一个 warp 里就会有线程闲置，造成浪费。

以块大小 100 为例，它会被”补”成 4 个 warp 的调度宽度，最后一个 warp 只有 4 个线程真正干活，那 28 个位置空转。虽然不算致命，但这种浪费会随着 kernel 执行时间放大。

因此我选块大小就三个选项：128、256、512。这几个值都是 32 的倍数，又没有超过 1024 的总线程上限，调度起来干净利落。

4.2 二维和三维块：块形状不是摆设

对二维数据来说，块尺寸最好是 (16, 16) 或 (32, 8)，而不是 (256, 1)。原因有两层。

第一层是边界覆盖。如果块是 (256, 1)，那么一行 1000 像素的图像，一个块只能覆盖一行里的 256 个像素，为了覆盖整行，网格里的 x 方向就要排 4 个块；如果块是 (16, 16)，一个块就能覆盖 16 行、每行 16 像素，块的形状跟图像的二维结构更匹配。

第二层是共享内存和算法局部性。很多二维算法，比如卷积和矩阵乘法，需要让同一个块里的线程协作读取一块数据到共享内存。块的形状越接近数据瓦片，共享内存的使用就越自然。

三维块的道理也一样。(8, 8, 4) 这种块，其实就是把一个 8×8×4 的小立方体当成一个团队，成员之间可以通过共享内存交换数据，处理体素邻域时非常顺手。

4.3 一个不迷信理论的占用率估算过程

所谓占用率，是指一个 SM（流多处理器）上能同时驻留的线程数和硬件上限的比例。占用率越高，越容易隐藏访存延迟。

假设某款 GPU 的 SM 最大支持 2048 个线程并发驻留：

• 块大小 256，那么一个 SM 最多可同时放 8 个 Block；
• 块大小 512，一个 SM 最多放 4 个 Block；
• 块大小 1024，一个 SM 最多放 2 个 Block。

从占用率数字看，256 和 512 都不错。但块越大，调度灵活性越差。原因很简单：SM 上的资源同时受寄存器数量和共享内存限制，如果一个块的寄存器用量很高，块太大容易让资源早早用完，反而拖低实际驻留块数。

所以，块大小 256 是一个很保险的中间值。不要一上来就追求理论最高占用率，而是先跑通功能，再用官方 profiling 工具查看真实占用率和访存瓶颈，按数据说话。

4.4 合并访问：三维结构下的内存友好姿势

GPU 访存性能和线程的组织方式强相关。最关键的一条规则是：同一 warp 里的 32 个线程，最好在同一时刻访问相邻的显存地址。

放到一维场景，就是tid直接作为数组下标；放到二维场景，就是让x对应数组最内侧的维度，让相邻线程的x相邻；放到三维场景，同理，x对应体数据最后那个维度。

如果反过来，线程的 y 和 x 用反了，比如典型的不合并转置写法：

dst[x, y] = src[y, x]

功能上没错，但线程访问 src 时，同一 warp 里相邻线程的x不在连续地址上，而是每次跳过一整行，访存效率会断崖式下降。这就是为什么很多矩阵转置 kernel 都要借助共享内存做分块转置，本质是为了把跨行访问转成“先合并读，再合并写”的模式。

5. 常见问题与排查实录

5.1 启动配置写错：元组里的逗号

Numba 里启动 kernel 的语法，我见过不少新手被元组坑到。最稳的写法是两对括号：

kernel[(blocks,), (threads,)](args)

而不是偷懒写成kernel[blocks, threads]。尤其当你 block 是二维元组时，写成：

kernel[(bx, by), (16, 16)](args)

这是最清晰的格式。如果去掉最外层的括号，不同 Numba 版本对参数的解析方式不一致，很容易出现“明明代码没改，换个环境就报错”的情况。

5.2 索引顺序颠倒：图像条纹和错乱的元凶

二维图像处理里最常见的症状是输出图像出现横向条纹，或者数据看起来像被 90 度旋转后重新排列。这通常是shape[0]和shape[1]的映射写反了。

排查步骤不复杂。先在 Host 端用 NumPy 模拟同样的索引映射：

y = np.arange(H)[:, None] x = np.arange(W)[None, :] expected = mat[y, x] * scale

然后在 kernel 里也输出同样的逻辑。两边对比，如果 Host 对的、GPU 不对，基本可以锁定是 Kernel 内索引公式的问题。把公式一条条拆开，先验证x，再验证y，不要一次改三行。

5.3 内核没有报错，但结果不对

CUDA kernel 是异步执行的。在 Python 里调用kernel[grid, block](args)后，错误不会立刻从 kernel 调用行抛出来，而是被积压到后续的同步点。

所以自查习惯是：每个 kernel 跑完、需要读取结果前，都显式调用：

cuda.synchronize()

如果这时抛出一个异常，再去看操作日志定位。尤其像越界写这类错误，往往不会马上崩溃，而是在后续某个内存分配或数组拷贝时突然给出非法地址错误。加了同步之后，错误出现的位置会明显提前，排查成本低不少。

5.4 三维索引超过维度上限的处理

数据显示为(100000, 500, 3)这类尺寸时，gridDim.y 或 z 很容易逼近甚至超过 65535 的限制。此时优先考虑把多维数据压平成一维，然后让每个线程通过整数除法还原三维坐标：

idx = cuda.grid(1) if idx < vol.size: z = idx // (H * W) rem = idx % (H * W) y = rem // W x = rem % W

直接用整数除法还原坐标，虽然多出几次除法运算，但换来的是完全不受网格维度上限限制，对绝大多数一次性初始化类任务非常划算。

5.5 浮点结果和 CPU 不一样

GPU 上的浮点运算顺序和 CPU 不一定完全相同，累加操作的顺序差异会导致结果末位不同。不要用==判断测试是否通过，而是要放宽到：

np.testing.assert_allclose(gpu_result, cpu_result, rtol=1e-5, atol=1e-5)

如果误差超过这个范围，才需要怀疑内核逻辑有问题，否则都属于正常的浮点误差。

6. 最后分享一点个人使用习惯

我实际写 Python CUDA 时，并不会每段代码都把二维三维结构用满。大部分场景，比如体素阈值、批量元素级运算，我直接用一维 grid、一维 block，然后在 kernel 里把坐标手工拆成二维或三维，因为这样更容易调试，也更容易套 grid-stride loop。

只有真正需要邻域协作的算法，比如卷积、滤波、矩阵乘法分块，我才会启用二维或三维 Block。因为那种场景下，共享内存的边界天然跟 Block 的边界对上，线程协作沟通成本最低。

块大小方面，我的习惯是一律从 256 线程开始写，跑通功能、确认索引正确以后再试 128 或 512，同时配合性能分析工具看实际数据，不在一开始就陷入块大小微调。三维数据则默认 (8, 8, 4)，除非共享内存特别大或特别小，再调整形状。

如果你正在写属于自己的第一个三维 kernel，我的建议是：先用一个小尺寸数据把映射验证好，再上真实数据；先保证结果正确，再谈性能优化。Grid/Block/Thread 这套层次结构，说穿了就是一套坐标系统，把坐标换算吃透，后面的路会顺很多。