ECE408 Lecture 5 Locality and Tiled Matrix Multiply

这次回顾ECE408 Lecture 5，这次介绍了Tiled矩阵乘法。

课程主页：

https://wiki.illinois.edu/wiki/display/ECE408

搬运视频：

https://www.youtube.com/playlist?list=PL6RdenZrxrw-UKfRL5smPfFFpeqwN3Dsz

本讲目标

学习评估全局内存访问的性能影响；
为MP3做准备
- tiled矩阵乘法；
了解tiling的好处；

在内存带宽为150 GB/s的设备上性能如何？

所有线程访问其输入矩阵元素所在的全局内存
- 每个浮点乘加 (2 fp ops) 两次内存访问（8字节）；
- 每个FLOP使用4B内存；
- 150 GB/s将代码限制在37.5 GFLOPS；
  - 150 / 4 = 37.5；
实际代码运行速度约为25 GFLOPS；
需要大幅减少内存访问以接近超过1,000 GFLOPS的峰值；

即主要性能瓶颈是全局内存的读取：

通过重用避免带宽（BW）瓶颈

$M$和$N$的每个元素在计算$P$时使用$\mathrm{Width}$次；
为避免BW限制，通过利用每个SM共享内存来利用数据重用；
将数据划分为适合共享内存的块；
每个线程块使用以下策略：
- 将图块从全局内存加载到共享内存；
- 从共享内存中对tile执行计算，每个线程都可以有效地访问任何数据元素；
- 将结果从共享内存复制到全局内存；

共享内存tiling基本思想

SM内存架构

架构示意图：

性能说明：

寄存器（~1 个周期）；
共享内存（~5 个周期）；
高速缓存/常量内存（~5 个周期）；
全局内存（~500 个周期）；

实现

整体思路

找到数据重用率较高的全局数据块；
将tile从全局内存加载到共享内存；
块中的线程从共享内存访问它们的数据；
移动到下一个block/tile；

共享内存

添加__shared__即为共享内存：

__global__
void MatrixMulKernel(float* d_M, float* d_N, float* d_P, int Width)
{
    // 在block中的线程中共享
    __shared__ float subTileM[TILE_WIDTH][TILE_WIDTH];
    __shared__ float subTileN[TILE_WIDTH][TILE_WIDTH];

Tiled乘法

将内核的执行分成多个阶段，以便每个阶段的数据访问都集中在$M$和$N$的一个tile上；
对于每个tile：
- 阶段1：将$M$和$N$的tile加载到共享内存中；
- 阶段2：计算tile $P$的部分点积；

图示：

上图的含义为，先处理蓝色的tile，然后处理橘色的tile，以此类推，我们从下面的例子中加深理解。

例子

以Block(0, 0)为例，这里假设$\text{BLOCK_WIDTH}=2$。

步骤1：加载数据至共享内存

一开始加载的是$M$的Block(0, 0)，$N$的Block(0, 0)：

步骤2：计算共享内存的矩阵乘法

这一步计算矩阵乘法即可：

步骤3：加载其他数据至共享内存

完成后我们加载$M$的Block(0, 1)，$N$的Block(1, 0)：

步骤4：计算共享内存的矩阵乘法

根据之前的例子，给出具体实现。

具体实现

内容回顾

在给出具体代码之前，首先回顾上一讲的代码：

__global__
void MatrixMulKernel(float *d_M, float *d_N, float *d_P, int Width)
{
    // Calculate the column index of d_P and d_N
    int Col = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Calculate the row index of d_P and d_M
    int Row = blockIdx.y * blockDim.y + threadIdx.y;
    
    if ((Row < Width) && (Col < Width)) {
        float Pvalue = 0;
        // each thread computes one element of d_P
        for (int k = 0; k < Width; ++k)
        	Pvalue += d_M[Row*Width+k] * d_N[k*Width+Col];
        d_P[Row*Width+Col] = Pvalue;
    }
}

我们使用共享内存和tile，优化的是for (int k = 0; k < Width; ++k)这一层循环，因为元素会有重复访问的问题；具体做法是将d_M[Row*Width+k], d_N[k*Width+Col]分块加载到共享内存（共享内存大小为BLOCK_WIDTH x BLOCK_WIDTH），因此循环变为for (int q = 0; q < Width/TILE_WIDTH; ++q)，q为分块后的坐标，而块内部的坐标，用tx, ty表示（对应于threadIdx.x, threadIdx.y），因此，二维坐标为：

M_shared_memory: (Row, q * TILE_WIDTH + tx)
N_shared_memory: (q * TILE_WIDTH + ty, Col)

最后转换为一维坐标即可。

现在将之前的讨论总结如下：

阶段1：加载tile

块中的所有线程将$M$和$N$的某个元素和加载到共享内存中；
将加载的元素分配给每个线程，以便合并每个warp中的访问；

代码：

tx = threadIdx.x;
ty = threadIdx.y;
subTileM[ty][tx] = M[Row*Width + (q*TILE_WIDTH + tx)]
subTileN[ty][tx] = N[(q*TILE_WIDTH+ty) * Width + Col]

图示：

阶段2：计算部分矩阵乘法

计算第$q$步的矩阵乘法：

PValue += subTileM[ty][k] * subTileN[k][tx];

图示：

不正确的版本

根据之前讨论，这里先给出一个正确版本的代码（后续会给出不正确的原因）：

__global__ void MatrixMulKernel(float* M, float* N, float* P, int Width)
{
    __shared__ float subTileM[TILE_WIDTH][TILE_WIDTH];
    __shared__ float subTileN[TILE_WIDTH][TILE_WIDTH];
    int bx = blockIdx.x; int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;
    
    // Identify the row and column of the P element to work on
    int Row = by * TILE_WIDTH + ty;
    int Col = bx * TILE_WIDTH + tx;
    float Pvalue = 0;
    
    // Loop over the M and N tiles required to compute the P element
    // The code assumes that the Width is a multiple of TILE_WIDTH!
    for (int q = 0; q < Width/TILE_WIDTH; ++q) {
        // Collaborative loading of M and N tiles into shared memory
        subTileM[ty][tx] = M[Row*Width + (q*TILE_WIDTH+tx)];
        subTileN[ty][tx] = N[(q*TILE_WIDTH+ty)*Width+Col];

        for (int k = 0; k < TILE_WIDTH; ++k) {
        	Pvalue += subTileM[ty][k] * subTileN[k][tx];
        }
	}
    P[Row*Width+Col] = Pvalue;
}

这个代码之所以有问题，是因为，如果跑得快的进行先进入$k$的循环，那么共享内存可能没有加载完，因此结果不正确，所以我们需要一个方法让各个进程之间保持一致，即同步。

基于Barriers的批量同步步骤

有时我们需要所有线程在任何线程继续之前赶上我们代码中的某个点；
barrier是一个同步点：
- 每个线程调用一个函数进入barrier；
- 线程在屏障函数中阻塞（睡眠），直到所有线程都调用它；
- 在最后一个线程调用barrier函数后，所有线程继续越过barrier；

图示：

Barrier同步

Barrier同步是CUDA中的API函数调用__syncthreads()；
同一块中的所有线程都必须到达__syncthreads()才能继续前进；
可用于协调tiled算法
- 确保加载了tile的所有元素；
- 确保对元素的某些计算是完整的；

添加了同步后的代码

__global__ void MatrixMulKernel(float* M, float* N, float* P, int Width)
{
    __shared__ float subTileM[TILE_WIDTH][TILE_WIDTH];
    __shared__ float subTileN[TILE_WIDTH][TILE_WIDTH];
    int bx = blockIdx.x; int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;
    
    // Identify the row and column of the P element to work on
    int Row = by * TILE_WIDTH + ty;
    int Col = bx * TILE_WIDTH + tx;
    float Pvalue = 0;
    
    // Loop over the M and N tiles required to compute the P element
    // The code assumes that the Width is a multiple of TILE_WIDTH!
    for (int q = 0; q < Width/TILE_WIDTH; ++q) {
        // Collaborative loading of M and N tiles into shared memory
        subTileM[ty][tx] = M[Row*Width + (q*TILE_WIDTH+tx)];
        subTileN[ty][tx] = N[(q*TILE_WIDTH+ty)*Width+Col];
		__syncthreads();
        for (int k = 0; k < TILE_WIDTH; ++k) {
        	Pvalue += subTileM[ty][k] * subTileN[k][tx];
        }
        __syncthreads();
	}
    P[Row*Width+Col] = Pvalue;
}

tiles的使用使得性能瓶颈转移

回顾我们的示例 GPU：1,000 GFLOP/s，150 GB/s 内存带宽
16x16 tiles重复使用每个操作数16次
- 将全局内存访问减少16倍；
- 150GB/s 带宽支持(150/4)x16 = 600 GFLOPS；
32x32块重复使用每个操作数32次
- 将全局内存访问减少32倍；
- 150GB/s带宽支持(150/4)x32 = 1,200 GFLOPS！
- 内存带宽不再是瓶颈！

SM限制也是一个因素

共享内存大小
- 取决于实现；
- Maxwell中每个SM 64kB（每个块最大48kB）；
给定TILE_WIDTH为 16（256 个线程/块）
- 每个线程块使用2x256x4B = 2kB的共享内存；
- 活动的block限制为64/2=32个；
- 每个SM最多有2048个线程；
- 将块限制为2048/256=8个；

另一个不错的选择：32x32 tiles

给定TILE_WIDTH为32（1,024 个线程/块），
每个线程块使用2x1024x4B = 8kB 共享内存，
将活动块限制为64/8=8个；
每个SM最多有2048个线程；
将块限制为2048/1024=2个；

如何查询当前显卡的信息？使用设备查询

系统中的设备数量

```cpp
int dev_count;
cudaGetDeviceCount(&dev_count);


- 设备能力

  - ```cpp
    cudaDeviceProp dev_prop;
    for (i = 0; i < dev_count; i++) {
    	cudaGetDeviceProperties(&dev_prop, i);
    	// decide if device has sufficient resources and capabilities
    }

cudaDeviceProp是内置的C结构类型
- dev_prop.dev_prop.maxThreadsPerBlock
- dev_prop.sharedMemoryPerBlock