温州市网站建设_网站建设公司_表单提交_seo优化

NVCC版本：12

关于GPU你需要了解的那些事

从硬件视角看GPU

之前在CUDA并行规约那篇文章中提到过，在进行CUDA开发时，我们是以Grid，Block和Thread三级的层次结构来组织线程的，那么这三者是如何对应到具体的硬件实现的呢？

从硬件视角下看，一张显卡里有一个GPU，GPU内部有多个流式多处理器（Streaming Multiprocessor，以下简称SM），如下图所示：

image

接下来把视角转向SM内部，每个SM有多个处理器，线程就是在这些处理器上具体执行的。

除此之外，每个SM还有一块共享内存区域，之前文章里提到的共享内存（Shared Memory，以下简称SMEM）就是在这个区域，这个区域只能是SM内部的处理器访问。SM内部的每个处理器又有着只能是自己访问的寄存器（REGS）。

从这里也可以总结出GPU上内存访问的速度排序，寄存器（REGS）是最快的，但是只能线程自身访问；其次是SMEM，但是只能是Block内的线程访问；最慢的是GMEM，GMEM就是我们通过cudaMalloc申请到的内存，所有线程均可访问。

那么上述的线程层级架构又是怎么和这个硬件架构相对应的呢？ 而且Warp在其中又是怎么体现的呢？这就得从线程调度的角度来看一个内核被启动的过程了。

在启动内核时，我们会指定GridDim和BlockDim，这就使得内核有了一定数量的线程块（Block）需要执行，每个Block里有若干个Thread。在进行调度时，GPU会以Block为单位，把一个Block分给一个SM，这时候一个SM可能会被分到多个Block。

接下来就是SM的工作了，一个Block里连续的32个线程为一个Warp，SM会以Warp为单位进行调度，即SM会选择32个连续的线程，然后放到32个处理器上运行。上述过程如下图所示。

image

全局内存访问合并

同一个Warp里的线程有很多有意思的特性，对这些特性加以利用就能够达成不错的优化效果，全局内存访问合并（Global Memory Coalescing）就是其中之一。

这个特性是：如果一个Warp里的线程访问的内存恰好是连续的32个4B的浮点数，那么GPU就只会做一次长度为128B的访存操作， 把128B的数据读取之后分发给32个线程。这里参考资料作者的精美的手绘图可以很形象地说明这一点：

image

一些约定

在开始正式实现前，首先需要把一些容易混淆的设定给明确了。

本文默认所有矩阵都是行优先存储的

本文中的x都是指行下标，y都是指列下标

如下图所示：

image

（注：本文所有的内核实现都只是在大小为4096的方阵下进行了正确性验证，如果要适配任意形状需要考虑很多corner case，这有点偏离主线了，所以本文就暂时不做这方面的适配工作了。）

V42：cuBLAS

这里先放出cuBLAS实现的性能数据，供后续比较和参考

image

V1：Naive Kernel

对于矩阵乘法

，一个最朴素的想法就是让每个线程都计算C中的一个元素，所以只需要一个Block，使用Thread的x和y表示要计算的C的元素坐标，然后2个for循环计算即可。

这个想法没问题，只是在实现的时候，由于一个Block里面最多有1024个线程，所以需要进行一次分块。

具体而言，可以把C分成若干个

的块（Tile），每个Tile交给一个Block进行计算，如下图所示：

image

对于每个线程而言，首先需要根据计算出当前线程需要处理的C的坐标，然后用2个for循环计算结果并写回即可。

至于布局，每个Block里自然是

个线程，而Grid则是需要用向上取整的除法来分配尽量多的Block。

（注：理清楚每个线程需要计算哪些C的元素是不被后面更复杂的分块绕晕的关键）

最终得到的源代码如下所示：

__global__ void MatmulKernelV1(const scalar_t *a, const scalar_t *b, scalar_t *out, uint32_t M, uint32_t N, uint32_t P) {

uint32_t x = blockIdx.x * blockDim.x + threadIdx.x;

uint32_t y = blockIdx.y * blockDim.y + threadIdx.y;

if (x < M && y < P) {

scalar_t tmp = 0;

for (uint32_t k = 0; k < N; k++) {

// out[i][j] = a[i][k] * b[k][j]

tmp += a[x * N + k] * b[k * P + y];

}

out[x * P + y] = tmp;

}

}

void MatmulCoreV1(const scalar_t *a, const scalar_t *b, scalar_t *out, uint32_t M, uint32_t N, uint32_t P) {

dim3 grid(std::ceil(M / 32.0), std::ceil(P / 32.0), 1);

dim3 block(32, 32, 1);

MatmulKernelV1<<<grid, block>>>(a, b, out, M, N, P);

}

实验数据

最终的性能数据如下所示：

image

性能只有cuBLAS的0.39%，可以说是相当拉垮了。

理论分析

这里先插入一段理论分析，来分析一下Naive Kernel可能的性能瓶颈在哪。

首先计算一下理论最快的运行时间，进行一个4096方阵的乘法所需要浮点运算次数为

（因为C有

个元素，每个元素需要进行4096次乘法和加法），大约为137GFLO；

而内存读取最低需要

，约134MB，写入需要

，约67MB。

实验用的显卡浮点数计算性能为870GFLOPS，显存带宽为29GB/s，所以理论上计算最快需要157ms，访存共需要6.9ms，也就是说，理论上来讲，矩阵乘法应该是计算瓶颈的。

但是我们的Naive Kernel似乎并不是这样的，下面来详细分析一下。

在计算次数上，如果不考虑计算下标的开销，那它的计算次数就是和理论最低值相等的；

在内存访问上，实际上每个线程都会访问

次全局内存（GMEM），如果这些访问没有经过任何优化，那么这个内核一共就会有

的访存，需要耗时18.9s，已经是计算的120倍了，所以很显然，目前的首要任务是优化内存访问。

（注：这里内存访问事实上并没有计算的那么多，因为有一些Warp层的自动优化，这个后面马上会提到）

V2：Global Memory Coalescing

这里可以像防止Bank Conflict那样，通过调整每个线程负责的区域来实现Coalescing。

我们首先分析一下V1里面每个线程都在计算C的哪一个元素。通过代码可以知道，线程计算的C的x坐标就是threadIdx.x，y坐标就是threadIdx.y，并且threadIdx是x先变化的，所以第一个线程计算的是(0, 0)，第二个是(0, 1)，以此类推，如下图所示（用背景色来区分T0和T1加载的数据）：

image

可以发现，一个Warp里计算的其实是C中的某一列，那么同一时刻，Warp里的线程访问的A一定不是连续的，所以访问A的部分一个Warp需要

次访存，而访问B的部分，由于一个Warp里的线程在同一时刻访问的都是同一个B，所以这里只会有一次访存开销，那么访问B总共就会有4096次访存。

这里如果我们让一个Warp计算C中的一行会怎么样呢？那么访问A就只需要4096次访存，但是访问B的时候，在同一时刻，线程们访问的数据是连续的，此时就可以触发Global Memory Coalescing，把32次访存压缩为1次，如下图所示：

image

实验数据

image

这里仅仅是对换一下x和y，性能就提升了接近8倍。

但是和cuBLAS相比，还是有不小的差距，目前也还只有cuBLAS性能的2%。

关于实现方式

具体实现时，只需要把x和y对换一下就行了。

参考资料作者在这里的实现是取消了threadIdx.y这个维度，然后把x维度的大小扩展为了1024，之后在线程内部根据threadIdx.x以及BlockSize来计算当前线程对应到C的坐标，如下所示：

constexpr uint32_t BLOCK_SIZE = 32;

__global__ void MatmulKernelV2(const scalar_t *a, const scalar_t *b, scalar_t *out, uint32_t M, uint32_t N, uint32_t P) {

// 相当于首先把OUT分成若干个32*32的Block，V1和V2都是如此，它们的区别在于Block内部的分配方式

// 这里blockIdx.x * BLOCK_SIZE, blockIdx.y * BLOCK_SIZE 就是在定位Block的起始x和y

const uint32_t x = blockIdx.x * BLOCK_SIZE + threadIdx.x / BLOCK_SIZE;

const uint32_t y = blockIdx.y * BLOCK_SIZE + threadIdx.x % BLOCK_SIZE;

......

个人认为这种实现肯定是不如直接对换x和y的，因为这种实现引入了除法和求余数这种非常昂贵的操作，但是实际测试下来两种实现的性能是几乎一致的。

于是我看了下两种实现对应的PTX汇编，如下图所示

image

image

发现两者指令数量几乎都差不多，可能是因为这里BlockSize为2的整数次幂的关系吧，这种特性使得编译器可以对求余数指令做优化，一般的求余数指令是需要用除法+减法来实现的。

如下图所示，可以看到，在把BlockSize换成34之后，确实多出来了一条sub指令。

image

在把BlockSize改成31然后实际运行后，确实出现了性能的损失。

image