GPU 编程实战——简单并行模式

数组的逐元素操作

逐元素计算与 CuPy

在并行编程中，最常见且最直接的模式就是逐元素操作。在许多实际问题里，我们需要对大型数据集的每个元素执行相同的算术运算——加、减、乘、除，或更复杂的操作。GPU 天生擅长此类任务：数以千万计的线程可以同时工作，每个线程独立处理一个（或多个）元素，无需线程间协调，使得逐元素计算既高效又天然可扩展。这些操作广泛应用于图像处理、科学模拟、机器学习、金融建模等领域。掌握这一模式，既能加速基本算术运算，也能为更复杂的算法打下简洁可行的基础。

CuPy 为我们提供了类似 NumPy 的接口，并在底层自动调用 GPU 加速。大多数情况下，只需替换导入和数组创建方式，就能将代码从 CPU 平滑迁移到 GPU。所有 NumPy 中的广播规则和语法在 CuPy 中无缝复用，让我们瞬间拥有成千上万 CUDA 核心的算力。

运行示例：大规模逐元素算术

下面通过实战演示：

import cupy as cp

# 在 GPU 上创建两个一千多万个元素的数组
size = 10_000_000
a = cp.random.rand(size).astype(cp.float32)
b = cp.random.rand(size).astype(cp.float32)

# 逐元素运算
sum_result  = a + b                  # 相加
diff_result = a - b                  # 相减
prod_result = a * b                  # 相乘
epsilon     = 1e-8
div_result  = a / (b + epsilon)      # 相除（加小常数以防除零）

所有运算均在 GPU 上完成，每个线程独立处理一个元素。相比于在大数组上使用 NumPy，这些计算速度会快得多，且优势随数据规模增长而显著。

广播的威力

NumPy 和 CuPy 的一大强大功能就是广播（broadcasting），它允许对形状不同但维度兼容的数组进行逐元素运算。例如，将一个一维数组加到二维数组的每一行：

# 在 GPU 上创建一个 (size, 100) 的二维数组
array_2d = cp.random.rand(size, 100).astype(cp.float32)
# 将一维数组 a 加到二维数组的每一行
broadcasted_result = array_2d + a[:, None]

CuPy 会在 GPU 上完成所有运算，既保留了 NumPy 的便捷性，又大幅提升了速度。

• 逐元素操作 是并行编程的核心构建块，与 GPU 硬件高度契合。
• CuPy 提供了与 NumPy 几乎相同的语法和广播能力，让 GPU 加速代码编写变得轻而易举。
• 我们可以通过对比 CPU 与 GPU 的输出，快速验证正确性与性能提升。

掌握了这一模式后，我们就为更高级的并行模式做好了准备。

并行归约求和

并行归约概览

在许多场景下，我们需要将大量数据归约为一个标量值。在 CPU 上，通常用循环或内置的 sum() 就能搞定；而在 GPU 上，如果让所有线程同时向同一个全局变量累加，既会产生严重的竞争，也会造成同步瓶颈。并行归约（parallel reduction）通过分阶段“折半”累加来解决这一问题：每个线程先在局部对元素求和，然后再合并这些部分和，直到最后仅剩一个结果。这种做法最大化利用了 GPU 的共享内存（shared memory）高速缓存，显著减少全局内存访问和竞争。

编写归约 Kernel

下面示例使用 PyCUDA，在每个线程块内先对所属区间进行部分求和，最终在主机或通过多次 Kernel 调用完成全局归约。

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

# 数据规模
size = 1_048_576  # 1M 元素

# 在主机上生成随机数组并拷贝到设备
host_array   = np.random.rand(size).astype(np.float32)
device_array = gpuarray.to_gpu(host_array)

# CUDA C 归约 kernel
kernel_code = """
extern "C"
__global__ void reduce_sum(float *g_idata, float *g_odata, int n)
{
    extern __shared__ float sdata[];
    unsigned int tid = threadIdx.x;
    unsigned int i   = blockIdx.x * blockDim.x * 2 + threadIdx.x;
    float sum = 0;

    // 每个线程最多读取两个元素进行累加（展开）
    if (i < n) 
        sum += g_idata[i];
    if (i + blockDim.x < n)
        sum += g_idata[i + blockDim.x];

    sdata[tid] = sum;
    __syncthreads();

    // 在共享内存中迭代折半归约
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // 将每个区块的部分和写回全局内存
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
"""
mod        = SourceModule(kernel_code)
reduce_sum = mod.get_function("reduce_sum")

# 配置线程块和网格
threads_per_block = 512
blocks_per_grid   = (size + threads_per_block * 2 - 1) // (threads_per_block * 2)
shared_mem_size   = threads_per_block * np.dtype(np.float32).itemsize

# 为每个区块的部分和分配输出数组
partial_sums = gpuarray.empty(blocks_per_grid, dtype=np.float32)

# 启动归约 kernel
reduce_sum(
    device_array, partial_sums, np.int32(size),
    block=(threads_per_block, 1, 1),
    grid=(blocks_per_grid, 1),
    shared=shared_mem_size
)

# 如果 blocks_per_grid > 1，可重复上述过程，或直接拷回主机并用 CPU 完成最后的归约
final_sum = partial_sums.get().sum()

print("GPU 归约结果:", final_sum)
print("CPU 参考结果:", host_array.sum())

同步点与优化

• 在每次归约迭代中，我们使用 __syncthreads() 确保所有线程完成对共享内存的读写，防止竞态。
• 通过“展开”（unrolling）让每个线程在第一步读取多达两个元素，减少循环次数并提高对全局内存的利用率。
• 采用共享内存折半归约，大幅降低全局内存访问次数。

若要处理更大规模数据，只需增加区块数量，并可递归地对 partial_sums 再次归约，直到只剩一个区块。与串行归约相比，这种方法在大数组上能带来显著的加速效果。

使用共享内存的分块矩阵乘法

我们之前了解到 CUDA 有多种内存空间，其中共享内存（shared memory）对优化尤为重要。它是每个流式多处理器（SM）上由用户管理的一块小型、超高速的内存区域，比全局内存快得多，并且同一线程块内的所有线程都可以访问。共享内存允许块内线程协作加载数据、重复使用中间结果，从而减少全局内存访问、提升计算密度。由于每个块可用的共享内存通常受限（现代 GPU 通常在 48 KB 或 96 KB 左右），最佳用例是那些被频繁重用的数据，例如在矩阵乘法中的分块（tiles）。

矩阵乘法是科学计算、图形学和机器学习中的基础操作。朴素的 GPU 实现会在全局内存中多次读取同一个元素，导致极高的内存流量并成为性能瓶颈。通过分块——将大矩阵按小子矩阵划分，每个子矩阵加载到共享内存——我们确保每个元素在对应子块中只需从全局内存读取一次，然后由本块内所有线程重用。该策略大幅削减全局内存访问。

分块矩阵乘法示例

下面演示如何对两个 N×N 的方阵进行分块乘法，并利用共享内存加速。

准备数据

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

N = 1024      # 矩阵大小
TILE_SIZE = 32  # 分块大小，最好能整除 N

# 在主机上生成随机矩阵，并拷贝到设备
A_host = np.random.rand(N, N).astype(np.float32)
B_host = np.random.rand(N, N).astype(np.float32)
A_gpu = gpuarray.to_gpu(A_host)
B_gpu = gpuarray.to_gpu(B_host)
C_gpu = gpuarray.zeros((N, N), dtype=np.float32)

CUDA Kernel（共享内存分块）

kernel_code = f"""
#define TILE_SIZE {TILE_SIZE}
__global__ void matmul_tiled(float *A, float *B, float *C, int N)
{{
    __shared__ float sA[TILE_SIZE][TILE_SIZE];
    __shared__ float sB[TILE_SIZE][TILE_SIZE];
    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float sum = 0.0f;

    for (int t = 0; t < N / TILE_SIZE; ++t) {{
        // 协作式加载到共享内存
        sA[threadIdx.y][threadIdx.x] = A[row * N + t * TILE_SIZE + threadIdx.x];
        sB[threadIdx.y][threadIdx.x] = B[(t * TILE_SIZE + threadIdx.y) * N + col];
        __syncthreads();

        // 计算子块内的部分和
        for (int k = 0; k < TILE_SIZE; ++k) {{
            sum += sA[threadIdx.y][k] * sB[k][threadIdx.x];
        }}
        __syncthreads();  // 同步后加载下一对子块
    }}

    C[row * N + col] = sum;
}}
"""
mod = SourceModule(kernel_code)
matmul_tiled = mod.get_function("matmul_tiled")

启动 Kernel

block = (TILE_SIZE, TILE_SIZE, 1)
grid  = (N // TILE_SIZE, N // TILE_SIZE, 1)

matmul_tiled(
    A_gpu, B_gpu, C_gpu, np.int32(N),
    block=block,
    grid=grid
)

验证结果

C_host = C_gpu.get()
C_ref  = np.dot(A_host, B_host)
print("结果一致吗？", np.allclose(C_host, C_ref, atol=1e-3))

每个线程块将对应的 A 和 B 子块从全局内存加载到共享内存，之后本块内线程可多次复用这些数据，而无需再访问全局内存。相比于每个线程在全局内存中多次读取同一行或列的元素，分块策略显著减少了慢速内存访问次数，大多数访问都命中高速的共享内存，从而大幅提升了矩阵乘法的性能。

Python 中的并行映射模式

虽然逐元素算术很强大，但很多实际应用需要更复杂的变换：非线性函数、自定义查表、条件逻辑，甚至多步链式操作。Python 内置的 map() 函数或列表推导式（如 [f(x) for x in xs]）提供了一种简洁表达的方式，对集合中的每个元素应用函数。但这些模式是串行执行的，无法利用现代 GPU 的并行能力。当数据量增大或函数计算开销提升时，串行执行的性能瓶颈就会十分明显。

并行映射模式

并行映射模式（Parallel Map）为此提供了一种通用解决方案：在成千上万的 GPU 线程上，将用户自定义函数并行地应用到每个数据元素上。每个线程接受自己的输入值、计算结果并写入输出数组。它与 Python 的 map() 功能相似，但具备 GPU 级别的规模和速度，能为适合的工作负载带来巨大的性能提升。

并行映射模式的优势在于灵活性——无论是简单的数学变换、数据归一化，还是任何无需元素间交互的操作，都可以用同样的方法实现。它还让 GPU 编程更加“Python 式”，将函数式编程的表达力带入高性能计算领域。

将映射模式映射到 GPU Kernel

下面通过实战示例演示如何用 PyCUDA 将自定义函数并行应用到数组的每个元素上——就像 Python 的 map()，但在 GPU 上充分并行处理。

定义自定义函数

假设我们要对每个元素应用机器学习中常用的 Sigmoid 激活函数：

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

我们不在 Python 中循环，而是在 CUDA Kernel 中实现这个函数。

准备数据

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

N = 10_000_000
input_host = np.random.randn(N).astype(np.float32)
input_gpu  = gpuarray.to_gpu(input_host)
output_gpu = gpuarray.empty_like(input_gpu)

编写并行映射 Kernel

kernel_code = """
__global__ void parallel_map(const float *input, float *output, int n)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n)
    {
        float x = input[idx];
        // 应用 sigmoid 函数
        output[idx] = 1.0f / (1.0f + expf(-x));
    }
}
"""
mod          = SourceModule(kernel_code)
parallel_map = mod.get_function("parallel_map")

启动 Kernel

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

parallel_map(
    input_gpu, output_gpu, np.int32(N),
    block=(threads_per_block, 1, 1),
    grid =(blocks_per_grid,   1   )
)

拷回主机并验证

output_host = output_gpu.get()
sigmoid_np  = 1.0 / (1.0 + np.exp(-input_host))
print("结果吻合吗？", np.allclose(output_host, sigmoid_np, atol=1e-6))

这个模式可以轻松扩展到多输入、复杂表达式或带参数的函数。无论是 ReLU、tanh、阈值处理，还是任意数学函数，都可将逻辑写入 Kernel，然后并行执行千万级元素。使用并行映射，我们只需关注变换本身，而无需为循环或性能担忧。

基础直方图计算

直方图是数据科学、图像处理及其他科学应用中的重要工具，可帮助我们了解数据分布、分析图像像素强度或统计数据集中的频率。对于小规模数组，NumPy 或 CPU 方法通常足够快速。但当数据规模达到数百万或数十亿，或在需要实时性能的场景（如流媒体视频、大规模分析）中，CPU 往往难以胜任。GPU 擅长并行处理大量操作，非常适合加速直方图计算，前提是我们能妥善解决并行更新带来的竞争问题。

使用 CuPy 实现 GPU 直方图

准备数据

假设我们有一组整数值（如像素强度）：

import cupy as cp
import numpy as np

size  = 10_000_000
nbins = 256  # 灰度值 0–255

# 随机生成数据并在 GPU 上分配直方图数组
data      = cp.random.randint(0, nbins, size, dtype=cp.int32)
histogram = cp.zeros(nbins, dtype=cp.int32)

编写直方图 Kernel

每个线程处理数据的一部分，并使用原子操作在共享内存（shared memory）中更新相应的 bin，以减少对全局内存的竞争。最后，再将每个块的部分直方图合并到全局直方图中：

kernel_code = f"""
extern "C" __global__
void histogram(const int *data, int *hist, int N, int nbins)
{{
    extern __shared__ int s_hist[];
    int tid = threadIdx.x;
    int gid = blockIdx.x * blockDim.x + threadIdx.x;

    // 初始化共享直方图
    for (int i = tid; i < nbins; i += blockDim.x)
        s_hist[i] = 0;
    __syncthreads();

    // 在共享内存中累加
    if (gid < N) {{
        int bin = data[gid];
        if (bin >= 0 && bin < nbins)
            atomicAdd(&s_hist[bin], 1);
    }}
    __syncthreads();

    // 将共享直方图合并到全局直方图
    for (int i = tid; i < nbins; i += blockDim.x)
        atomicAdd(&hist[i], s_hist[i]);
}}
"""
module      = cp.RawModule(code=kernel_code)
hist_kernel = module.get_function("histogram")

配置并启动 Kernel

threads_per_block = 256
blocks_per_grid   = (size + threads_per_block - 1) // threads_per_block
shared_mem_size   = nbins * np.dtype(np.int32).itemsize

hist_kernel(
    (blocks_per_grid,), (threads_per_block,),
    (data, histogram, size, nbins),
    shared_mem=shared_mem_size
)

拷回主机并验证

# 参考 CPU 版本直方图
hist_cpu = np.histogram(cp.asnumpy(data), bins=nbins, range=(0, nbins))[0]
# GPU 直方图
hist_gpu = histogram.get()

print("GPU 与 CPU 直方图一致吗？", np.array_equal(hist_cpu, hist_gpu))

上述代码中，每个线程块先在共享内存中构建自己的直方图片段，再通过原子加到全局直方图。这种方法在处理海量数据和多达数百个 bins 时，既能保证正确性，又具备良好可扩展性。

直方图计算是典型的数据密集型任务，通过共享内存与原子操作的精心设计，GPU 加速可显著提升性能。借助 CuPy 的 RawKernel，我们能够高效实现这一经典并行模式，从而以高速和大规模处理数据集、图像或流式数据。

总结

总而言之，我们已经掌握了多个常用的 GPU 并行模式及其实践实现：

1. 逐元素操作：利用 CuPy 的广播机制和类 NumPy 语法，在 GPU 上对大数组执行高吞吐量算术运算。
2. 并行归约：使用 PyCUDA 编写高效的求和归约 kernel，通过共享内存和精心同步，实现大规模数据的快速汇总，轻松扩展到各种归约任务。
3. 分块矩阵乘法：利用 CUDA 共享内存对矩阵进行分块（tiling），大幅减少全局内存访问，从而显著提升矩阵运算的吞吐量。
4. 并行映射：将 Python 的函数式 map() 模式推广到 GPU，编写自定义 kernel 并行应用任意元素级函数，实现复杂变换。
5. 直方图计算：结合共享内存与原子操作，在每个线程块内安全并发更新直方图子段，再合并到全局直方图，轻松处理数百万线程的并发更新。

通过这些模式，我们不仅提高了对高性能并行原语的理解，还学会了如何在实际工作中，应对大规模数据处理和复杂算法的挑战。