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

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

定义矩阵乘法:

$C_{m\times k}=A_{m\times n}\cdot B_{n\times k}$

其中

$C_{i,j}=\sum_{k=1}^nA_{i,k}\times B_{k,j}$

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

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

代码

设两矩阵大小分别是 $m\times n$ 、$n\times k$ , blocksize 是一个 block 的大小,这个值是可以自己设定的,而且有不少学问,这里按下不表,暂且假设是一个小于 $m$ 和 $k$ 的数。这里略去分配显存和拷贝数据的过程,直接看 Device 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);
    ……
}

简析

DimGridDimBlock 可知,这里的线程组织方式是二维的,定位某个线程的方法是:

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 memory 上,但 global memory 速度慢于 block 内的 shared memory。瓦片化的思想是将某些频繁使用的值放进 shared memory 中,减少对 global memory 的访问次数,以此提高执行速度。但由于 shared memory 的大小有限,无法每次都把所有的数据载上去,因而每次只将 “一片” 元素值载入,故称「瓦片化处理」。

代码

__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,否则可能发生意想不到的后果。

Tags:c++cuda代码code
上一篇
打赏
下一篇

添加新评论