使用下面一种或多种优化方法完成 CUDA 的矩阵乘法$C=A\times B$

• 使用 global memory 合并访存
• 采用分块乘法，使用 shared memory
• 请找出最佳的执行配置参数：grid 和 block

其中 $A$，$B$，$C$ 是$2^{12}\times 2^{12}$的方阵，按下面定义元素：

• a[i][j] = (i - 0.1 * j + 1) / (i + j + 1)
• b[i][j] = (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1)

实验环境

实验在老师提供的计算集群的一个节点上进行。单节点的显卡配置如下：

$ nvdia-smi
Mon Dec  2 08:38:49 2019
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 410.48                 Driver Version: 410.48                    |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|===============================+======================+======================|
|   0  Tesla V100-PCIE...  On   | 00000000:3B:00.0 Off |                    0 |
| N/A   30C    P0    24W / 250W |      0MiB / 16130MiB |      0%      Default |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                       GPU Memory |
|  GPU       PID   Type   Process name                             Usage      |
|=============================================================================|
|  No running processes found                                                 |
+-----------------------------------------------------------------------------+

实验原理

优化 CUDA 架构上的程序，一般从以下几个方面考虑：

• 选择好的并行算法，发掘更多的数据并行性
• 保持 SM 尽可能忙碌，尽量利用所有的 SM 参与计算
  • 加大数据量
  • 减小线程块大小
• 优化存储器的使用
  • 全局存储器合并访问
  • 使用更快的 constant memory 或 shared memory

实验过程

• 使用 global memory 合并访存
• 采用分块乘法，使用 shared memory
• 请找出最佳的执行配置参数：grid 和 blocks

这次实验和上一个实验CUDA 矩阵向量乘的多种优化其实非常相似：矩阵向量乘法中一个线程对应答案向量中的一个元素，矩阵矩阵乘法中一个线程对应答案矩阵的一个元素（使用二维分块）。不过，还有这些代码细节需要注意：

• 矩阵 B 的访问天然就是连续的，要使得对矩阵 A 的访问连续，可以考虑在 A 的列优先表达上进行计算；
• 矩阵 A 和 B 的每个位置都被访问了多次，因此都可以使用 shared memory 进行优化
• 由于测试的矩阵是长宽相等的正方形，因此在两个维度上的计算是大抵相似的，分块时也按正方形分块会有最少的访存（global 内存）次数
• 一个线程块里最多 1024 个线程，这里由于要正方形分块，因此分块的宽度不能超过 32（$32\times 32=1024$）
• 二维 Shared Memory 大小增加了一位用于避免列维度上多个线程访问同一个 Bank 产生 Bank conflict（不妨#define double float，更加明显，因为现在的卡的 bank 至少都有 32 个 banks，每个 bank 带宽是 32bit）

矩阵矩阵乘法的实现也可以通过重复调用之前的矩阵向量乘来实现，但是无法通过 shared memory 对矩阵的访存进行优化（只优化了对向量的访问），因而不如矩阵分块的方法优秀。此外，也可以通过调用 <cublas_v2.h> 或者 <cublasLt.h> 实现，但是不是本次实验重点。

template <size_t BLOCK_SIZE>
void __global__ MatMatMul(
	const double *Ac,
	const double *B,
	double *C,
	const size_t m,
	const size_t n,
	const size_t p)
{
	const size_t
		r = blockIdx.y * blockDim.y + threadIdx.y,
		c = blockIdx.x * blockDim.x + threadIdx.x;
	double res = 0;
	for (size_t t = 0; t < n; t += BLOCK_SIZE)
	{
		double __shared__
			sAc[BLOCK_SIZE + 1][BLOCK_SIZE + 1],
			sB[BLOCK_SIZE + 1][BLOCK_SIZE + 1];
		__syncthreads();
		sAc[threadIdx.y][threadIdx.x] = r < m && t + threadIdx.x < n ? Ac[(t + threadIdx.x) * m + r] : 0;
		sB[threadIdx.x][threadIdx.y] = c < p && t + threadIdx.y < n ? B[(t + threadIdx.y) * p + c] : 0;
		__syncthreads();
		for (size_t i = 0; i < blockDim.x; ++i)
			res += sAc[i][threadIdx.y] * sB[i][threadIdx.x];
	}
	if (r < m && c < p)
		C[r * p + c] = res;
}

最终测试结果如下：

可以看到，当执行参数设置成(128,128)grids, (32,32)blocks的时候，这里的矩阵乘法有最佳的性能表现。相对于分块大小为 1（相当于不分块），计算性能提高了将近一百六十倍，还是相当显著的。

源代码

MatMatMul.cu

#include <stdio.h>
#include <cuda_runtime.h>
template <size_t BLOCK_SIZE>
void __global__ MatMatMul(
	const double *Ac,
	const double *B,
	double *C,
	const size_t m,
	const size_t n,
	const size_t p)
{
	const size_t
		r = blockIdx.y * blockDim.y + threadIdx.y,
		c = blockIdx.x * blockDim.x + threadIdx.x;
	double res = 0;
	for (size_t t = 0; t < n; t += BLOCK_SIZE)
	{
		double __shared__
			sAc[BLOCK_SIZE + 1][BLOCK_SIZE + 1],
			sB[BLOCK_SIZE + 1][BLOCK_SIZE + 1];
		__syncthreads();
		sAc[threadIdx.y][threadIdx.x] = r < m && t + threadIdx.x < n ? Ac[(t + threadIdx.x) * m + r] : 0;
		sB[threadIdx.x][threadIdx.y] = c < p && t + threadIdx.y < n ? B[(t + threadIdx.y) * p + c] : 0;
		__syncthreads();
		for (size_t i = 0; i < blockDim.x; ++i)
			res += sAc[i][threadIdx.y] * sB[i][threadIdx.x];
	}
	if (r < m && c < p)
		C[r * p + c] = res;
}
int main()
{
	const size_t
		nRow = 1 << 12,
		nCol = 1 << 12;
	double
		*h_Ac,
		*h_B,
		*d_Ac,
		*d_B,
		*d_C;

	cudaMalloc(
		&d_Ac,
		sizeof(double) * nRow * nCol);
	cudaMalloc(
		&d_B,
		sizeof(double) * nCol * nRow);
	cudaMalloc(
		&d_C,
		sizeof(double) * nRow * nRow);

	cudaHostAlloc(
		&h_Ac,
		sizeof(double) * nRow * nCol,
		cudaHostAllocWriteCombined);
	cudaHostAlloc(
		&h_B,
		sizeof(double) * nCol * nRow,
		cudaHostAllocWriteCombined);

	for (size_t i = 0; i < nRow; ++i)
		for (size_t j = 0; j < nCol; ++j)
		{
			h_Ac[j * nRow + i] = (i - 0.1 * j + 1) / (i + j + 1);
			h_B[i * nCol + j] = (j - 0.2 * i + 1) * (i + j + 1) / (i * i + j * j + 1);
		}

	cudaMemcpy(
		d_Ac,
		h_Ac,
		sizeof(double) * nRow * nCol,
		cudaMemcpyHostToDevice);
	cudaMemcpy(
		d_B,
		h_B,
		sizeof(double) * nCol * nRow,
		cudaMemcpyHostToDevice);

	cudaFreeHost(h_Ac);
	cudaFreeHost(h_B);

	for (int tile_width = 32; tile_width; tile_width >>= 1)
	{
		cudaEvent_t beg, end;
		cudaEventCreate(&beg);
		cudaEventCreate(&end);
		cudaEventRecord(beg, 0);

		const dim3
			gridDim((nRow + tile_width - 1) / tile_width, (nRow + tile_width - 1) / tile_width),
			blockDim(tile_width, tile_width);
#define cal_kernel(TILE)        \
	do                          \
	{                           \
		if (tile_width == TILE) \
			MatMatMul<          \
				TILE><<<        \
				gridDim,        \
				blockDim>>>(    \
				d_Ac,           \
				d_B,            \
				d_C,            \
				nRow,           \
				nCol,           \
				nRow);          \
	} while (0)
		cal_kernel(32);
		cal_kernel(16);
		cal_kernel(8);
		cal_kernel(4);
		cal_kernel(2);
		cal_kernel(1);

		cudaDeviceSynchronize();
		cudaEventRecord(end, 0);
		cudaEventSynchronize(beg);
		cudaEventSynchronize(end);
		float elapsed_time;
		cudaEventElapsedTime(
			&elapsed_time,
			beg,
			end);
		printf(
			"(%d,%d)grids, (%d,%d)blocks, %fms elapsed.\n",
			gridDim.x,
			gridDim.y,
			blockDim.x,
			blockDim.y,
			elapsed_time);
	}

	cudaFree(d_Ac);
	cudaFree(d_B);
	cudaFree(d_C);
}

MatMatMul.o16434

各种参数和分块的运行时间。

(128,128)grids, (32,32)blocks, 108.868576ms elapsed.
(256,256)grids, (16,16)blocks, 125.305084ms elapsed.
(512,512)grids, (8,8)blocks, 178.711777ms elapsed.
(1024,1024)grids, (4,4)blocks, 615.897156ms elapsed.
(2048,2048)grids, (2,2)blocks, 2492.207520ms elapsed.
(4096,4096)grids, (1,1)blocks, 16997.095703ms elapsed.

MatMatMul.pbs

调度脚本。

#PBS -N MatMatMul
#PBS -l nodes=1:ppn=32:gpus=1
#PBS -j oe
#PBS -q gpu
source /public/software/profile.d/cuda10.0.sh
cd $PBS_O_WORKDIR
nvcc MatMatMul.cu -run