ECE408 Lecture 6 Generalized Tiling and DRAM Bandwidth

这次回顾ECE408 Lecture 6，这次介绍了Tiling的一般形式以及DRAM带宽。

课程主页：

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

搬运视频：

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

处理任意大小的矩阵

第5讲中的分块矩阵乘法内核只能处理维度是分块维度的倍数的矩阵；
但是，实际应用程序需要处理任意大小的矩阵；
我们可以将行和列填充（添加元素）成块大小的倍数，但是会有很大的空间和数据传输时间开销；
我们可以在代码中添加显式检查来处理边界；

例子

考虑TILE_WIDTH = 2, width = 3的矩阵乘法。

考虑载入Block(0, 0)的第2个tile，边界情况需要特殊处理：

计算：

考虑载入Block(1, 1)的第1个tile，边界情况需要特殊处理：

计算：

2x2例子中的主要情形

计算有效$P$元素但使用无效输入的线程
- Block(0, 0)的第2个tile；
计算无效$P$元素的线程
- Block(1, 1), Thread(1, 0), 不存在的行；
- Block(1, 1), Thread(0, 1), 不存在的列；
- Block(1, 1), Thread(1, 1), 不存在的行/列；

一个“简单”的解决方案

无效的输入元素，“加载”0
- 基本原理：0 值将确保乘加步骤不会影响输出元素的最终值；
无效的输出元素，不要更新全局内存
- 仍然可以执行pvalue计算（部分点积），但不会写入内核末尾的全局内存；

代码

// 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 < (ceil((float)Width/TILE_WIDTH)); ++q) {
    // Collaborative loading of M and N tiles into shared memory
    // 保证输入有效, 否则为0
    if (Row < Width && (q*TILE_WIDTH+tx) < Width)
    	subTileM[ty][tx] = M[Row*Width + q*TILE_WIDTH+tx];
    else
    	subTileM[ty][tx] = 0;
    
    // 保证输入有效, 否则为0
    if (Col < Width && (q*TILE_WIDTH+ty) < Width)
    	subTileN[ty][tx] = N[(q*TILE_WIDTH+ty)*Width+Col];
    else
    	subTileN[ty][tx] = 0;
    __syncthreads();
    
    for (int k = 0; k < TILE_WIDTH; ++k)
        Pvalue += subTileM[ty][k] * subTileN[k][tx];
        __syncthreads();
    }
}
// 保证输出有效, 否则不更新
if (Row < Width && Col < Width)
    P[Row*Width+Col] = Pvalue;

一些要点

对于每个线程，如下情形略有不同
- 加载$M$元素；
- 加载$N$元素；
- 计算/存储输出元素；
对于大矩阵，控制divergence的影响应该很小；
非方阵矩阵怎么样？

局部性

这里老师主要介绍了局部性的概念，即加载二维矩阵中的某个元素，会一次性加载该元素同一行的周围几个元素，如果后续会访问这些已被加载的元素，则代码的速度很快，原因是二维矩阵是以一维数组的形式存储在内存中，所以同一行的元素之间距离较近，而同一列的元素之间距离较远：

代码回顾

现在从局部性的角度重新分析之前的代码：

__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;
    }
}