weixin_55035144的博客

在上面的代码例子中，假设对于16×16的块来说，每个块使用16×16×4 = 1KB的存储空间来存放Mds，同理需要1KB空间存放Nds，每个块需要使用2KB的共享存储空间，但是因为线程的硬件限制，每个SM最多容纳768个线程，使得最多只有3个块，那么只能使用共享存储器的6KB的空间，就浪费了10KB的空间了，不过硬件的限制是因显卡的更新换代而不同的，比如GT200中每个SM最大支持1024个线程。用Md_(0,1),Md(1,1)和Nd_(0,0)和Nd(0,1)来计算Pd_(0,1)，Pd(1,1)。

在内存存储中，默认矩阵是按照行优先储存的，即矩阵的每一列在内存中是连续的。行优先矩阵储存中行数据是不连续的。而对于列主序的矩阵，是按照列优先储存的，即矩阵的每一行在内存中是连续的。

1.向量外积的矩阵乘法要取A中一列，B中一行。但A是行主序存储的，要把它改成列主序存储，来使得。这也能避免Bank Conflict。A的长度为TM，B的长度为TN。一次外积共TM * TN次乘累加。一共需要循环BK次，取BK次循环，每次循环中进行TM * TN次乘累加。得到的结果累加为C矩阵的结果。

为了高效计算矩阵相乘，并进行并行化放到GPU中运算。需要将矩阵切分成子矩阵，用子矩阵相乘的结果组合为原矩阵相乘的结果：上图是拆分矩阵的方法，M表示矩阵，X方向和Y方向的两个维度都是未知的。P表示横条或竖条，X方向或Y方向有一个方向的维度是极小的。B表示block块，X方向和Y方向的两个维度都是极小的。为了减小单个子矩阵计算量，要拆开A的整行和B的整列。不能让A的整行和B的整列作为子矩阵放入缓存进行计算。因此下图中第二列的Fig8和Fig10拆得最好，把A按列拆，使A的行不再完整，把B按行拆，

下图左边是一个非常简单的多层内存模型，只有寄存器/ cache/RAM。在这种简单模型结构下考虑优化GEBP，Cmc,n+=Amc,kcBkc,nCmc,n+=Amc,kcBkc,n，其中3个假设基于以上三点假设，上图中GEBP的RAM和cache之间的数据搬移开销为mckc+kcn+2mcnmemops而Cj:=ABj+Cj的计算量为2mckcnflops，那么计算量和数据搬移的比例。问题变成，现实中kckc的选择还受一些其他因素制约，我们将在6.3章节看到。类似地可以分析GEPB和GEDOT操作。

以下列代码为例：分配的GPU内存。它引用的存储器类型是，因为它并不关心数据类型。一旦我们有我们刚刚分配的内部的存储器的地址，我们需要以某种方式返回给调用者。但所有CUDA API函数都返回错误代码（或成功代码），并且任何其他返回的数据都必须由参数引用传回。既然你想存储一个指针值以存储结果的变量，传递一个指针。使用类型。

可以将 A 和 B 矩阵先搬运到 Shared Memory（SM 中低延迟的 on-chip memory，block 内线程共享，附 NVIDIA GPU 内存结构图）中降低访存的开销，这的确是一个很好的思路，但是这只能将访存代价从几百 cycle 降低到几十 cycle，并不改变问题的本质。相对于之前的情况，计算指令的占比大大提高了。问题的关键在于主体循环由两条 Load 指令与一条 FMA 指令构成，计算指令只占总体的 1/3，计算访存比过低，最终导致了访存延迟不能被隐藏，从而性能不理想。

其中 A 的维度为 Mxw，B 的维度为 wxN，C 的维度为 MxN。为了保持内核简单，M 和 N 是 32 的倍数，因为当前设备的warp size (w) ，即warp内thread数量是 32。因此，对于wxw的瓦片tile而言，A是列矩阵，B是行矩阵，C是它们的。这里启动 N/w x M/w 块的网格，其中每个线程块根据 A 的同一tile和 B 的同一tile计算 C 中。wxw-thread 块中的每个线程计算 C 的tile中的一个元素，将 A 的行乘以 B 的列，然后将其写入 C。

对于矩阵C中的每一个tile，矩阵A与矩阵B中的tile只需要读取一次，这样就可以达到 O(N) 的计算强度。理想情况下性能应当受限于算术吞吐量，确切来说对于一个大型方阵，即 M=N=K ，矩阵乘法的数学操作复杂度为 O(N3) ，而数据量为 O(N2) ，同时计算强度为 N （这里的计算强度可以被理解为同一个数据被重复利用的次数）。GEMM就是指计算 C=A∗B+C ，其中 A 、B和C都是矩阵，A是 M×K 的矩阵，B是 K×N 的矩阵， C是 M×N 的矩阵。这就是向量内积和向量外积的最大区别。

从上图中我们可以看到三种处理方法。第一种是将A和B矩阵分块，第二种方法是将C和B矩阵分块，第三种方法是将C和A矩阵分块。GEMM的子任务是GEPP或GEMP；最小粒度的任务是GEBP或GEPB或点乘。

CUDA缓存包括L1缓存和L2缓存。常规的路径是一级和二级缓存，需要使用常量和只读缓存的需要在代码中显式声明。但是提高性能，主要还是要取决于访问模式。控制全局加载操作是否通过一级缓存可以通过编译选项来控制，当然比较老的设备可能就没有一级缓存。编译器禁用一级缓存的选项是：1编译器启用一级缓存的选项是：1当一级缓存被禁用的时候，对全局内存的加载请求直接进入二级缓存，如果二级缓存缺失，则由DRAM完成请求。

矩阵乘法的优化需要将矩阵切分成子矩阵，用子矩阵相乘的结果组合为原矩阵相乘的结果：上图是拆分矩阵的方法，M表示矩阵，X方向和Y方向的两个维度都是未知的。P表示横条或竖条，X方向或Y方向有一个方向的维度是极小的。B表示block块，X方向和Y方向的两个维度都是极小的。为了减小单个子矩阵计算量，要拆开A的整行和B的整列。不能让和作为子矩阵。因此下图中第二列的Fig8和Fig10拆得最好，把A按列拆，使A的行不再完整，把B按行拆，使B的列不再完整。

我们一次计算C矩阵的一个元素，这个时候需要遍历A矩阵的一行和B矩阵的一列并做乘加运算。一次计算4个元素(我们在寄存器中累加C的元素，并对a的元素使用寄存器),用指针来寻址B中的元素。NEON指令集优化, 并且为了保持较小问题规模所获得的性能，我们分块矩阵C（以及相应的A和B）一次计算C中的4x4小块(我们在寄存器中累加C的元素，并对a的元素使用寄存器)在MMult_4x4_6的基础上用指针来寻址B中的元素。一次计算C中的4x4小块,将16个循环合并一个。一次计算4个元素(将4个循环合并为1个)

因为GPU存在Host和Device内存，所以先申请host内存h_a,h_b，存放a，b的一维矩阵的内容（也可以生成随机数），并申请host内存h_c存放c的计算结果。如果没有block的概念，要同时进行同步、通信、协作时，整体的核心都要产生等待的行为，如要进行扩展时，扩展的越多等待也越多。block有多大，用blockDim表示它有多少个thread，具体分为blockDim.x，blockDim.y，blockDim.z。先分配 源地址空间a,b，目的地址空间c，并生成a，b的随机数。

OpenCL与CUDA存储模型。

而若一个 thread 并不只计算一个结果，而是计算 4x4=16个结果，就要从A和B中分别取出4个数据，共8个数据。访存比变为16/8=2，是上面的4倍。C矩阵为128*128大小的矩阵。C矩阵被分成了四份，每份的尺寸都为4*4，使用同一个线程计算这四份4*4大小区域的FMA计算。计算一次 FMA（乘累加）为一次运算，而各读取 A 和B中一个元素为1+1=2次运算。注意，这里将全局内存中的A矩阵存入共享内存smemA中时进行了矩阵转置。

具体的过程还需要进一步研究。

MMult_4x4_3:一次计算C矩阵的16个元素:

为了避免bank conflict，可以使用padding的方式将数组中的每一行各自占据一个bank，这样每个线程读取不同行的数据时就不会发生bank conflict。padding的方式是在每行数据后面添加一定数量的空数据，使得每行的字节数恰好是bank的整数倍，这样就能保证每行数据各自占据一个bank。具体来说，当一个线程访问共享内存中的一个元素时，GPU会根据该元素在共享内存中的地址计算出其所在的bank，然后将该元素从该bank中读取出来。

非空流的操作可以被空流中的操作所堵塞，即kerne1执行完后，才会执行kernel 2和3.其中kernel1操作加入到非空流stream[0]中， kernel为2一个空流操作，kernel3为一个非空流stream[1]中，由于非空堵塞流中的操作会堵塞空流中的操作，所以kernel2会一直等待kernel1执行完毕之后，才会得到执行，而kernel2为空流，同时会堵塞kernel3,所以执行顺序为kernel1->kernel 2->kernel 3串行执行的效果。1：空流对非空流的堵塞行为。

【代码】CUDA-cudaEvent_t记录程序运行时间。