CUDA：矩阵乘法与瓦片化处理

并行计算的另一个入门例子：矩阵乘法，附带粗略地说一说瓦片化处理的原因和在这个例子里的应用。

定义矩阵乘法：

其中

大矩阵乘法在 CPU 上的运算是很耗时的，其特点决定了它更适合并行处理。那么如何在 CUDA 上完成运算呢？

确实应该放图上来，但是实在是太懒不想做图……

代码

设两矩阵大小分别是 $m\times n$ 、$n\times k$ , blocksize 是一个 block 的大小，这个值是可以自己设定的，而且有不少学问，这里按下不表，暂且假设是一个小于 $m$ 和 $k$ 的数。这里略去分配显存和拷贝数据的过程，直接看 De­vice Code：

__global__
void MatMUL(int m,int n,int k,float *A,float *B,float *C) {
    int Row = blockIdx.y*blockDim.y + threadIdx.y;
    int Col = blockIdx.x*blockDim.x + threadIdx.x;

    if ((Row < m) && (Col < k)) {
        float c = 0.0;
        for (int i = 0; i < n; i++) {
            c += A[Row*n + i] * B[i*k + Col];
        }
        C[k*Row + Col] = c;
    }
}

__host__
int main(){
    ……
    dim3 DimGrid(k / blocksize + 1, m / blocksize + 1, 1);
    dim3 DimBlock(blocksize, blocksize, 1);

    //非瓦片化
    MatMUL <<<DimGrid, DimBlock >>> (m, n, k, d_A, d_B, d_C);
    ……
}

简析

从 DimGrid 与 DimBlock 可知，这里的线程组织方式是二维的，定位某个线程的方法是：

int Row = blockIdx.y*blockDim.y + threadIdx.y;
int Col = blockIdx.x*blockDim.x + threadIdx.x;

对应 $C$ 矩阵中的第 $Row$ 行 $Col$ 列的元素，其值即通过一个循环累加得到。注意虽然 A 与 B 矩阵是二维矩阵，但是它们的存储方式仍然是从上到下从左到右线性存储的，所以可以通过 A[Row*n + i] 这样的方式来访问矩阵元素值。

瓦片化处理

上面的代码中，GPU 执行过程中读写均发生在 global mem­ory 上，但 global mem­ory 速度慢于 block 内的 shared mem­ory。瓦片化的思想是将某些频繁使用的值放进 shared mem­ory 中，减少对 global mem­ory 的访问次数，以此提高执行速度。但由于 shared mem­ory 的大小有限，无法每次都把所有的数据载上去，因而每次只将 “一片” 元素值载入，故称「瓦片化处理」。

代码

__global__
void MatMUL_Tile(int m, int n, int k, float *A, float *B, float *C) {
#define T 7
    __shared__ float ds_A[T][T];
    __shared__ float ds_B[T][T];

    int Row = blockIdx.y*T + threadIdx.y;
    int Col = blockIdx.x*T + threadIdx.x;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    float c = 0.0;

    if ((Row < m) && (Col < k))
    {
        for (int t = 0; t < n / T + 1; t++)
        {
            if ((Row < m) && (t*T + tx < n))    ds_A[ty][tx] = A[Row*n + t*T + tx];
            else                                ds_A[ty][tx] = 0;
            if ((t*T + ty < n) && (Col < k))    ds_B[ty][tx] = B[(t*T + ty)*k + Col];
            else                                ds_B[ty][tx] = 0;
            __syncthreads();
            for (int i = 0; i < T; i++)
            {
                c += ds_A[ty][i] * ds_B[i][tx];
            }
            __syncthreads();
        }
        C[Row*k + Col] = c;
    }
}

T 代表每个瓦片的尺寸，它与 block 的大小一致。__syncthreads(); 的作用是让 block 内的线程保持一样的执行进度，在遍历过程中瓦片中的值会被后面的操作覆盖，同步所有的线程可以防止读取、写入混乱。

要特别注意边界条件，保证线程覆盖所有的待计算元素。当瓦片大小与矩阵大小不成整数比时，为覆盖所有元素，瓦片势必有一部分在矩阵之外，这部分瓦片的元素应该赋值为 0，否则可能发生意想不到的后果。