CUDA 深度解析

Appearance

并行算法�?GPU 优化策略 ​

GPU 优化不只是调参，更需要从算法层面重新思考问题。本文梳�?GPU 并行算法的核心模式和系统性优化策略�?

GPU 优化的层�? ​

优化层次（从高到低，收益从大到小）：

1. 算法选择        �?选择适合 GPU 的算法（10-100x�?
2. 内存访问模式    �?合并访问、共享内存（2-10x�?
3. 计算强度        �?减少内存访问，增加计算复用（2-5x�?
4. Occupancy      �?提高 SM 利用率（1.2-2x�?
5. 指令级优�?     �?使用快速数学函数（1.1-1.5x�?

核心并行模式 ​

1. Map（映射） ​

最简单的并行模式：每个线程独立处理一个元素�?

cpp
// 向量加法：完美的 Map 操作
__global__ void vectorAdd(float* a, float* b, float* c, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) c[i] = a[i] + b[i];
}

2. Reduce（归约） ​

�?N 个元素归约为 1 个结果，是最重要的并行模式之一�?

cpp
// 朴素归约（有 Bank Conflict �?Warp Divergence�?
__global__ void reduceNaive(float* data, float* result, int N) {
    __shared__ float smem[256];
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;
    
    smem[tid] = (idx < N) ? data[idx] : 0.0f;
    __syncthreads();
    
    // 有分歧：stride 为奇数时 Warp 内分�?
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        if (tid % (2 * stride) == 0) {
            smem[tid] += smem[tid + stride];
        }
        __syncthreads();
    }
    if (tid == 0) result[blockIdx.x] = smem[0];
}

// 优化归约（无分歧，无 Bank Conflict�?
__global__ void reduceOptimized(float* data, float* result, int N) {
    __shared__ float smem[256];
    int tid = threadIdx.x;
    int idx = blockIdx.x * (blockDim.x * 2) + tid;
    
    // 每个线程加载两个元素（减�?Block 数量�?
    smem[tid] = (idx < N ? data[idx] : 0.0f) +
                (idx + blockDim.x < N ? data[idx + blockDim.x] : 0.0f);
    __syncthreads();
    
    // 从大步长开始，无分�?
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) smem[tid] += smem[tid + stride];
        __syncthreads();
    }
    
    // Warp 内归约（无需 __syncthreads�?
    if (tid < 32) {
        volatile float* vsmem = smem;
        vsmem[tid] += vsmem[tid + 32];
        vsmem[tid] += vsmem[tid + 16];
        vsmem[tid] += vsmem[tid + 8];
        vsmem[tid] += vsmem[tid + 4];
        vsmem[tid] += vsmem[tid + 2];
        vsmem[tid] += vsmem[tid + 1];
    }
    
    if (tid == 0) result[blockIdx.x] = smem[0];
}

Warp Shuffle 归约（现代最优方案） ​

cpp
__device__ float warpReduceSum(float val) {
    // Warp �?Shuffle 归约，无需共享内存
    for (int offset = 16; offset > 0; offset >>= 1)
        val += __shfl_down_sync(0xffffffff, val, offset);
    return val;
}

__global__ void reduceWarpShuffle(float* data, float* result, int N) {
    float sum = 0.0f;
    
    // 每个线程处理多个元素（Grid-stride loop�?
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N;
         i += blockDim.x * gridDim.x) {
        sum += data[i];
    }
    
    // Warp 内归�?
    sum = warpReduceSum(sum);
    
    // 每个 Warp 的结果写入共享内�?
    __shared__ float warpSums[32];
    int lane = threadIdx.x % 32;
    int warpId = threadIdx.x / 32;
    
    if (lane == 0) warpSums[warpId] = sum;
    __syncthreads();
    
    // 第一�?Warp 归约所�?Warp 的结�?
    if (warpId == 0) {
        sum = (threadIdx.x < blockDim.x / 32) ? warpSums[lane] : 0.0f;
        sum = warpReduceSum(sum);
        if (lane == 0) atomicAdd(result, sum);
    }
}

3. Scan（前缀扫描�? ​

cpp
// Blelloch 并行前缀扫描（Work-Efficient�?
__global__ void exclusiveScan(int* data, int* output, int N) {
    __shared__ int smem[512];
    int tid = threadIdx.x;
    
    smem[tid] = (tid < N) ? data[tid] : 0;
    __syncthreads();
    
    // Up-sweep（归约树�?
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        int idx = (tid + 1) * stride * 2 - 1;
        if (idx < blockDim.x)
            smem[idx] += smem[idx - stride];
        __syncthreads();
    }
    
    // 清零根节�?
    if (tid == 0) smem[blockDim.x - 1] = 0;
    __syncthreads();
    
    // Down-sweep（分发树�?
    for (int stride = blockDim.x / 2; stride >= 1; stride /= 2) {
        int idx = (tid + 1) * stride * 2 - 1;
        if (idx < blockDim.x) {
            int tmp = smem[idx - stride];
            smem[idx - stride] = smem[idx];
            smem[idx] += tmp;
        }
        __syncthreads();
    }
    
    if (tid < N) output[tid] = smem[tid];
}

4. Histogram（直方图�? ​

cpp
// 使用共享内存原子操作优化直方�?
__global__ void histogram(int* data, int* hist, int N, int numBins) {
    __shared__ int localHist[256];
    
    // 初始化局部直方图
    if (threadIdx.x < numBins)
        localHist[threadIdx.x] = 0;
    __syncthreads();
    
    // 局部统计（共享内存原子，比全局内存�?10x�?
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        atomicAdd(&localHist[data[idx] % numBins], 1);
    }
    __syncthreads();
    
    // 合并到全局直方�?
    if (threadIdx.x < numBins)
        atomicAdd(&hist[threadIdx.x], localHist[threadIdx.x]);
}

关键优化技�? ​

1. Grid-Stride Loop ​

让每个线程处理多个元素，提高灵活性和效率�?

cpp
__global__ void processLargeArray(float* data, int N) {
    // Grid-stride loop：线程数不必等于数据�?
    for (int i = blockIdx.x * blockDim.x + threadIdx.x;
         i < N;
         i += blockDim.x * gridDim.x) {
        data[i] = process(data[i]);
    }
}

// 启动时可以用固定�?Grid 大小
int numSMs;
cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
processLargeArray<<<numSMs * 32, 256>>>(d_data, N);

2. 向量化内存访�? ​

使用 float4/int4 等向量类型，一次加�?128 位：

cpp
__global__ void vectorizedLoad(float* data, float* output, int N) {
    int idx = (blockIdx.x * blockDim.x + threadIdx.x) * 4;
    
    if (idx + 3 < N) {
        // 一次加�?4 �?float�?28 位，一�?Cache Line 事务�?
        float4 val = reinterpret_cast<float4*>(data)[idx / 4];
        val.x *= 2.0f;
        val.y *= 2.0f;
        val.z *= 2.0f;
        val.w *= 2.0f;
        reinterpret_cast<float4*>(output)[idx / 4] = val;
    }
}

3. 循环展开 ​

cpp
// 手动展开（减少循环控制开销�?
#pragma unroll 4
for (int k = 0; k < K; k++) {
    sum += A[row * K + k] * B[k * N + col];
}

// 完全展开（编译时已知循环次数�?
#pragma unroll
for (int k = 0; k < TILE_SIZE; k++) {
    sum += tileA[ty][k] * tileB[k][tx];
}

4. 快速数学函�? ​

cpp
// 标准精度（慢�?
float y = sinf(x);
float z = sqrtf(x);
float w = x / y;

// 快速近似（精度略低，但�?2-4x�?
float y = __sinf(x);      // �?4 ULP 误差
float z = __fsqrt_rn(x);  // 舍入到最�?
float w = __fdividef(x, y); // 快速除�?

// 编译器全局开启快速数�?
// nvcc --use_fast_math my_kernel.cu

5. 寄存器复用与指令级并�? ​

cpp
// 差：依赖链长，无法并�?
float a = load(0);
float b = a * 2.0f;
float c = b + 1.0f;
float d = c * c;

// 好：独立计算，编译器可以并行发射
float a = load(0);
float b = load(1);
float c = load(2);
float d = load(3);
// 四条加载指令可以并行执行
float r0 = a * 2.0f;
float r1 = b * 3.0f;
float r2 = c * 4.0f;
float r3 = d * 5.0f;

Tiling（分块）优化 ​

分块是提高计算强度的核心技术，通过共享内存减少全局内存访问�?

不分块的矩阵乘法�?
  每个输出元素需要读�?2K 个全局内存元素
  算术强度 = 2K FLOP / (2K × 4 Byte) = 0.25 FLOP/Byte

分块（Tile=16）的矩阵乘法�?
  每个 Tile 的数据被 16 个线程复�?
  算术强度 = 2K FLOP / (2K/16 × 4 Byte) = 4 FLOP/Byte
  提升 16x�?

性能优化检查清�? ​

�?内存访问
  �?Warp 内线程访问连续内存（合并访问�?
  �?热点数据放入共享内存
  �?避免共享内存 Bank Conflict
  �?使用向量化加载（float4�?

�?计算效率
  �?避免 Warp Divergence
  �?使用 Warp Shuffle 代替共享内存归约
  �?循环展开�?pragma unroll�?
  �?使用快速数学函�?

�?并发与流水线
  �?使用�?Stream 重叠计算与传�?
  �?使用固定内存（Pinned Memory�?
  �?考虑 CUDA Graph 减少 CPU 开销

�?资源利用
  �?Occupancy �?50%（用 Nsight 测量�?
  �?Block 大小�?32 的倍数
  �?控制寄存器使用量（避免溢出）

下一篇：Nsight 性能分析工具 →

基于 NVIDIA CUDA 官方文档整理

CUDA Deep Dive © 2024