矩阵乘法的优化

4.7.1.1 matrixMul 0

下面给出矩阵乘的基本实现。如图 4-13 所示，基本的矩阵乘法使用了带状划分，每个线 程负责读取 A 中的一行和 B 中的一列，并计算出 C 中相应位置上的值。于是，整个 kernel 从

全局存储器对 A 矩阵进行了 B.width 次读取、对 B 矩阵进行了 A.height 次读取，才能完成对 矩阵 C 的计算。

为讲解方便起见，假设数组在每个维度上的尺寸都是 BLOCK_SIZE 的整数倍。本例中，

BLOCK_SIZE 为 16，block 尺寸为（16,16），grid 尺寸为（8,5）。数组 A、B、C 的大小如图 4-13 所示，图中的每个小矩形块（也称 tile 或 blocking）是 16*16 的数组元素块。

图 4-13 C = A * B

以下是主机端代码，A、B 矩阵均按行存储，且各维度均是 BLOCK_SIZE 的整数倍。

#define BLOCK_SIZE 16

void MatMul(const Matrix A, const Matrix N, Matrix C) {

//为 A、B 加载做好准备，分配显存空间并初始化 Matrix d_A;

d_A.width = A.width; d_A.height = A.height;

size_t size = A.width * A.height * sizeof(float);

cudaMalloc((void**)&d_A.elements, size);

cudaMemcpy(d_A.elements, A.elements, size,cudaMemcpyHostToDevice);

略,d_B 的分配

//为结果矩阵 C 分配空间 Matrix d_C;

d_C.width = C.width; d_C.height = C.height;

size = C.width * C.height * sizeof(float);

cudaMalloc((void**)&d_C.elements, size);

//启用 kernel

dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);

dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);

MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);

// 将结果矩阵 C 从显存中拷贝回内存

cudaMemcpy(C.elements, Cd.elements, size, cudaMemcpyDeviceToHost);

//释放显存空间 cudaFree(d_A.elements);

cudaFree(d_B.elements);

cudaFree(d_C.elements);

}

以下是 matrixMul 0 版本的 kernel 实现，每个线程都负责 C 中一个位置的元素值的计算，

例如 tid 线程负责 C 中(r, c)位置的值，for()循环完成 A 中第 r 行与 B 中第 r 列对应元素的乘加 操作。整个计算过程如图 4-14 所示。

//kernel 定义

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C) {

float Cvalue = 0;

int row = blockIdx.y * blockDim.y + threadIdx.y;

int col = blockIdx.x * blockDim.x + threadIdx.x;

// 每个线程负责计算 C 中的一个元素，并将值累加到 Cvalue for (int e = 0; e < A.width; ++e)

Cvalue += A.elements[row * A.width + e]

* B.elements[e * B.width + col];

C.elements[row * C.width + col] = Cvalue;

}

图 4-14 矩阵相乘中 C 中元素的计算

4.7.1.2 matrixMul 1

访存共享存储器的延迟远小于全局存储器，并且使用共享存储器进行线程间通信，还可 以通过更灵活的算法获得更好的性能。以下代码对矩阵进行了棋盘划分，使用了共享存储器来

实现矩阵乘法。与上一个版本相比，这一个版本还能避免按列读取 B 中元素，引起的非合并 访问造成的性能下降。

//---matrixMul_kernel.cu--- __global__ void

matrixMul( float* C, float* A, float* B, int wA, int wB) {

在这个实现中，每个 block 负责数组 C 中的一个方块 Csub 的计算，其中每个线程负责计 算 Csub 的一个元素。如图 4-15 所示，Csub 是两个矩形矩阵的乘：A 的 sub-matrix，维度是(A.width, block_size)，与 Csub 有相同的行索引；B 的 sub-matrix，维度是(block_size, A.width)，与 Csub 有相同的列索引。为了适应设备资源，两个矩形矩阵被分为大小为 block_size 的正方形矩阵，

Csub 的计算就是这些小正方形矩阵的乘积的和。这些乘积如下计算：首先从 global memory 加载两个相应的正方形矩阵到 shared memory，一个线程负责加载一个元素；接下来每个线程 负责计算乘积中的一个元素。每个线程累积每个这样乘积的结果到一个寄存器中，一旦结束就 将结果写回 global memory。

图 4-15 分块计算矩阵乘法

这种棋盘划分方式利用了高速的 shared memory，同时也节约了大量 global memory 带宽，

因为 A 只被读 B.width/block_size 次，B 也只被读 A.height/block_size 次。

具体是这样实现的：例如 block(0,0)需要由 A 的一行 tile 和 B 的一列 tile 计算得到，如图 4-15 中矩形框标注的部分。对于每一个 block：

 第一次循环，从 global memory 中取出 A 的 u 块和 B 的 u 块，分别放入 shared memory 中的 AS[16][16]、BS[16][16]中，然后 block 中的每个 thread 根据 AS、BS、tx、ty 计 算出一个 Csub。

 第二次 for 循环中，从 global memory 中取出 v 和 v 块并放入 AS、BS 中，每个 thread 在原来的 Csub 基础上再进行加和。

 依此类推。

 最后一次循环中，w 和 w 被放入 AS、BS，每个 thread 做最后的 Csub 加和，这时的 Csub 就是每个 thread 计算出来的最终结果了，在这里相当于计算出了（0，0）tile。

最后，每个 block 中的 16*16 个线程，对应完成 C 中 16*16 个位置上（即 1 个 tile）

的输出，输出到显存中。

这个代码实现中不存在 bank conflict 问题。因为在进行两个子块乘加的 for()循环中，每 half-warp 访问的是 AS[]一行的元素，分别分布在 16 个 bank，对 BS[]的访问也分布在 16 个不 同的 bank，不存在 bank conflict 问题。

但这个版本仍有许多地方需要改进，例如数组大小必须是 BLOCK_SIZE 的整数倍，每个 tile 会被重复取出多次。访问显存的高延迟，使得我们必须思考如何将“好不容易”取到的数 据得到充分使用。限于篇幅，这里仅给出一些提示：

（1）可以使用 cudaMallocPitch()来解决数组各维大小不是 2 的幂次方的问题。它会在 global memory 中做好数组宽度、起始地址的工作，自动以最佳倍数来分配内存。

（2）可以使用 cudaMemcpy2D()进行二维数组的复制，显存上数组 pitch 与内存上数组 pitch 可以不同。

（3）如果是浮点运算，可以借助 Kahan’s Summation Formula 提高计算精度。

4.7.1.3 matrixMul 2

除了共享存储器，我们也应该充分利用 GPU 的寄存器资源。register 是矢量计算单元（一 条指令可取 4 个），shared memory 是标量单元（一条指令只可取一个），因此理论上使用寄存 器更加快速。事实上，每个 SM 都拥有 8K 或 16K 的寄存器资源，线程少的话每个线程分得的 register 还是很可观的。根据这个思想，下面给出了改进的 multrixMul 2 版本矩阵乘法。

在这个版本的代码中，A、B、C 均按列存储，并且数组各维大小均是 BLOCK_SIZE 整数 倍，A、B 矩阵均为 64*64，而 block 的维度是(16,4)。B 被划分为 16*16 的子块，并借助共享 存储器进行计算；对 A、C 按列做 vector，这样保证对显存中 C 的 stride-1 访问。

__global__ void sgemmNN( const float *A, int lda, const float *B, int ldb, float* C, int ldc, int k) {

//计算每个线程的起始指针

A += blockIdx.x * 64 + threadIdx.x + threadIdx.y*16;

B += threadIdx.x + ( blockIdx.y * 16 + threadIdx.y ) * ldb;

C += blockIdx.x * 64 + threadIdx.x + (threadIdx.y + blockIdx.y * ldc ) * 16;

__shared__ float bs[16][17];

float c[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};

const float *Blast = B + k;

bs[threadIdx.x][threadIdx.y+i] = B[i*ldb];

B += 16;

c[12] += A[0]*bs[i][12];c[13] += A[0]*bs[i][13];

c[14] += A[0]*bs[i][14];c[15] += A[0]*bs[i][15];

} 成了 Asub(16*64)与 bs[][]的乘加，并将结果累加于 Csub(16*64)中。而 while()内的一轮循环就 完成了 Asub(16*64 子矩阵)与 B 中一个 16*16 子矩阵的乘加。请注意，此时 64 个线程的 c[]

组成的 16*64 的 Csub 并不是最终结果，还需要以 while(B<Blast)来遍历整个 Bsub(16*64)。在 本例中，在矩阵 A、B 都是 64*64 的情况下，while()循环会执行 4 次。接下来的执行过程如图 4-16（d）所示，对 A 中每一个 Asub(16*64)与 B 中相应的(16*16)矩阵块进行计算。循环 4 次 以后最后累加得到的即为最终的 Csub，此时将结果传回全局存储器即可。

这与早期矢量处理器上的两个算法——Agarwal、Gustavson 为 IBM 3090 Vector Facility，

以及 Anderson et al.为 Cray X1 设计的算法思想相同。A、C vector 存储于矢量 register 或 cache 中，B blocking 存储于另外一个可被 A 中不同元素共享的内存上，这说明了 CUDA 与一些传 统超级计算机在架构上存在相似性。

（a）

（b） （c）

（d）

图 4-16 利用寄存器进行矩阵乘法