测试案例
优化方案
- 使用 global memory
- 使用合并访存
- 使用 constant memory 存放向量
- 使用 shared memory 存放向量和矩阵
测试环境
基本原理
优化 CUDA 架构上的程序，一般从以下几个方面考虑：
- 选择好的并行算法，发掘更多的数据并行性
- 保持 SM 尽可能忙碌，尽量利用所有的 SM 参与计算
- 加大数据量
- 减小线程块大小
- 优化存储器的使用
- 全局存储器合并访问
- 使用更快的 constant memory 或 shared memory
方案对比测试
由于都是 CUDA 架构上的核函数对比性能，下面的计时都只测了用于核函数计算的时间，而不包含数据拷贝的部分（否则运行时间都在 300ms 左右，基本上都是拷贝的时间而没有参考价值了）。当然，由于没有计入拷贝等预处理的时间，那些需要计算转置或者预读取的算法在这里会有优势一些。
1.使用 global memory
这是最基础的矩阵向量乘法。这里假设线程块都是一维组织的，每个 CUDA 线程计算矩阵的一行与向量乘积，这样各线程之间没有读写冲突，不需要使用原子操作。
void __global__ MatVecMulGlobalMemory(const lf *A, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0; //将结果先存在寄存器里，减少对向量C的访存
for (size_t j = 0; j < nCol; ++j)
res += A[i * nCol + j] * B[j];
C[i] = res;
}
}2.使用合并访存
所谓合并访存，指的是相邻的线程访问段对齐的地址。比如在之前的代码中，j == 0时线程 0 访问A[0]，线程 1 访问A[nCol]，线程 2 访问A[2 * nCol]…它们并不相邻，因此不满足合并访问的要求。在这里我们把原来的矩阵A𝐴求出转置AT𝐴𝑇，此时j == 0时线程 0 访问At[0]，线程 1 访问At[1]，线程 2 访问At[2]…此时满足了合并访问的要求。
void __global__ MatVecMulGlobalMemoryAlign(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0;
for (size_t j = 0; j < nCol; ++j)
res += At[j * nRow + i] * B[j];
C[i] = res;
}
}3.使用 constant memory 存放向量
注意到向量在计算过程中不会改变，且每个线程访问相同地址，因此考虑把它放在 constant memory 中。
NVIDIA 硬件提供了 64KB 的常量内存，并且常量内存采用了不同于标准全局内存的处理方式。在这里我们大小为214214的单精度浮点数向量b𝑏大小恰好为 64KB，正好可以完整保存。如果向量超过了 constant memory 的 64KB 上限，那就需要分批进行，多次传输和启动内核。
lf __constant__ d_Bc[(1 << 16) / sizeof(lf)]; //64KB
void __global__ MatVecMulConstantMemory(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0;
for (size_t j = 0; j < nCol; ++j)
res += At[j * nRow + i] * d_Bc[j];
C[i] = res;
}
}
运行时间为1.536000ms，在上一步的基础上略微提高。使用常量内存可以提升运算性能的原因主要有两个：
- 对常量内存的单次读操作可以广播到同个半线程束的其他1515个线程，这种方式产生的内存流量只是使用全局内存时的116116。
- 硬件将主动把常量数据缓存在 GPU 上。在第一次从常量内存的某个地址上读取后，当其他半线程束请求同一个地址时，那么将命中缓存，这同样减少了额外的内存流量。
4.使用 shared memory 存放向量和矩阵
对于 block 内内存来说，向量都是共享的，因此我们可以使用比 constant memory 更快的 shared memory 来存储，此时相比较使用常量内存，我们免掉了向量比较大的时候多次数据拷贝和启动核函数的开销，也没有使用全局变量，增加了代码的可扩展性。当然，shared memory 更小，因此需要对向量进行分块处理。
另外需要更正的一个问题是，并不需要使用 shared memory 去存矩阵，因为在这个矩阵向量乘的过程中，每个矩阵元素只被访问了一次。此外，shared memory 的大小也并不足以存下完整的矩阵（甚至是向量）。
void __global__ MatVecMulSharedMemory(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
extern lf __shared__ Bs[];
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
lf res = 0;
for (size_t jBeg = 0, jEnd = blockDim.x < nCol ? blockDim.x : nCol;
jBeg < nCol;
jBeg += blockDim.x, jEnd += blockDim.x)
{
__syncthreads(); //防止有的进程还在读Bs
if (jBeg + threadIdx.x < nCol)
Bs[threadIdx.x] = B[jBeg + threadIdx.x];
__syncthreads();
if (i < nRow)
for (size_t j = jBeg; j < jEnd; ++j)
res += At[j * nRow + i] * Bs[j - jBeg];
}
if (i < nRow)
C[i] = res;
}
运行时间为1.821696ms，反而有一定的下降。分析一下原因：
- 使用常量内存时，将原来的数据拷贝到常量内存的时间并没有被算进去，而这里读内存的过程是在核函数内部的。想来如果在更大的矩阵上进行计算，使用 shared memory 应该会有更好的表现。
- 老师集群上的显卡性能过于强悍（在今年十一月 SC 超算大会刚发布 Tesla V100S 前，Tesla V100 一直都是市面能买到的最强算力），内存读写性能比以往的显卡都要强很多，因此对本来已经很快的 global memory 的优化效果没有那么明显了，而对应的核函数却比朴素的算法更复杂，一定程度上增加了运行的时间常数。
5.使用shuffle
我试图使用 shuffle 去优化，但是结果不尽人意。原因是，__shfl_sync的广播操作耗时出人意料地高，而__shfl_xor_sync倒是避开了广播，但是又导致了At的访存不连续。
测试结果
五个方案耗时分别为：
3.728736ms
1.580352ms
1.521280ms
1.773952ms
1.981728ms
可以看到，由于现在的硬件性能已经大大强于数年前，做存储器的优化效果已经比较小（并不是说没有）。因此，对这个问题来说，最主要的是要选一个优秀的并行算法，并对程序代码做好访存分析和优化。
当然也不是说存储器结构就不再重要，还是要具体问题具体分析。上面很多算法都是要对矩阵或者向量进行预处理的，而并没有把对应的代价（时间、内存空间、可扩展性等）计入在内，实际上在运用到生产环境的时候这些仍然是必须要考虑的。
完整代码
MatVecMul.cu
#include <stdio.h>
#include <math.h>
#include <cuda_runtime.h>
typedef float lf;
void __global__ MatVecMulGlobalMemory(const lf *A, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0; //将结果先存在寄存器里，减少对向量C的访存
for (size_t j = 0; j < nCol; ++j)
res += A[i * nCol + j] * B[j];
C[i] = res;
}
}
void __global__ MatVecMulGlobalMemoryAlign(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0;
for (size_t j = 0; j < nCol; ++j)
res += At[j * nRow + i] * B[j];
C[i] = res;
}
}
lf __constant__ d_Bc[(1 << 16) / sizeof(lf)]; //64KB
void __global__ MatVecMulConstantMemory(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < nRow)
{
lf res = 0;
for (size_t j = 0; j < nCol; ++j)
res += At[j * nRow + i] * d_Bc[j];
C[i] = res;
}
}
void __global__ MatVecMulSharedMemory(const lf *At, const lf *B, lf *C, const size_t nRow, const size_t nCol)
{
extern lf __shared__ Bs[];
const size_t i = blockDim.x * blockIdx.x + threadIdx.x;
lf res = 0;
for (size_t jBeg = 0, jEnd = blockDim.x < nCol ? blockDim.x : nCol;
jBeg < nCol;
jBeg += blockDim.x, jEnd += blockDim.x)
{
__syncthreads(); //防止有的进程还在读Bs
if (jBeg + threadIdx.x < nCol)
Bs[threadIdx.x] = B[jBeg + threadIdx.x];
__syncthreads();
if (i < nRow)
for (size_t j = jBeg; j < jEnd; ++j)
res += At[j * nRow + i] * Bs[j - jBeg];
}
if (i < nRow)
C[i] = res;
}
void __global__ MatVecMulWarp(
const lf *At,
const lf *B,
lf *C,
const size_t nRow,
const size_t nCol)
{
const size_t
i = blockDim.x * blockIdx.x + threadIdx.x,
warp_size = 32,
lane_id = i % warp_size;
lf res = 0;
for (size_t jBeg = 0, jEnd = warp_size < nCol ? warp_size : nCol;
jBeg < nCol;
jBeg += warp_size, jEnd += warp_size)
{
const lf val = jBeg + lane_id < nCol ? B[jBeg + lane_id] : 0;
for (size_t j = 0; j < warp_size; ++j)
res += (i < nRow ? At[(jBeg + (lane_id ^ j)) * nRow + i] : 0) * __shfl_xor_sync(0xffffffff, val, j);
}
if (i < nRow)
C[i] = res;
}
int main()
{
const size_t
nRow = 1 << 14,
nCol = 1 << 14;
lf
*h_A,
*h_At,
*h_B,
*h_C,
*d_A,
*d_At,
*d_B,
*d_C;
cudaHostAlloc(
(void **)&h_A,
sizeof(lf) * nRow * nCol,
cudaHostAllocWriteCombined);
cudaHostAlloc(
(void **)&h_At,
sizeof(lf) * nRow * nCol,
cudaHostAllocWriteCombined);
cudaHostAlloc(
(void **)&h_B,
sizeof(lf) * nCol,
cudaHostAllocWriteCombined);
cudaHostAlloc(
(void **)&h_C,
sizeof(lf) * nRow,
cudaHostAllocWriteCombined);
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;
for (size_t j = 0; j < nCol; ++j)
{
h_B[j] = log(sqrt(j * j - j + 2));
for (size_t i = 0; i < nRow; ++i)
h_At[j * nRow + i] = i - 0.1 * 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);
cudaMalloc((lf **)&d_C, sizeof(lf) * 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);
cudaMemcpyToSymbol(d_Bc, h_B, sizeof(lf) * nCol, cudaMemcpyHostToDevice);
cudaFreeHost(h_A);
cudaFreeHost(h_At);
cudaFreeHost(h_B);
cudaFreeHost(h_C);
for (int i = 0; i < 5; ++i)
{
cudaEvent_t beg, end;
cudaEventCreate(&beg);
cudaEventCreate(&end);
cudaEventRecord(beg, 0);
size_t blocks = 1 << 7, grids = (nRow + blocks - 1) / blocks;
if (i == 0)
MatVecMulGlobalMemory<<<grids, blocks>>>(d_A, d_B, d_C, nRow, nCol);
else if (i == 1)
MatVecMulGlobalMemoryAlign<<<grids, blocks>>>(d_At, d_B, d_C, nRow, nCol);
else if (i == 2)
MatVecMulConstantMemory<<<grids, blocks>>>(d_At, d_B, d_C, nRow, nCol);
else if (i == 3)
MatVecMulSharedMemory<<<grids, blocks, sizeof(lf) * blocks>>>(d_At, d_B, d_C, nRow, nCol);
else if (i == 4)
MatVecMulWarp<<<grids, blocks>>>(d_At, d_B, d_C, nRow, nCol);
cudaDeviceSynchronize();
cudaEventRecord(end, 0);
cudaEventSynchronize(beg);
cudaEventSynchronize(end);
lf elapsed_time;
cudaEventElapsedTime(&elapsed_time, beg, end);
printf("%fms\n", elapsed_time);
}
cudaFree(d_A);
cudaFree(d_At);
cudaFree(d_B);
cudaFree(d_C);
}
452
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
