第6章 矩阵乘法的示例
6.1 概述
计算两个维度分别为(wA, hA)和 (wB, wA) 的矩阵A和B的乘积C的任务以下列方式分为多个线程：
每个线程负责计算C 的一个平方子矩阵Csub ；
块内的每个线程负责计算Csub的一个元素。
选择Csub的维度block_size等于16，以便每块的线程数是warp大小的倍数（参见5.2），并且保持低于每块的最大线程数（参见附录A）。
如图6-1所示，Csub 等于两个矩形矩阵的乘积：维度为(wA, block_size)的子矩阵A，与Csub具有相同的行索引，维度为(block_size, wA)的B的子矩阵，与Csub具有相同的列索引。为了适应设备的资源，这两个矩形矩阵可根据需要划分为许多维度为block_size的平方矩阵，并且Csub计算为这些平方矩阵的乘积之和。其中每个乘积的执行过程是：首先使用一个加载每个矩阵的一个元素的线程，将两个相应的平方矩阵从全局内存加载到共享内存，然后让每个线程计算乘积的一个元素。每一线程将其中每个乘积的结果累计到寄存器中，执行完毕后，将结果写入全局内存。
通过以这种方式分块计算，我们可以有效利用快速的共享内存，并节省许多全局内存带宽，因为A和B仅从全局内存读取(wA / block_size)次。
尽管如此，编写此示例是为了清楚地说明各种CUDA编程原则，并非是为了为一般的矩阵乘法提供高性能的内核，所以不应如此构造。
每一线程块计算C的一个子矩阵Csub。块内的每一线程计算Csub的一个元素。
图6-1. 矩阵乘法
6.2 源码清单
6.3 源码攻略
源码包含下列两个函数：
Mul()，作为Muld()的包装器的主机函数。
Muld()，在设备上执行矩阵乘法的内核。
6.3.1 Mul()
Mul()接受下列输入：
指向A和B的元素的主机内存的两个指针，
A的高度和宽度，B的宽度，
指向应该写入C的主机内存的指针。
Mul()执行下列操作：
使用cudaMalloc()将足够的全局内存分配到库A、B和C中；
使用cudaMemcpy()将A和B从主机内存复制到全局内存；
调用Muld()在设备上计算C；
使用cudaMemcpy()将C从全局内存复制到主机内存；
使用cudaFree()释放为A、B和C分配的全局内存。
6.3.2 Muld()
除了指针指向设备内存而非主机内存之外，Muld()与Mul()具有相同的输入。
对于每个块，Muld()迭代处理所有需要计算Csub的A和B的子矩阵。在每次迭代中，此函数：
将A的一个子集和B的一个子集从全局内存加载到共享内存中；
同步以确保两个子矩阵都由块内的所有线程完全加载；
计算两个子集的乘积并将其加到上一次迭代期间获得的乘积中；
再次同步以确保在开始下一次迭代之前两个子集的乘积已经完成。
按照5.1.2.1和5.1.2.4所述，编写Muld()是为了最大化内存性能。
的确，假设wA和wB是16的倍数（如5.1.2.1所建议的），则确保了全局内存合并，因为a、b和c都是BLOCK_SIZE的倍数，BLOCK_SIZE等于6。
对于每个半warp，也没有任何共享内存库冲突，所有线程的ty和k都是相同的，tx在0到15之间变化，因此对于内存访问As[ty][tx]、Bs[ty][tx]和Bs[k][tx]，每个线程都访问一个不同的库，对于内存访问As[ty][k]，每个线程都访问同一个库。
