深入浅出CUDA：GPU并行计算的入门指南

在AI和大数据时代，GPU已成为加速计算的核心力量。CUDA（Compute Unified Device Architecture）作为NVIDIA推出的GPU编程框架，让开发者能利用GPU的并行能力处理复杂任务。

CUDA的核心基础概念

CUDA的并行计算模型建立在线程的层次组织上。这种结构允许高效分配任务到GPU的多个计算单元。

线程层次：Thread、Block和Grid

Thread（线程）：CUDA中最小的执行单位。每个线程执行相同的内核代码，但处理不同数据。通过内置变量threadIdx.x、threadIdx.y、threadIdx.z标识其在块内的位置。线程是并行计算的基本元素，适合处理独立的数据元素。
Block（线程块）：一组线程的集合，通常包含128到1024个线程（最大1024）。块内线程可以共享内存和同步执行，由blockIdx.x、blockIdx.y、blockIdx.z和blockDim.x、blockDim.y、blockDim.z定义位置和大小。块是GPU调度的基本单位，被分配到Streaming Multiprocessor (SM)上执行。块的维度可以是1D、2D或3D，适应不同任务如向量、矩阵或体积数据。
Grid（网格）：多个块的集合，可达数万个块，形成一个大网格。由gridDim.x、gridDim.y、gridDim.z定义大小。网格负责覆盖整个计算任务，例如在大型矩阵运算中，网格维度对应矩阵的行和列块数。

这些层次的关系是嵌套的：线程组成块，块组成网格。GPU硬件将块分配到多个SM上并发执行，线程间协作通过同步指令如__syncthreads()实现。实际调度中，块的大小影响资源利用率，而网格的大小确保任务覆盖。

以下是PyCUDA中一个简单示例，展示线程层次在向量加法中的应用：

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy as np

# 定义内核代码
mod = SourceModule("""
__global__ void vector_add(float *a, float *b, float *c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}
""")

# 准备数据
n = 1024
a = np.random.randn(n).astype(np.float32)
b = np.random.randn(n).astype(np.float32)
c = np.zeros_like(a)

# 分配设备内存
a_gpu = cuda.mem_alloc(a.nbytes)
b_gpu = cuda.mem_alloc(b.nbytes)
c_gpu = cuda.mem_alloc(c.nbytes)

# 拷贝数据到设备
cuda.memcpy_htod(a_gpu, a)
cuda.memcpy_htod(b_gpu, b)

# 定义块和网格：块大小256线程（1D），网格大小(n+255)//256
block_size = (256, 1, 1)
grid_size = ((n + 255) // 256, 1, 1)

# 调用内核
func = mod.get_function("vector_add")
func(a_gpu, b_gpu, c_gpu, np.int32(n), block=block_size, grid=grid_size)

# 拷贝结果回主机
cuda.memcpy_dtoh(c, c_gpu)

# 验证
print(np.allclose(c, a + b))  # True

Warp（线程束）：SIMT架构的核心

**Warp是32个连续线程的组（NVIDIA主流架构的标准）。**GPU以Warp为单位调度指令，这源于SIMT（Single Instruction Multiple Threads）模型：Warp内所有线程执行同一指令，但操作不同数据。SIMT允许高效的单指令多数据处理。

SIMT与CPU的区别在于架构设计

CPU优化低延迟，核心少但复杂，适合分支密集的串行任务
GPU优化高吞吐，核心多但简单，适合数据并行任务，通过海量线程和Warp切换掩盖内存延迟

如果Warp内线程执行路径不同（如条件分支），发生分支发散（Divergence）：硬件使用掩码串行执行分支路径，导致部分线程闲置，效率下降。开发者需设计统一执行路径的代码。

在PyCUDA中，Divergence可能在条件语句中出现。以下示例展示潜在Divergence：

mod = SourceModule("""
__global__ void conditional_add(float *a, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        if (a[idx] > 0) {  // 可能导致Divergence
            a[idx] += 1.0;
        } else {
            a[idx] -= 1.0;
        }
    }
}
""")

为避免Divergence，可用掩码替换：a[idx] += (a[idx] > 0 ? 1.0 : -1.0);，确保Warp统一执行。

CUDA在并行计算中的生态位

CUDA是GPU通用计算的事实标准，其价值在于整合硬件、软件和生态。

核心价值与优势

CUDA标准化GPU编程：定义线程层次、内存模型和异构计算范式（CPU控制逻辑、GPU执行并行）。它允许直接利用GPU进行通用计算，而非早期通过图形API间接实现。相比OpenCL，CUDA的优势在于NVIDIA硬件优化，如访问Tensor Core（矩阵加速单元）和NVLink（高带宽互联），性能更高。OpenCL提供跨平台兼容，但牺牲专属优化，适用于AMD/Intel GPU场景。

在适用性上，CUDA专注高性能：AI训练加速矩阵乘（GEMM）和卷积；科学计算处理线性代数、FFT和模拟；图形渲染利用并行像素处理；深度学习框架通过CUDA实现加速，如PyTorch的torch.cuda模块。

生态影响力

CUDA生态形成正向循环：硬件（如Ampere、Hopper、Blackwell）持续支持CUDA新特性（如TMA异步拷贝）；库如cuDNN（深度神经网络）、cuBLAS（线性代数）、cuFFT（傅里叶变换）和NCCL（多GPU通信）封装优化算法；开发者社区贡献开源项目，如CUTLASS（自定义GEMM模板）和Triton（内核DSL）。在大模型训练中，NCCL处理AllReduce同步，FlashAttention用CUDA优化注意力机制。CUDA影响力扩展到量子计算（cuQuantum）和光刻模拟（cuLitho）。PyCUDA/Torch等绑定让Python开发者无缝集成。

CUDA计算的基本流程

CUDA计算采用异构模式：主机（CPU）协调，设备（GPU）执行。流程包括初始化、数据传输、内核调用、同步和释放。流程强调最小化交互：批量操作减少拷贝。

初始化与内存分配：主机用cudaGetDeviceCount()和cudaGetDeviceProperties()查询设备，cudaMalloc()分配显存。这步检查SM版本和能力。
数据传输（H2D）：cudaMemcpy()拷贝主机到设备。异步版cudaMemcpyAsync()与Stream结合，重叠传输。
Kernel调用：主机用<<<gridDim, blockDim>>>启动__global__内核。内核在GPU执行，如向量加法。
同步与结果回传（D2H）：cudaDeviceSynchronize()等待完成，cudaMemcpy()回传。
资源释放：cudaFree()释放显存，避免泄漏。

提升CUDA计算效率的实用技巧

优化聚焦并行度、内存和指令效率。

线程块与网格的合理划分

Block大小选Warp倍数（128-1024），如256；网格大小(n + block_size - 1) / block_size覆盖任务。原理：平衡调度和资源。高Occupancy（活跃Warp比率）需考虑寄存器和Shared用量。用Nsight Compute检查。

Warp调度优化

避免Divergence：用三元运算符统一路径。__syncwarp()同步Warp。确保Coalesced访问：线程连续地址。

核函数优化思路

简化代码：少分支，float优先。融合内核减少启动（~微秒级）。用Tensor Core加速GEMM。

避免冗余数据交互

批量拷贝，用cudaHostAlloc()Pinned内存加速。Stream重叠：多Stream并发Kernel/拷贝。

Profiler迭代：Nsight Systems看时间线，Nsight Compute看瓶颈。

以下PyCUDA矩阵乘示例展示优化：

mod = SourceModule("""
__global__ void matrix_mul(float *a, float *b, float *c, int width) {
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int bx = blockIdx.x;
    int by = blockIdx.y;

    float sum = 0.0;
    for (int k = 0; k < width; ++k) {
        sum += a[by * width + k] * b[k * width + bx];
    }
    c[by * width + bx] = sum;
}
""")

# 优化版用Shared tiling（后续内存节扩展）

CUDA的分层内存模型与优化

内存类型详解

内存层次从快到慢：寄存器、Shared、L1/L2、Global

寄存器：线程私有，~1周期，~256KB/SM。局部变量。用多Spill到Local（慢）。
共享内存：块内共享，~几周期，~100-192KB/SM。声明__shared__，用于Tiling。避免Bank Conflict（32 Banks）。
全局内存：GB级，数百周期，全局可见。优化Coalesced。
常量内存：64KB，只读广播。
纹理内存：2D/3D优化，带缓存。

内存优化策略

频繁数据放寄存器/Shared，大数据Global。统一内存简化但页面错误慢。Shared Tiling减少Global访存：分块加载。

PyCUDA Tiling示例：

mod = SourceModule("""
#define TILE_WIDTH 16
__global__ void matrix_mul_tiled(float *a, float *b, float *c, int width) {
    __shared__ float tileA[TILE_WIDTH][TILE_WIDTH];
    __shared__ float tileB[TILE_WIDTH][TILE_WIDTH];

    int tx = threadIdx.x, ty = threadIdx.y;
    int bx = blockIdx.x, by = blockIdx.y;

    float sum = 0.0;
    for (int t = 0; t < (width + TILE_WIDTH - 1) / TILE_WIDTH; ++t) {
        if (by * TILE_WIDTH + ty < width && t * TILE_WIDTH + tx < width) {
            tileA[ty][tx] = a[(by * TILE_WIDTH + ty) * width + (t * TILE_WIDTH + tx)];
        } else {
            tileA[ty][tx] = 0.0;
        }
        if (bx * TILE_WIDTH + tx < width && t * TILE_WIDTH + ty < width) {
            tileB[ty][tx] = b[(t * TILE_WIDTH + ty) * width + (bx * TILE_WIDTH + tx)];
        } else {
            tileB[ty][tx] = 0.0;
        }
        __syncthreads();

        for (int k = 0; k < TILE_WIDTH; ++k) {
            sum += tileA[ty][k] * tileB[k][tx];
        }
        __syncthreads();
    }
    if (by * TILE_WIDTH + ty < width && bx * TILE_WIDTH + tx < width) {
        c[(by * TILE_WIDTH + ty) * width + (bx * TILE_WIDTH + tx)] = sum;
    }
}
""")

Tiling将矩阵分块到Shared，减少Global访问，提高带宽利用。

Ref

https://developer.nvidia.cn/blog/even-easier-introduction-cuda-2/