本文记录了我的cuda学习经历，和大多数人一样，通过优化矩阵乘法的过程来了解一些基本的概念。仓库链接：

Refences

其中Fermi架构是Compute Capability 2.0的架构。从白皮书里能了解到硬件相关的一些基本概念。比如streaming multiprocessor，有时候也简称multiprocessor或者SM。一个SM里有32个cuda core，有两个warp调度器。一个warp是由32个thread组成。和硬件结合后就比较容易理解，为什么一个block里最好至少放64个thread，因为有两个warp scheduler存在，至少可以放两个warp的thread进行工作。

Programming Guide里比较详细地介绍了编程模型（Programming Model），也比较详细地介绍了一些Runtime API。CUDA也提供了更底层的Driver API，但一般Runtime API已经够用了，而且使用起来更容易。除此之外，我看的比较多的还有不同版本的Compute Capability的介绍，其中包括每个SM最多能同时处理的block数量，每个block最大的线程数……这些在实际调用kernel的时候都需要考虑。

Best Practices Guide介绍的优化技巧和硬件就比较相关了。特别是需要了解设备中的存储结构，因为很大部分情况下是在想办法降低访存的延时。比如Memory Optimizations这一章介绍的内容就非常值得细看。

System Requirment

  • Ubuntu 20.04
  • NVIDIA Driver Version 550.67
  • CUDA Version 12.4
  • Eigen repository-url
  • cmake version 3.25.1
  • gcc (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0

How To Build

在开始之前需要有一台带有Nvidia显卡的主机，然后安装上驱动，最新的驱动可以从官网下载得到。直接运行然后按照指示进行安装即可，网上的教程需要用户手动去禁用nouveau，但现在这些操作驱动安装程序都可以完成，所以不需要额外的准备工作了（至少我安装的时候是这样的）。

驱动安装完成后再同样按照官网的指导步骤安装CUDA Toolkit。安装完成后再修改.zshrc或.bashrc将bin路径和lib路径添加到分别添加到PATH和LD_LIBRARY_PATH中。

1
2
export PATH=/usr/local/cuda/bin${PATH:+:${PATH}}
export LD_LIBRARY_PATH=/usr/local/cuda/lib64${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}

另外本项目依赖Eigen，因此还需要另外安装Eigen库。拉取源码后直接install即可。

Check Compute Capability

在开始前需要确定显卡的Compute Capability，可通过官网查询。也可以先编译check_cc来获取当前设备的Compute Capability。（默认查询的是device 0 的设备信息）

1
2
3
4
mkdir build
cd build
cmake ..
make check_cc

我的设备用到了两种显卡，因此设置的CMAKE_CUDA_ARCHITECTURES是75和86，如果你的设备用的是其他的显卡，可以修改CMakeLists.txt，将CMAKE_CUDA_ARCHITECTURES修改成你需要的值。

Build All

切换到build目录下，make all即可。

Content

本项目做的主要就是用Nvidia GPU实现两个NxN的双精度矩阵乘法，在common.h中，设置了N的大小以及thread block的大小。

主要内容：

  1. baseline - 最基础的矩阵乘法实现方式与eigen，cublas的实现进行对比；
  2. shared_memory - 用共享内存实现，减少访问global memory的次数；
  3. coalesce - 尽可能用coalesce的形式访存；减少访问shared memory的bank冲突；
  4. other_practice - 用capture graph执行矩阵乘法；用memory mapped方式执行矩阵乘法。

baseline

按照矩阵乘法的规则，两个NxN的矩阵相乘，得到的也是一个NxN的矩阵。结果矩阵中的每个元素都是由一个行向量和一个列向量求内积得到的。最直接的想法就是用NxN个线程来完成计算，每个线程负责计算一组内积。

1
2
3
4
5
6
7
8
9
10
__global__ void basic(int N, double *a, double *b, double *c) {
  int row = blockIdx.y * blockDim.y + threadIdx.y;
  int col = blockIdx.x * blockDim.x + threadIdx.x;
  double sum = 0.0;
  for (int i = 0; i < N; i++) {
    // 注意数据是按照列优先存储的
    sum += a[row + i * N] * b[col * N + i];
  }
  c[col * N + row] = sum;
}

值得注意的是，如果使用Eigen库，那么矩阵的数据是优先按照列存放的，即矩阵中同一列的数据是连续地址存放的。

1
2
3
4
5
6
7
8
9
// 进行矩阵乘法。将数据拷贝到设备，再将结果拷贝回来
Eigen::MatrixXd result_cuda = Eigen::MatrixXd::Zero(N, N);
start = std::chrono::high_resolution_clock::now();
cudaMemcpy(d_mat1, mat1.data(), N * N * sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(d_mat2, mat2.data(), N * N * sizeof(double), cudaMemcpyHostToDevice);
basic<<<gridSize, blockSize>>>(N, d_mat1, d_mat2, d_result);
cudaMemcpy(result_cuda.data(), d_result, N * N * sizeof(double), cudaMemcpyDeviceToHost);
cudaDeviceSynchronize();
end = std::chrono::high_resolution_clock::now();

调用自己写的kernel实现矩阵乘法也不需要多考虑什么，直接在一条stream上，先把数据拷贝到device，再进行运算，算完之后再把结果拷贝到host。统计这整个过程的时间。

分别将Eigen库中的矩阵乘法，我们自己写的基础的矩阵乘法和cublas库实现的矩阵乘法进行对比。因为用到里Eigen库，所以编译的时候一定要编译Release的版本，不然用Eigen库实现的矩阵乘法需要很久。不光是比较时间，而且也需要确认计算结果的正确性。以Eiegn库的计算结果作为参照，要求GPU的计算结果与Eigen库的计算结果相同。最终得到如下结果：

再用Nsight Systems分析一下，可以看到这个自己写的kernel相比于cublas实现的kernel真的慢很多。

shared memory

在device上malloc的数据分配在global memory中，访问global memory相对来说是比较慢的，而访问shared memory会很快，shared memory类似于scratchpad memory，是一块可以由程序员自己管理的cache。

shared memory是在片内的，而且不会经过cache，因为它很小而且访问速度够快；global memory是在片外的，而且会经过L1，L2cache。

所以如果频繁地从global memory读数据是很费时间的。考虑矩阵乘法的过程，两个相乘的矩阵需要分别被读取N次才能计算得到最终结果。

shared memory是一个thread block中的线程共享的，那么就可以考虑让一个thread block中的线程“互帮互助”。一个block中共有BLOCK_SIZExBLOCK_SIZE个线程，每个线程从global memory拷贝一个数据到到shared memory，这些数据就可以由这个block中的线程共享。在一个block中，原本一个线程需要从global memory读BLOCK_SIZE个数据，而采用共享的方式之后，就可以每个block里每个线程只读1个数据。完整的计算下来，就只需要从global memory读N/BLOCK_SIZE次数据。

这样看似乎BLOCK_SIZE越大，加速效果会越明显。但实际上block越大，block中的线程同步也会更久。而且一个block中的shared memeory是有限的，register也是有限的。一般来说总的线程数一定的话，block分的小一点，多一点，更容易把每个block分到multiprocessor上去运行，提高并行度。所以BLOCK_SIZE的选取我没有去详细比较分析了，为了简便起见，我在这里固定用的是16x16的大小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
__global__ void shared_memory(double *a, double *b, double *c) {
  int result_row = blockIdx.y * BLOCK_SIZE + threadIdx.y;
  int result_col = blockIdx.x * BLOCK_SIZE + threadIdx.x;

  // 将结果清空
  c[result_row + result_col * N] = 0.0;

  int a_col_global, b_row_global;

  // 每个block一起load数据，放入s_a s_b中
  __shared__ double s_a[BLOCK_SIZE][BLOCK_SIZE];
  __shared__ double s_b[BLOCK_SIZE][BLOCK_SIZE];

  // 每个thread需要load N/BLOCK_SIZE次数据
  for (int i = 0; i < N / BLOCK_SIZE; ++i) {
    // 计算要搬运的数据在global下的索引
    a_col_global = i * BLOCK_SIZE + threadIdx.x;
    b_row_global = i * BLOCK_SIZE + threadIdx.y;

    // 搬运数据
    s_a[threadIdx.y][threadIdx.x] = a[result_row + a_col_global * N];
    s_b[threadIdx.y][threadIdx.x] = b[b_row_global + result_col * N];
    __syncthreads();

    // 计算部分和
    for (int j = 0; j < BLOCK_SIZE; ++j) {
      c[result_row + result_col * N] += s_a[threadIdx.y][j] * s_b[j][threadIdx.x];
    }
    __syncthreads();
  }
}

从运行结果看，用shared memory之后，乘法计算所需的时间明显降低。

查看nsight systems的分析结果：

coalesce

对比采用shared memory加速的结果和cublas实现的结果，实际上还是有很多优化空间。

有个新的概念是coalesced memory access，是为了最大化数据传输的带宽，尽可能使数据“联合”访问。

一个warp中的线程如果是访问global memory中的一块连续地址，那就是可以联合访问的。block的维度有一维，二维，三维，这是为了更好地与具体应用进行映射，而block的性能，只和block的size有关。也就是说8x2的block和4x4的block本质上是一样的，都是16个线程。每个线程有自己的ID，类似于二维数组，二维block中的线程ID计算：id=x+yDxid=x+yD_x，其中xx, yy分别是两个维度的索引值，DxD_x表示x方向的维度。然后以连续的32个ID的thread作为一个warp。在设计kernel的时候需要考虑让一个warp中的thread访问连续的一片global memory。

shared memory分32个bank，每个bank是4字节，比如字节地址0~3属于bank0，4~7属于bank1，8~11属于bank2……不同bank的数据可以同时被访问，同一bank的数据就不能一起访问。比如字节地址0和字节地址128的数据都属于bank0，就不能一起访问。

在矩阵乘法过程中，从global memory读数据，然后写入shared memory。不仅需要考虑从global memory的连续地址读取数据，而且在写入shared memory的时候也需要考虑减少bank冲突。

考虑到矩阵乘法axb的时候，a中取一行数据，和b中的一列数据计算内积。b本身就是按列存放数据的，warp中的线程也按照列来组织，这样能够保证warp中的线程访问的是shared memory中的连续的一片数据。但是要取a中的一行数据，就要跨行取数了，或者是当我从global memory读数据后，就将矩阵转置，再存到shared memory里。这里不管怎样都会涉及跨行的问题，而且如果block size是128字节的整数倍，那就肯定会有bank访问冲突。这里有个技巧就是在shared memory中额外多分配一点空间，从而让跨行的数据不再是128字节的整数倍，人为地让地址错开。我最后采用的做法就是将数据转置后存入shared memory，而在计算部分和的时候每个warp就可以从连续的地址加载数据了。在矩阵的维度足够大的时候，bank冲突是无法避免的，只不过通过这种方式能够充分利用访问shared memory的带宽，减少无效的数据访问。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
__global__ void coalesce(double *a, double *b, double *c) {
  // 数据按列存放，所以x按列方向增长，y按行方向增长
  int result_row = blockIdx.x * BLOCK_SIZE + threadIdx.x;
  int result_col = blockIdx.y * BLOCK_SIZE + threadIdx.y;

  // 将结果清空
  c[result_col * N + result_row] = 0.0;

  // 每个block一起load数据，放入s_a s_b中，同一列的数据在地址上不能对齐，对齐的话会bank访问冲突
  __shared__ double s_a[BLOCK_SIZE][BLOCK_SIZE + 1];
  __shared__ double s_b[BLOCK_SIZE][BLOCK_SIZE + 1];

  // 每个thread需要load N/BLOCK_SIZE次数据
  for (int i = 0; i < N / BLOCK_SIZE; ++i) {
    // (y, x) thread 负责拷贝global中的 (y, x)的数据
    int s_a_row_in_global = result_row;
    int s_a_col_in_global = i * BLOCK_SIZE + threadIdx.y;
    int s_b_col_in_global = result_col;
    int s_b_row_in_global = i * BLOCK_SIZE + threadIdx.x;

    // s_a中的数据需要跨行存放，因为一个warp读的是列数据
    s_a[threadIdx.x][threadIdx.y] = a[s_a_row_in_global + s_a_col_in_global * N];
    s_b[threadIdx.y][threadIdx.x] = b[s_b_row_in_global + s_b_col_in_global * N];
    __syncthreads();

    for (int j = 0; j < BLOCK_SIZE; ++j) {
      c[result_col * N + result_row] += s_a[threadIdx.x][j] * s_b[threadIdx.y][j];
    }
    __syncthreads();
  }
}

从运行结果可以看出，在考虑上访存的过程后，同样维度的矩阵乘法进一步得到加速，而且与cublas实现的性能比较接近了。

查看nsight systems的分析结果：

other practice

除了一些优化的措施外我也做了一些其它的尝试。

graph

根据官方文档里介绍的，可以将一些要在某条stream上执行的任务capture，来组成一个capture graph。再结合event，graph不再是单条stream的顺序结构，而是可以结合多个stream，形成一个有向图。在graph的节点中还可以添加条件判断，做一些简单的控制逻辑。graph在执行的过程中虽然也是stream，但它在实例化的过程中可以提前完成一些工作，而且经过实例化后可以多次执行，相比于单纯的stream应该会更快。

所以我也实际用graph试了一下，实际效果并不会有明显的提速，但我觉得这个在其他地方应该会有应用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// 创建capture graph
cudaStreamBeginCapture(s, cudaStreamCaptureModeRelaxed);
cudaMemcpyAsync(d_mat1, mat1.data(), N * N * sizeof(double), cudaMemcpyHostToDevice, s);
cudaMemcpyAsync(d_mat2, mat2.data(), N * N * sizeof(double), cudaMemcpyHostToDevice, s);
coalesce<<<gridSize, blockSize, 0, s>>>(d_mat1, d_mat2, d_result);
cudaMemcpyAsync(result_graph.data(), d_result, N * N * sizeof(double), cudaMemcpyDeviceToHost, s);
cudaStreamEndCapture(s, &graph);

// graph实例化
cudaGraphInstantiate(&exec, graph, NULL, NULL, 0);

// 运行graph
start = std::chrono::high_resolution_clock::now();
cudaGraphLaunch(exec, s);
cudaDeviceSynchronize();
end = std::chrono::high_resolution_clock::now();

memory mapped

在矩阵乘法的例子里，数据在device和host之间传输的时间开销并不明显。但我也在想能不能省去这部分拷贝数据过程，然后看到了有一个地址映射的操作。可以将host的一块数据映射到device端，相当于让kernel直接处理host的memory数据。

1
2
3
4
5
6
7
8
9
10
11
12
13
// 分配结果存放的空间，获取map后的device端地址
Eigen::MatrixXd reuslt_mapped = Eigen::MatrixXd::Zero(N, N);
double *d_mat1_mapped, *d_mat2_mapped, *d_result_mapped;
cudaHostRegister(reuslt_mapped.data(), N * N * sizeof(double), cudaHostRegisterDefault);
cudaHostGetDevicePointer(&d_mat1_mapped, mat1.data(), 0);
cudaHostGetDevicePointer(&d_mat2_mapped, mat2.data(), 0);
cudaHostGetDevicePointer(&d_result_mapped, reuslt_mapped.data(), 0);

// 进行矩阵乘法
start = std::chrono::high_resolution_clock::now();
coalesce<<<gridSize, blockSize>>>(d_mat1_mapped, d_mat2_mapped, d_result_mapped);
cudaDeviceSynchronize();
end = std::chrono::high_resolution_clock::now();

实际运行后发现这么做会更慢，需要尽量减少host和device之间的数据拷贝！