GPU计算：这次怎么落地的

GPU计算 CUDA 性能优化 并行计算
作者: 醉月思 发布日期: 2021-03-06 分类: 技术实践

GPU计算相关的坑，多半出在边界条件上。

最开始接触GPU计算是因为一个图像处理任务：对4K视频做实时边缘检测。

为什么要折腾GPU

最开始接触GPU计算是因为一个图像处理任务：对4K视频做实时边缘检测。单线程在CPU上处理一帧要200ms，完全达不到实时要求。第一次尝试用Python + NumPy的多线程优化，勉强降到50ms，但还是不够。这时候才认真考虑GPU。

选择CUDA而不是OpenCL主要是因为当时手头只有NVIDIA显卡，而且CUDA生态相对成熟，文档、社区、工具链都齐全。后来证明这个决定是对的，踩坑时容易找到前人的经验。

环境搭建：第一次上机的坑

环境搭建本身不复杂，但有几个地方容易踩坑。当时用的是Ubuntu 20.04 + RTX 3060，CUDA版本选了11.8（那时候12.x刚出不稳定）。

# 先检查NVIDIA驱动是否正常
nvidia-smi

# 如果驱动版本过低，可能需要手动更新
sudo apt update
sudo apt install nvidia-driver-525

# 安装CUDA Toolkit（不安装驱动，避免冲突）
wget https://developer.download.nvidia.com/compute/cuda/11.8.0/local_installers/cuda_11.8.0_520.61.05_linux.run
sudo sh cuda_11.8.0_520.61.05_linux.run --toolkit --silent --override

# 配置环境变量
echo 'export PATH=/usr/local/cuda/bin:$PATH' >> ~/.bashrc
echo 'export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH' >> ~/.bashrc
source ~/.bashrc

# 验证安装
nvcc --version

第一个坑是驱动和CUDA版本的匹配。CUDA 11.8要求驱动版本至少520，但系统自带的驱动只有470。强制用老版本驱动跑新CUDA会报莫名其妙的问题，比如编译通过但运行时core dump。

第二个坑是环境变量配了但没生效。编译时找不到cudart，链接时报错。解决方案是重启终端或者在当前shell里手动source一下。如果还是不行，检查LD_LIBRARY_PATH是否正确，有时候软链接没建立好。

Hello CUDA：从矩阵乘法开始

第一个CUDA程序通常是"Hello World"或者向量加法，但我选择矩阵乘法是因为它更能看出GPU的威力，也更暴露问题。

// matrix_mul.cu
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>

__global__ void matrixMulNaive(float *A, float *B, float *C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

int main() {
    int N = 1024;
    size_t size = N * N * sizeof(float);

    float *h_A, *h_B, *h_C;
    float *d_A, *d_B, *d_C;

    // 分配主机内存
    h_A = (float*)malloc(size);
    h_B = (float*)malloc(size);
    h_C = (float*)malloc(size);

    // 初始化数据
    for (int i = 0; i < N * N; i++) {
        h_A[i] = 1.0f;
        h_B[i] = 1.0f;
    }

    // 分配设备内存
    cudaMalloc(&d_A, size);
    cudaMalloc(&d_B, size);
    cudaMalloc(&d_C, size);

    // 数据传输到GPU
    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    // 配置执行参数
    dim3 blockSize(16, 16);
    dim3 gridSize((N + blockSize.x - 1) / blockSize.x,
                   (N + blockSize.y - 1) / blockSize.y);

    // 记录时间
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);

    // 执行内核
    matrixMulNaive<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);

    // 等待完成并记录时间
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);

    // 传输结果回主机
    cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

    printf("GPU time: %.2f ms\n", milliseconds);

    // 验证结果（前几个元素应该是N）
    printf("C[0][0] = %f (expected %f)\n", h_C[0], (float)N);
    printf("C[1][1] = %f (expected %f)\n", h_C[N + 1], (float)N);

    // 释放资源
    free(h_A); free(h_B); free(h_C);
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    return 0;
}

编译和运行：

nvcc -O2 matrix_mul.cu -o matrix_mul
./matrix_mul

1024x1024的矩阵乘法，这个naive版本在我的RTX 3060上大概需要30ms左右。相比之下，用OpenBLAS优化后的CPU版本大概需要150ms。看起来还不错，但这里面的问题不少。

第一次性能分析：问题出在哪

用nvprof分析一下性能瓶颈：

nvprof ./matrix_mul

结果发现，大部分时间花在两个地方：全局内存访问和低算术强度。具体来说，矩阵乘法中每个线程要从全局内存读2N次（A和B各N个元素），写1次，但只做N次乘加运算。当N=1024时，算术强度大约是N/3 ≈ 341 FLOP/byte，理论上应该达到很高的计算强度，但实际上内存带宽利用率很低。

这时候才理解：GPU的优势在于大量线程并行和shared memory的利用，而naive版本完全没有利用shared memory，每个线程都在频繁访问全局内存。

Shared Memory优化：第一次加速

把A和B矩阵分块加载到shared memory，减少全局内存访问：

__global__ void matrixMulShared(float *A, float *B, float *C, int N) {
    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;

    // 每个block处理一个tile，tile大小和block一致
    int aBegin = N * BLOCK_SIZE * by;
    int aEnd   = aBegin + N - 1;
    int aStep  = BLOCK_SIZE;

    int bBegin = BLOCK_SIZE * bx;
    int bStep  = BLOCK_SIZE * N;

    float cSub = 0.0f;

    // 共享内存，存放tile
    __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

    for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep) {
        // 协作加载tile到shared memory
        As[ty][tx] = A[a + N * ty + tx];
        Bs[ty][tx] = B[b + N * ty + tx];

        __syncthreads();  // 等待所有线程加载完成

        // 从shared memory计算
        for (int k = 0; k < BLOCK_SIZE; ++k) {
            cSub += As[ty][k] * Bs[k][tx];
        }

        __syncthreads();  // 等待所有线程计算完成
    }

    // 写回结果
    int cIndex = N * BLOCK_SIZE * by + BLOCK_SIZE * bx;
    C[cIndex + N * ty + tx] = cSub;
}

这里BLOCK_SIZE需要和grid/block配置对应：

#define BLOCK_SIZE 16

// 配置执行参数
dim3 blockSize(BLOCK_SIZE, BLOCK_SIZE);
dim3 gridSize((N + BLOCK_SIZE - 1) / BLOCK_SIZE,
               (N + BLOCK_SIZE - 1) / BLOCK_SIZE);

matrixMulShared<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);

优化后的性能提升明显：1024x1024的矩阵乘法从30ms降到了8ms左右。但还没到极限，因为还有进一步优化的空间。

踩过的坑：Bank Conflict

优化到一半发现性能提升不如预期，用nvprof一查，发现bank conflict严重。shared memory被分成32个banks，当warp中的32个线程同时访问同一个bank的不同地址时，就会冲突。

在矩阵乘法的代码里，As[ty][k]的访问模式会导致bank conflict，因为当ty不变、k变化时，线程会访问不同的bank，但如果多个线程访问同一个bank的不同位置，就会冲突。

解决方法是调整shared memory的布局，用padding避免冲突：

__shared__ float As[BLOCK_SIZE][BLOCK_SIZE + 1];  // +1 padding
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE + 1];

// 访问时保持不变
As[ty][tx] = A[a + N * ty + tx];
Bs[ty][tx] = B[b + N * ty + tx];

加上padding后，As[ty][k]的访问会让不同线程落在不同的bank上，减少冲突。这个改动又带来了约20%的性能提升。

向量化和Warp优化

另一个容易忽略的点是warp级别的优化。CUDA以warp（32个线程）为单位执行，让warp内的线程执行相同的指令。如果warp内的线程出现分支，就会性能下降。

在矩阵乘法中，边界检查if (row < N && col < N)会导致warp divergence。解决方案是让grid大小是block大小的整数倍，然后在kernel内部过滤无效线程：

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

// 不做边界检查，但只计算有效元素
float sum = 0.0f;
for (int k = 0; k < N; k++) {
    if (row < N && col < N && k < N) {  // 但这里还是有分支
        sum += A[row * N + k] * B[k * N + col];
    }
}

更好的办法是让grid大小完全对齐，不需要边界检查：

dim3 gridSize((N + blockSize.x - 1) / blockSize.x * blockSize.x,
               (N + blockSize.y - 1) / blockSize.y * blockSize.y);

但这会浪费一些计算资源。实际项目中需要权衡：如果N很大，浪费的比例很小；如果N很小，可能需要重新考虑是否真的需要GPU加速。

内存合并访问：另一个容易被忽略的点

GPU的全局内存访问最好是合并的（coalesced），即warp内的线程访问连续的内存地址。在naive版本中，A[row * N + k]的访问模式是不合并的，因为row相同、k变化，线程访问的是跨stride的地址。

Shared memory版本已经部分解决了这个问题，因为tile是从连续内存加载的。但如果优化到更深层次，比如使用更复杂的tiling策略，就要注意内存访问模式。

一个简单的验证方法是使用nvprof的gld_transactions和gst_transactions指标，检查内存事务数量。理论上，合并访问的事务数量应该远小于非合并访问。

流水线和异步执行

真正的性能优化到后面，会发现计算本身不是瓶颈，而是数据传输和kernel启动的开销。CUDA支持流（stream）和异步执行，可以重叠计算和数据传输。

// 创建多个流
cudaStream_t streams[4];
for (int i = 0; i < 4; i++) {
    cudaStreamCreate(&streams[i]);
}

// 分块处理
int chunkSize = N / 4;
for (int i = 0; i < 4; i++) {
    int offset = i * chunkSize;

    // 在不同流中执行数据传输和计算
    cudaMemcpyAsync(d_A + offset, h_A + offset, chunkSize * N * sizeof(float),
                   cudaMemcpyHostToDevice, streams[i]);
    cudaMemcpyAsync(d_B + offset, h_B + offset, chunkSize * N * sizeof(float),
                   cudaMemcpyHostToDevice, streams[i]);

    dim3 blockSize(16, 16);
    dim3 gridSize((chunkSize + blockSize.x - 1) / blockSize.x,
                   (N + blockSize.y - 1) / blockSize.y);

    matrixMulShared<<<gridSize, blockSize, 0, streams[i]>>>(
        d_A + offset, d_B, d_C + offset, N);

    cudaMemcpyAsync(h_C + offset, d_C + offset, chunkSize * N * sizeof(float),
                   cudaMemcpyDeviceToHost, streams[i]);
}

// 同步所有流
for (int i = 0; i < 4; i++) {
    cudaStreamSynchronize(streams[i]);
    cudaStreamDestroy(streams[i]);
}

这个改动对大规模矩阵乘法效果明显，因为数据传输和计算可以重叠。但对小规模数据可能适得其反，因为流管理本身有开销。

实践中的几个判断

折腾这么多之后，总结几个在项目中做GPU加速时的判断：

什么时候值得用GPU

  1. 计算密集度高，不是I/O bound
  2. 数据规模够大，能喂满GPU
  3. 问题本身有足够并行度
  4. 考虑到GPU加速的迁移成本和调试难度

如果问题不符合这些条件，用CPU可能更省事。有时候优化CPU代码（用SIMD、OpenBLAS、MKL等）能获得更好的性价比。

选择CUDA还是其他方案

CUDA生态成熟，但绑定NVIDIA。如果需要跨平台，可以考虑：

  • OpenCL：支持更多硬件，但开发体验不如CUDA
  • HIP：AMD的方案，类似CUDA，可以转换
  • 高层框架：cuBLAS、cuDNN等，避免自己写kernel
  • Python包装：Numba、PyCUDA、CuPy等

生产环境通常先用高层框架，发现性能瓶颈再考虑手写kernel。

优化的边界在哪里

早期优化投入产出比高，但后期优化可能花几天时间才提升10%。需要权衡：

  • 瓶颈是否真的在GPU kernel
  • 是否有算法层面可以优化的地方
  • 增加复杂度是否值得

有时候换一个算法（比如用快速傅里叶变换代替直接卷积）比优化现有实现更有效。

最后的成果

经过这些优化，1024x1024的矩阵乘法从最初的30ms降到了3ms左右，比naive版本快了10倍。更重要的是，这些经验和思考方式可以迁移到其他GPU加速任务上。

代码最终版本大概长这样：

#define BLOCK_SIZE 16
#define TILE_SIZE 16

__global__ void matrixMulOptimized(float *A, float *B, float *C, int N) {
    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;

    __shared__ float As[TILE_SIZE][TILE_SIZE + 1];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE + 1];

    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float cValue = 0.0f;

    for (int t = 0; t < (N + TILE_SIZE - 1) / TILE_SIZE; ++t) {
        if (row < N && t * TILE_SIZE + tx < N)
            As[ty][tx] = A[row * N + t * TILE_SIZE + tx];
        else
            As[ty][tx] = 0.0f;

        if (t * TILE_SIZE + ty < N && col < N)
            Bs[ty][tx] = B[(t * TILE_SIZE + ty) * N + col];
        else
            Bs[ty][tx] = 0.0f;

        __syncthreads();

        for (int k = 0; k < TILE_SIZE; ++k) {
            cValue += As[ty][k] * Bs[k][tx];
        }

        __syncthreads();
    }

    if (row < N && col < N)
        C[row * N + col] = cValue;
}

这不是最优版本（cuBLAS的矩阵乘法更快），但已经展示了从naive到优化的大部分思路。真正生产环境会用cuBLAS这样的优化库，除非有特殊需求才手写kernel。

小结

从第一次看到GPU加速比数字的兴奋，到后来发现各种问题的沮丧，再到逐步理解GPU工作原理后的从容，这是一段典型的技术成长之路。

GPU加速不是魔法，它只是一个工具，用得好不好关键在于理解它的约束和优势。代码写得再巧妙，如果不符合GPU的执行模型，也跑不出性能；反之，理解了原理，很多"优化"其实是让代码更听话。

最后想说的是，不要为了用GPU而用GPU。CPU在某些场景下依然是更好的选择，有时候换一个算法或者数据结构，比硬推到GPU上更有效。技术选型的核心永远是"解决实际问题"，而不是"用上最酷的技术"。

折腾到现在，遇到GPU相关的任务，心里已经有了基本的判断：什么时候值得用，大概能做到什么程度，可能会遇到什么问题。这种判断比某个具体的kernel代码更有价值，因为它能指导后续的技术决策。

可用性说明：本文发布于 2021 年 3 月，距今已超过五年。文中涉及的软件版本、接口、下载地址、命令参数和操作界面可能已经发生变化，部分方案在当前环境下可能失效。请结合官方最新文档核对后再操作，生产环境使用前务必先行验证。

版权声明: 本文首发于 指尖魔法屋-GPU计算：这次怎么落地的（https://blog.thinkmoon.cn/post/183-gpu-computing-cuda-optimization-guide/） 转载或引用必须申明原指尖魔法屋来源及源地址！

相关阅读

根据标签、分类和主题为你挑选的相近文章

GraphQL实战踩坑记录
这次搞 GraphQL，起因很直接：团队内部前端反馈查用户详情接口返回了 30 多个字段，实际只用了 5 个；另一个接口要三次请求才能拿到完整数…
📁 技术实践
📅 2021-05-28