CUDA 编程入门笔记

入门一下 CUDA 编程，防止 remark

首先可以去 https://leetgpu.com/ 去做题，这是个 GPU kernel 学习、练习、评测和性能竞赛平台。

一般来说，

CPU 上执行的代码：host code
GPU 上执行的代码：device code

然后 __global__ 表示被 CPU 调用的 GPU 函数，__device__ 表示被 GPU 调用的 GPU 函数，__host__ 表示被 CPU 调用的 CPU 函数。

先来看一道向量加法的题目：https://leetgpu.com/challenges/vector-addition

然后直接来看 CUDA 代码：

#include <cuda_runtime.h>

__global__ void vector_add(
    const float* A,
    const float* B,
    float* C,
    int N
) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

// A, B, C are device pointers
extern "C" void solve(
    const float* A,
    const float* B,
    float* C,
    int N
) {
    int threadsPerBlock = 256;
    int blocksPerGrid =
        (N + threadsPerBlock - 1) / threadsPerBlock;

    vector_add<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, N);

    cudaDeviceSynchronize();
}

solve 就是 CPU 会调用的函数，vector_add 是运行在 GPU 上，被 CPU 调用的 kernal，核函数。这个函数会被并行执行。

每次 kernel 启动（也就是说调用 fff<<<blocksPerGrid, threadsPerBlock>>>(...)）会产生一个 Grid，而线程的组织结构是一个 Grid 下面有多个 Block，一个 Block 下面有多个 Threshold（线程），上面代码里的 blockDim.x 表示一个 Block 有多少个线程，blockIdx.x 表示当前线程在第几个 Block 里，threadIdx.x 表示当前线程在所在 Block 里的编号。这里计算 i（原数组下标）的方式其实就和二维数组 a[][] 求 a[i][j] 相对 a 的内存地址偏移量很相似，blockDim.x 对应 a 的第二维的大小，blockIdx.x 对应 i，threadIdx.x 对应 j。这里可以看出每个线程会计算一次标量加法。

当然你想让一个线程做多个操作也是可以的，这有个好处就是让 kernel 能适应不同规模的输入，如下所示：

__global__ void vector_add(
    const float* A,
    const float* B,
    float* C,
    int N
) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    int stride = blockDim.x * gridDim.x;

    for (; i < N; i += stride) {
        C[i] = A[i] + B[i];
    }
}

但是这道向量加法的题目每次只做一次加法可以让并行度最大，负载均衡，且“内存访问连续（合并访存）”，具体还不太清楚，后面会写一下。如果一个线程做太多操作的话，还会大量占用寄存器，使得同时能并行的线程数量下降，这是不好的。

再解释代码里两个问题：

1. cudaDeviceSynchronize(); 是干什么的？

这是表示让 CPU 等待，直到 GPU 当前任务执行完毕。

2. extern "C" 是干什么的？

这是 C++ 的链接规则，和 GPU 没关系，它要求编译器使用 C 风格函数名，防止编译时发生 name mangling，函数名字变复杂，以免评测系统无法稳定地找到名为 solve 的函数。

再来看一道矩阵乘法的代码：https://leetgpu.com/challenges/matrix-multiplication

#include <cuda_runtime.h>

__global__ void matrix_multiplication_kernel(const float* A, const float* B, float* C, int M, int N, int K) {
    int i = blockIdx.x * blockDim.x + threadIdx.x, j = blockIdx.y * blockDim.y + threadIdx.y;
    if (i < M && j < K) {
        float sum = 0.0f;
        for (int k = 0; k < N; ++k) {
            sum += A[i * N + k] * B[k * K + j];
        }
        C[i * K + j] = sum;
    }
}

// A, B, C are device pointers (i.e. pointers to memory on the GPU)
extern "C" void solve(const float* A, const float* B, float* C, int M, int N, int K) {
    dim3 threadsPerBlock(16, 16);
    dim3 blocksPerGrid((M + threadsPerBlock.x - 1) / threadsPerBlock.x,
                       (K + threadsPerBlock.y - 1) / threadsPerBlock.y);

    matrix_multiplication_kernel<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, M, N, K);
    cudaDeviceSynchronize();
}

dim3 是 CUDA 提供的一个表示三维尺寸的类型，包含：.x .y .z，这里传两个参数相当于 .x = 16, .y = 16, .z = 1。前面一道题，用两个 int 去启动 vector_add，就相当于只用到第一维的 dim3。（CUDA 不支持更高维）

然后这个代码看上去没什么问题，但是交上去你会发现 TLE 了，为什么呢？因为 CUDA 在线程块中会优先把 threadIdx.x 连续的线程组成 warp。也就是说，同一个 warp 中通常是 threadIdx.y 相同，而 threadIdx.x 单调递增，跟 C++ 的数组在内存里的储存顺序是相反的。所以，为了让相邻线程访问相邻的全局内存地址，合并显存访问，这里需要将 i，j 以及 blocksPerGrid 的两维反向一下。GPU 编程里面往往更关心线程间的访问连续性而不是线程内的访问连续性。

改成下面的代码就可以通过了：

#include <cuda_runtime.h>

__global__ void matrix_multiplication_kernel(const float* A, const float* B, float* C, int M, int N, int K) {
    int i = blockIdx.y * blockDim.y + threadIdx.y, j = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < M && j < K) {
        float sum = 0.0f;
        for (int k = 0; k < N; ++k) {
            sum += A[i * N + k] * B[k * K + j];
        }
        C[i * K + j] = sum;
    }
}

// A, B, C are device pointers (i.e. pointers to memory on the GPU)
extern "C" void solve(const float* A, const float* B, float* C, int M, int N, int K) {
    dim3 threadsPerBlock(16, 16);
    dim3 blocksPerGrid((K + threadsPerBlock.x - 1) / threadsPerBlock.x,
                       (M + threadsPerBlock.y - 1) / threadsPerBlock.y);

    matrix_multiplication_kernel<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, M, N, K);
    cudaDeviceSynchronize();
}