转载-CUDA矩阵乘法的优化
写在前面
本文转载自吴坎的博客。
实验简介
使用下面一种或多种优化方法完成 CUDA 的矩阵乘法\(A\times B=C\)
- 使用 global memory 合并访存
- 采用分块乘法，使用 shared memory
- 请找出最佳的执行配置参数：grid 和 blocks
其中 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 对矩阵的访存进行优化（只优化了对向量的访问），因而不如矩阵分块的方法优秀。
void __global__ MatMatMul(lf *At, lf *B, lf *C, size_t m, size_t n, size_t p)
{
extern lf __shared__ shared[];
lf res = 0, *Ats = shared, *Bs = shared + blockDim.x * blockDim.y;
size_t
r = blockIdx.y * blockDim.x + threadIdx.y,
c = blockIdx.x * blockDim.x + threadIdx.x;
for (size_t t = 0; t < n; t += blockDim.x)
{
__syncthreads();
Ats[threadIdx.y * blockDim.x + threadIdx.x] = r < m && t + threadIdx.x < n ? At[(t + threadIdx.x) * m + r] : 0;
Bs[threadIdx.x * blockDim.y + 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 += Ats[i * blockDim.y + threadIdx.y] * Bs[i * blockDim.x + threadIdx.x];
}
if (r < m && c < p)
C[r * p + c] = res;
}
最终测试结果如下：
| Grids | Blocks | Elapsed Time(ms) |
|---|---|---|
| (128,128) | (32,32) | 102.511581 |
| (256,256) | (16,16) | 96.044540 |
| (512,512) | (8,8) | 109.401154 |
| (1024,1024) | (4,4) | 385.759155 |
| (2048,2048) | (2,2) | 2366.376953 |
| (4096,4096) | (1,1) | 16237.383789 |
可以看到，当执行参数设置成(256,256)grids, (16,16)blocks的时候，这里的矩阵乘法有最佳的性能表现。相对于分块大小为 1（相当于不分块），计算性能提高了一百六十多倍，还是相当显著的。
源代码
MatMatMul.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
typedef float lf;
void __global__ MatMatMul(lf *At, lf *B, lf *C, size_t m, size_t n, size_t p)
{
extern lf __shared__ shared[];
lf res = 0, *Ats = shared, *Bs = shared + blockDim.x * blockDim.y;
size_t
r = blockIdx.y * blockDim.x + threadIdx.y,
c = blockIdx.x * blockDim.x + threadIdx.x;
for (size_t t = 0; t < n; t += blockDim.x)
{
__syncthreads();
Ats[threadIdx.y * blockDim.x + threadIdx.x] = r < m && t + threadIdx.x < n ? At[(t + threadIdx.x) * m + r] : 0;
Bs[threadIdx.x * blockDim.y + 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 += Ats[i * blockDim.y + threadIdx.y] * Bs[i * blockDim.x + threadIdx.x];
}
if (r < m && c < p)
C[r * p + c] = res;
}
int main()
{
const size_t
nRow = 1 << 12,
nCol = 1 << 12;
lf
*h_A = (lf *)malloc(sizeof(lf) * nRow * nCol),
*h_At = (lf *)malloc(sizeof(lf) * nRow * nCol),
*h_B = (lf *)malloc(sizeof(lf) * nCol * nRow),
*d_A,
*d_At,
*d_B,
*d_C;
for (size_t i = 0; i < nRow; ++i)
for (size_t j = 0; j < nCol; ++j)
{
h_A[i * nCol + j] = (i - 0.1 * j + 1) / (i + j + 1);
h_At[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);
}
cudaMalloc((lf **)&d_A, sizeof(lf) * nRow * nCol);
cudaMalloc((lf **)&d_At, sizeof(lf) * nRow * nCol);
cudaMalloc((lf **)&d_B, sizeof(lf) * nCol * nRow);
cudaMalloc((lf **)&d_C, sizeof(lf) * nRow * nRow);
cudaMemcpy(d_A, h_A, sizeof(lf) * nRow * nCol, cudaMemcpyHostToDevice);
cudaMemcpy(d_At, h_At, sizeof(lf) * nRow * nCol, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, sizeof(lf) * nCol, cudaMemcpyHostToDevice);
for (int tile_width = 32; tile_width; tile_width >>= 1)
{
cudaEvent_t beg, end;
cudaEventCreate(&beg);
cudaEventCreate(&end);
cudaEventRecord(beg, 0);
dim3
grids((nRow - 1) / tile_width + 1, (nRow - 1) / tile_width + 1),
blocks(tile_width, tile_width);
MatMatMul<<<grids, blocks, sizeof(lf) * 2 * blocks.x * blocks.y>>>(d_At, d_B, d_C, nRow, nCol, nRow);
cudaDeviceSynchronize();
cudaEventRecord(end, 0);
cudaEventSynchronize(beg);
cudaEventSynchronize(end);
lf elapsed_time;
cudaEventElapsedTime(&elapsed_time, beg, end);
printf("(%d,%d)grids, (%d,%d)blocks, %fms elapsed.\n", grids.x, grids.y, blocks.x, blocks.y, elapsed_time);
}
cudaFree(d_A);
cudaFree(d_At);
cudaFree(d_B);
cudaFree(d_C);
free(h_A);
free(h_At);
free(h_B);
}
MatMatMul.o16434
各种参数和分块的运行时间。
(128,128)grids, (32,32)blocks, 102.511581ms elapsed.
(256,256)grids, (16,16)blocks, 96.044540ms elapsed.
(512,512)grids, (8,8)blocks, 109.401154ms elapsed.
(1024,1024)grids, (4,4)blocks, 385.759155ms elapsed.
(2048,2048)grids, (2,2)blocks, 2366.376953ms elapsed.
(4096,4096)grids, (1,1)blocks, 16237.383789ms 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 -o MatMatMul
./MatMatMul
浙公网安备 33010602011771号