GPU 编程实战——Kernel 优化入门

数据智能老司机
136 阅读13分钟

选择线程块（Block）和网格（Grid）大小

在开发更高级的 GPU kernel 时，选择合适的线程块大小和网格大小是影响性能的关键因素。虽然 CUDA 的海量并行能力能加速大规模数据处理，但硬件约束——如每块可用的共享内存、寄存器使用量，以及每个多处理器（SM）可并行运行的线程总数——意味着并非所有的块/网格配置都能带来最佳吞吐率。不合理的设置可能导致 GPU 资源利用不足、产生内存瓶颈，甚至因超出硬件限制而无法启动 kernel。

最佳的 kernel 往往在保持足够大块以让所有 SM 都保持满负荷运行，与避免过大块耗尽资源之间找到平衡。我们的目标是最大化占用率——即在任意时刻实际被用于计算的 GPU 资源比例。

CUDA 将线程组织成两级层次：

  • 线程块（Thread block） ：一组可以通过共享内存和 __syncthreads() 同步协作的线程，所有线程块内的线程都运行在同一个 SM 上。
  • 网格（Grid） ：由 kernel 启动时所有线程块的全体构成，覆盖整个数据集。

线程块大小（每块线程数）和网格大小（块数）的选择决定了工作如何切分并映射到硬件。一般做法是使用 32 的整数倍线程数（warp 大小），常见设置为 128、256 或 512 线程/块。但理想值取决于 kernel 所用的寄存器数量、共享内存需求以及待处理数据量。

确定最优配置

最优块/网格大小取决于三个主要因素：

  1. 问题规模（待处理元素或任务总数）
  2. 硬件限制（每块最大线程数、可用共享内存、寄存器总数）
  3. Kernel 资源使用（每块或每线程的共享内存和寄存器消耗）

我们可以借助 NVIDIA 提供的占用率计算器（Occupancy Calculator）来估算：给定 kernel 的资源需求，该工具会模拟每个 SM 上可并行运行多少个线程块，从而指导我们选择合适的配置。

示例：探索块大小与网格大小

假设我们要为一个大规模一维数组启动一个相加 kernel：

import numpy as np
import cupy as cp

N = 10_000_000  # 元素总数
a = cp.random.rand(N).astype(cp.float32)
b = cp.random.rand(N).astype(cp.float32)
c = cp.empty_like(a)

我们要计算 c[i] = a[i] + b[i]，需选择块和网格维度。

选择块大小并启动 kernel

首先，尝试 256 或 512 线程/块（均为 warp 大小 32 的倍数），同时确保不超过 GPU 最大线程数（常为 1024）：

threads_per_block = 256
blocks_per_grid   = (N + threads_per_block - 1) // threads_per_block

kernel_code = """
extern "C" __global__
void vector_add(const float* a, const float* b, float* c, int n) {
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}
"""
mod        = cp.RawModule(code=kernel_code)
vector_add = mod.get_function("vector_add")

vector_add(
    (blocks_per_grid,), (threads_per_block,),
    (a, b, c, N)
)

评估占用率

要使用 NVIDIA 的占用率计算器，你需要知道：

  • 线程块大小（threads per block）
  • 每块共享内存使用量
  • 每线程寄存器使用量（可用 nvcc --ptxas-options=-v 获取）

对于共享内存和寄存器占用较小的简单 kernel，即使是大块大小也能实现高占用率；若资源需求增加，则占用率会下降。

性能对比与实验

我们可以测量不同块大小下的执行时间，观察吞吐率变化：

import time

for block_size in [64, 128, 256, 512, 1024]:
    grid_size = (N + block_size - 1) // block_size
    c_test = cp.empty_like(a)

    start = time.time()
    vector_add(
        (grid_size,), (block_size,),
        (a, b, c_test, N)
    )
    cp.cuda.Stream.null.synchronize()
    end = time.time()

    print(f"块大小 {block_size}，耗时：{end - start:.5f} 秒")

通常最佳性能出现在 128–512 线程/块之间，但也会因硬件和 kernel 复杂度而异。

通过实验与占用率计算器工具相结合，我们能找到最大化吞吐率的块/网格配置，实现高效、可扩展的并行执行。这也是优化 GPU kernel 性能的核心步骤。

共享内存数据重用

缓冲原理

在 GPU 编程中，全局内存（global memory）访问常常成为吞吐量的瓶颈。每当线程从全局内存读取数据，都要付出显著的延迟成本，往往远超实际计算所需时间。如果同一线程块内的多个线程需访问相同数据（如相邻线程操作重叠区域，或线程块处理矩阵/图像子块），反复的全局内存读取会迅速累积开销。为此，CUDA 提供了共享内存——每个线程块可管理的一小块超高速本地内存。我们只需将所需数据从全局内存加载到共享内存一次，线程块内的所有线程即可反复访问，减轻全局内存流量、提升 kernel 吞吐。

这种缓冲方式非常适合“同一数据被多次重用”的场景，如模板计算（stencil）、矩阵分块、卷积等操作。

典型示例：1D 模板计算（Stencil）

在一维模板计算中，每个输出元素依赖自己及其邻域的值。如果让每个线程分别从全局内存读取这些邻值，将导致大量重复访问。使用共享内存后，每个线程块协作加载包含中心元素与“halo”邻域在内的一段数据，只需一次全局加载，之后各线程在共享内存中快速读取。

准备数据

import numpy as np
import pycuda.autoinit
import pycuda.driver as drv
import pycuda.gpuarray as gpuarray
from pycuda.compiler import SourceModule

N      = 1024 * 1024
radius = 3  # 左右各 3 个邻域
input_host  = np.random.rand(N).astype(np.float32)
output_host = np.zeros_like(input_host)

input_gpu  = gpuarray.to_gpu(input_host)
output_gpu = gpuarray.empty_like(input_gpu)

编写共享内存 Kernel

kernel_code = f"""
#define RADIUS {radius}
#define BLOCK_SIZE 256

__global__ void stencil_1d(const float *input, float *output, int N)
{{
    extern __shared__ float smem[BLOCK_SIZE + 2 * RADIUS];
    int tid        = threadIdx.x;
    int global_idx = blockIdx.x * BLOCK_SIZE + tid;
    int smem_idx   = tid + RADIUS;

    // 加载中心数据
    if (global_idx < N)
        smem[smem_idx] = input[global_idx];

    // 加载左侧 halo
    if (tid < RADIUS && global_idx >= RADIUS)
        smem[smem_idx - RADIUS] = input[global_idx - RADIUS];

    // 加载右侧 halo
    if (tid >= BLOCK_SIZE - RADIUS && (global_idx + RADIUS) < N)
        smem[smem_idx + RADIUS] = input[global_idx + RADIUS];

    __syncthreads();

    // 计算输出（跳过边界）
    if (global_idx >= RADIUS && global_idx < N - RADIUS) {{
        float sum = 0.0f;
        for (int j = -RADIUS; j <= RADIUS; ++j)
            sum += smem[smem_idx + j];
        output[global_idx] = sum / (2 * RADIUS + 1);
    }}
}}
"""
mod         = SourceModule(kernel_code)
stencil_1d  = mod.get_function("stencil_1d")

启动 Kernel

block_size = 256
grid_size  = (N + block_size - 1) // block_size
shared_mem = (block_size + 2 * radius) * np.dtype(np.float32).itemsize

stencil_1d(
    input_gpu, output_gpu, np.int32(N),
    block=(block_size, 1, 1), grid=(grid_size, 1),
    shared=shared_mem
)

验证输出

# CPU 参考实现
for i in range(radius, N - radius):
    output_host[i] = np.mean(input_host[i - radius : i + radius + 1])

# 拷回 GPU 结果
output_gpu_host = output_gpu.get()

print("结果匹配吗？",
      np.allclose(output_host[radius:N-radius],
                  output_gpu_host[radius:N-radius],
                  atol=1e-6))

若不使用共享内存，每个线程及其邻域都要从全局内存读取多次，无疑会产生大量额外流量。引入共享内存后，每个值在每个线程块中仅从全局内存加载一次，其后由多个线程重用。随着块大小和重用率提升，这种优化带来的性能增益非常可观，常常将内存绑定（memory‐bound）kernel 转变为计算绑定（compute‐bound）kernel。

减少全局内存访问

带宽利用不当

当我们的 kernel 越来越复杂时，会发现并非所有的内存访问都是等价的。当同一 warp 中的多个线程访问全局内存时，GPU 硬件会尝试将这些请求“合并”（coalesce）为尽可能少的内存事务。如果线程以规则、对齐、连续的方式访问数据，硬件就能一次性读取所需数据，充分利用带宽；若访问模式混乱无序，可能会导致每个线程都发起单独的内存事务，称为非合并访问（uncoalesced access），会严重浪费带宽，拖慢 kernel 性能。GPU 编程中最常见且最昂贵的低效之一，就是由于非合并的加载和存储造成带宽浪费。为了最大化利用硬件，我们应尽量确保同一 warp（通常 32 线程）内的线程访问连续的内存地址。

实现的方法很简单：在 warp 内将所有加载/存储请求对齐。当线程索引连续时，它们访问的地址也应连续，让内存控制器能将这些请求打包，单次传输一个对齐块，达到最大的吞吐——每个 warp 只需一次内存事务，而不是每个线程一次。我们只需合理安排数据结构和访问模式，确保线程编号相邻的线程访问相邻地址即可。

示例：合并访问 vs 非合并访问

下面通过一个简单的数组复制示例来对比两种访问模式，并展示如何调整访问以实现合并。

场景：在设备上将一个大一维数组复制到另一个数组

import numpy as np
import pycuda.autoinit
import pycuda.driver as drv
import pycuda.gpuarray as gpuarray
from pycuda.compiler import SourceModule

N = 1024 * 1024
input_host = np.random.rand(N).astype(np.float32)
input_gpu  = gpuarray.to_gpu(input_host)
output_gpu = gpuarray.empty_like(input_gpu)

非合并访问 Kernel

每个线程按固定步长（stride）访问，导致访问地址相距较远，无法合并：

kernel_uncoalesced = """
__global__ void copy_uncoalesced(const float *input, float *output, int N, int stride)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N/stride)
        output[idx] = input[idx * stride];
}
"""
mod_uncoalesced    = SourceModule(kernel_uncoalesced)
copy_uncoalesced  = mod_uncoalesced.get_function("copy_uncoalesced")

合并访问 Kernel

每个线程访问连续元素，满足合并访问条件：

kernel_coalesced = """
__global__ void copy_coalesced(const float *input, float *output, int N)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N)
        output[idx] = input[idx];
}
"""
mod_coalesced   = SourceModule(kernel_coalesced)
copy_coalesced  = mod_coalesced.get_function("copy_coalesced")

性能测量

使用相同的块和网格配置，比较两者耗时：

block_size = 256
grid_size  = (N + block_size - 1) // block_size
import time

# 非合并访问
stride = 32
output_uncoalesced = gpuarray.empty_like(input_gpu)
start = time.time()
copy_uncoalesced(
    input_gpu, output_uncoalesced, np.int32(N), np.int32(stride),
    block=(block_size, 1, 1), grid=(grid_size // stride, 1)
)
drv.Context.synchronize()
end = time.time()
print(f"非合并访问耗时：{end - start:.5f} 秒")

# 合并访问
output_coalesced = gpuarray.empty_like(input_gpu)
start = time.time()
copy_coalesced(
    input_gpu, output_coalesced, np.int32(N),
    block=(block_size, 1, 1), grid=(grid_size, 1)
)
drv.Context.synchronize()
end = time.time()
print(f"合并访问耗时：{end - start:.5f} 秒")

通常合并访问版本要快好几倍，因为硬件可以高效地合并内存请求。

数据对齐与结构

在更复杂的场景（如二维、三维数组或结构体数组）中，也应遵循相同原则：保证同一 warp 内相邻线程访问相邻内存地址。通常**结构体数组（SoA）比数组结构（AoS）**更有利于合并访问，因为 SoA 将相同字段存放在连续内存中，而 AoS 会将不同字段交错存储，破坏地址连续性。

通过这种简单调整，让数据和访问模式保持高度规律，我们就能彻底激发全局内存的带宽潜力，将性能低效的 kernel 转化为高效例程，并随着数据规模和计算需求的增长，实现良好的可扩展性。

使用 CuPy Profiler 测量占用率

在优化 GPU kernel 时，我们的主要目标是让硬件保持尽可能高的利用率。希望所有 SM 都被活跃线程饱满占据，最大化寄存器和内存资源的使用。占用率（occupancy）反映了实际运行线程数与理论最大线程数的比值。高占用率通常意味着性能更好；占用率过低会导致大量资源闲置，而占用率过高则可能耗尽寄存器或共享内存，引发寄存器溢出或降低速度。因此，理解并调节占用率对于编写高效 kernel 至关重要。

CuPy Profiler 概览

CuPy 提供了方便的分析 API（在 cupyx.profiler 模块中），可帮助我们在 Python 环境中直接收集 kernel 执行统计数据——包括占用率、每线程寄存器数、内存使用、活跃 warp 数等。借助它，我们无需离开熟悉的工作流，就能快速定位瓶颈并迭代优化。

收集 Kernel 统计信息

函数级别基准

使用 benchmark 可以多次运行某个函数，统计时间和吞吐率：

import cupy as cp
from cupyx.profiler import benchmark

size = 10_000_000
a = cp.random.rand(size, dtype=cp.float32)
b = cp.random.rand(size, dtype=cp.float32)

def vector_add():
    cp.add(a, b, out=a)

# 基准测试 vector_add，重复 10 次
result = benchmark(vector_add, n_repeat=10)
print(result)

输出包含运行时长、吞吐量等指标；在支持的环境下，还能扩展输出内核级别统计。

时间范围标记

若在 Nsight Systems 或 Nsight Compute 中运行，可用 time_range 标记分析区间：

from cupyx.profiler import time_range

with time_range('Vector Add Profile', color_id=0):
    cp.add(a, b, out=a)

这会在分析工具中插入标记，收集该范围内的内核执行数据。

解析与解读指标

在启用分析后，典型输出会展示如下关键指标：

  • 每 SM 活跃 Warp 数（Active Warps per SM）：表示同时运行的 Warp 数（每 Warp 32 线程），越高越接近最大并行度。
  • 每线程寄存器数（Registers per Thread）：寄存器占用越高，可并行块数越少，可能限制占用率。
  • 每块共享内存用量（Shared Memory per Block）：若使用过多共享内存，能同时调度的块数会减少，降低占用率。
  • 占用率（Occupancy）：通常以百分比或比值形式出现（如 0.75），表示实际并行度与理论最大并行度之比。

例如，Nsight Compute 可能报告：

Grid Size:             39063 x 1 x 1  
Block Size:            256 x 1 x 1  
Registers per Thread:  32  
Shared Memory per Block: 2.0 KB  
Active Warps per SM:   24  
Occupancy:             75%

通过这些数据，我们可以回答：

  • 是否每线程寄存器用得太多，限制了块并发？
  • 是否每块共享内存占用过高，减少了活动块？
  • 活跃 Warp 太少，导致吞吐不足？

若占用率因寄存器瓶颈过低，可尝试减少内核中使用的变量或开启编译器标志降低寄存器使用；若共享内存限制占用率，可减少每块共享内存用量或缩小分块尺寸。通过反复测量与调整，即可将 kernel 性能提升至最佳。

总结

本章聚焦于 Kernel 优化的实战技巧以及如何最大化利用 GPU 硬件。我们首先探讨了线程块（block）与网格（grid）尺寸的选择，及其对占用率、并行吞吐量与硬件资源利用率的影响。通过动手实验和 NVIDIA 占用率计算器，我们学会了在工作量、内存消耗和计算能力之间平衡，选择最佳的 block/grid 配置。

接着，我们通过共享内存缓冲技术演示了如何减少全局内存访问，利用高速共享缓存显著提升性能。随后，我们检查了全局内存访问模式，发现未对齐或不规则的访问会严重浪费带宽并拖慢 Kernel。

最后，我们借助 CuPy 的分析工具收集了运行时的 Kernel 统计信息——包括占用率、寄存器和共享内存使用量，以及每个多处理器上的活跃 Warp 数——从而深入了解内核瓶颈。

综上，这些方法帮助我们编写高效、资源感知的 GPU 程序，充分发挥现代 CUDA 硬件的并行计算能力。

avatar
文章
阅读
粉丝