gpu上的矩阵乘法的设计与实现

《gpu上的矩阵乘法的设计与实现》由会员分享，可在线阅读，更多相关《gpu上的矩阵乘法的设计与实现（5页珍藏版）》请在金锄头文库上搜索。

1、 计 算 机 系 统 应 用 http:/www.c-s- 2011 年 第 20 卷 第 1 期 178 经验交流 Experiences Exchange GPU 上的矩阵乘法的设计与实现 梁娟娟，任开新，郭利财，刘燕君 (中国科学技术大学 计算机科学与技术学院，合肥 230027) 摘 要: 矩阵乘法是科学计算中最基本的操作， 高效实现矩阵乘法可以加速许多应用。 本文使用NVIDIA 的CUDA在 GPU 上实现了一个高效的矩阵乘法。测试结果表明，在 Geforce GTX 260 上，本文提出的矩阵乘法的速度是理论峰值的 97%，跟 CUBLAS 库中的矩阵乘法相当。 关键词: 矩阵乘

2、法；GPU；CUDA Design and Implementation of Matrix Multiplication on GPU LIANG Juan-Juan, REN Kai-Xin, GUO Li-Cai, LIU Yan-Jun (School of Computer Science and Technology, University of Science and Technology of China, Hefei 230027, China) Abstract: Matrix multiplication is a basic operation in scientifi

3、c computing. Efficient implementation of matrix multiplication can speed up many applications. In this paper, we implement an efficient matrix multiplication on GPU using NVIDIAs CUDA. The experiment shows that our implementation is as fast as the implementation in CUBLAS, and the speed of our imple

4、mentation can reach the peak speeds 97%, on Geforce GTX260. Keywords: matrix multiplication; GPU; CUDAGPU 是一种高性能的众核处理器，可以用来加速许多应用。 CUDA 是 NVIDIA 公司为 NVIDIA 的 GPU开发的一个并行计算架构和一门基于 C 的编程语言。在CUDA中程序可以直接操作数据而无需借助于图形系统的 API。现在已经有许多应用和典型算法使用CUDA 在 GPU 上实现出来。 1 引言 矩阵乘法是科学计算中的最基本的操作，在许多领域中有广泛的应用。对于矩阵乘法的研究有几个

5、方向。一个是研究矩阵乘法的计算复杂度，研究矩阵乘法的时间复杂度的下界，这方面的工作有 strassen算法1等。 另外一个方向是根据不同的处理器体系结构，将经典的矩阵乘法高效的实现出来，这方面的结果体现在许多高效的 BLAS 库。 许多高效的 BLAS库都根据体系结构的特点高效的实现了矩阵乘法，比如 GotoBLAS2, ATLAS3等。Fatahalian4等人使 用着色语言设计了在 GPU 上的矩阵乘法。CUBLAS库是使用 CUDA 实现的 BLAS 库，里面包含了高性能的矩阵乘法。 本文剩下的部分组织如下，第 2 节介绍了 CUDA的编程模型，简单描述了 CUDA 上编程的特点。第 3

6、节讨论了数据已经拷贝到显存上的矩阵乘法，首先根据矩阵分块的公式给出了一个朴素的矩阵乘法实现，分析朴素的矩阵乘法的资源利用情况，然后提出了一种新的高效的矩阵乘法。第 4 节讨论了大规模的矩阵乘法的设计和实现， 着重讨论了数据在显存中的调度。第 5 节是实验结果。第 6 节是总结和展望。 2 CUDA编程模型和矩阵乘法回顾 2.1 CUDA 编程模型 NVIDIA的GPU是由N个多核处理器和一块显存构成的。每个多核处理器由 M 个处理器核，1 个指令部件，一个非常大的寄存器堆，一小块片上的共享内 基金项目:国家自然科学基金(60833004);国家高技术研究发展计划(863)(2008AA0109

7、02) 收稿时间:2010-04-26;收到修改稿时间:2010-05-21 2011 年 第 20 卷 第 1 期 http:/www.c-s- 计 算 机 系 统 应 用 Experiences Exchange 经验交流 179存以及一些只读的缓存(包括纹理缓存和常量缓存)。 多处理器内所有核共用一个指令部件，在一个时钟周期内，每个处理器核执行同一条指令，但是数据可以是不一样的。 GT200是NVIDIA的第10代GPU， 它发布于2008年 6 月。 GT200 的多核处理器包含 8 个处理器核， 16K字节的片上共享存储器，16384 个 32 位的寄存器，1个指令部件，1 个特殊功

8、能单元部件，1 个双精度部件以及一些只读的缓存。GT200 是第一代支持双精度计算的显卡。 在 CUDA 中的核心概念是内核函数。内核函数是从 CPU 端调用，运行在 GPU 上的函数。CUDA 的运行管理系统产生一个 grid 来运行一个内核函数。1 个grid 包含一组 block，每个 block 是有一组线程组成。grid 维度和 block 的维度都是在调用内核函数的时候指定的。 在一个 grid 中， block 之间没有同步和通信机制， 因此 block 之间必须是独立的。 在一个 block 里面，线程之间可以用同步原语_syncthreads()来同步，每个block 可以分

9、配一定量的所有线程都可以访问的共享内存。一个 block 只能映射到一个多核处理器上执行，而不能将线程分布在不同的处理器中。一个多核处理器上可以同时运行多个 block。在运行一个 block 的线程时，CUDA 运行管理系统将 block 中的线程分成一个个 warp，每次取出一个 warp 来运行。一个 warp 分成前半个 warp 和后半个 warp。运行的最小单位是半个 warp，一个周期里面，半个 warp 的所有线程执行相同的指令。在 CUDA 2.3 的版本中，一个 warp 中含有 32 个线程。 一个内核函数可以使用的存储器资源包括全局内存，共享内存，纹理内存，常量内存以及

10、寄存器。全局内存， 纹理内存以及常量内存都是分配在显存中的，可以被所有的 block 访问。纹理内存和常量内存都是只读的，如果重复访问的话，它们会被缓存在片上的缓存上。全局内存是不会被缓存的。GPU 访问显存是使用内存事务来完成的。完成一个内存事务的代价是非常高的，通常需要 400-600 个周期。但是，如果半个 warp 的所有线程的访问满足：(1)所有线程访问 1个字节， 并且访问的位置在一个 32 字节对齐的同一个32 字节的段里面。(2)所有的线程访问 2 个字节，并且访问的位置在一个 64 字节对齐的同一个 64 字节的段里面。(3)所有的线程访问 4 个字节，并且访问的位置在一个

11、128 字节对齐的同一个 128 字节的段里面。三个条件中的一个， 则半个 warp 的访存会被合并成一个内存事务。 2.2 矩阵乘法回顾 给定矩阵乘法 C=AB，其中( )ijCc=是MN的矩阵，( )ijAa=是MK的矩阵，( )ijBb=是KN的矩阵。则ijc的定义为 1Kijikkj kca b= (1) 假定,Mwp Nhq Kkr=，将矩阵 C 分解为pq的分块矩阵( )ijc，每个ijc是一个wh的小矩阵，A 分解成为pk的分块矩阵( )ija，每个ija是一个wr的矩阵，B 分解成为kq的分块矩阵( )ijb，每个ijb是一个rh的元素。则分块矩阵乘法的定义为 1kijissj

12、 sca b= (2) 在等式中核心的操作是 ijissjijca bc=+ (3) 假定缓存中同时可以放下分块后的ijc，isa，sjb，则等式需要的访存数量是whwKKh+，计算量是2whK, 平 均 一 次 访 存 的 计 算 量 是2/ ()whKwhwKKh+，略去分母中的小项wh得到2/ (1/1/ )wh+，从而知道1/1/wh+越小，平均一次访存做的计算量是越多的。 3 算法设计和实现 本节假定矩阵 A， B， C 已经存放在显卡的显存中，矩 阵 是 以 行 主 序 的 方 式 存 储 的 。 假 定,Mwp Nhq Kkr=， A， B， C 分别是MK，KN，MN的矩阵。

13、3.1 普通实现 按照等式，每个ijc的计算是相互独立的。将矩阵C 分成许多wh的子矩阵，每个子矩阵使用不同的block 来计算。假定一个 block 的线程数目是T，这样每个线程需要计算 C 中/wh T个元素。将C中的元素固定分配给某个线程，这样可以将C中的元素存储在每个线程的寄存器中。为了充分利用 GPU 上的带宽，将矩阵A和B读取进来的数据先缓存在片上的共享存储器中。每次读取A的一个wr的子矩阵和B的一个rh的子矩阵到共享存储里，在共享存储器中完成一个小块的矩阵乘法。这样得到的算法如图 1 所示。 计 算 机 系 统 应 用 http:/www.c-s- 2011 年 第20卷 第 1

14、期 180 经验交流 Experiences Exchange 下面讨论算法的参数。为了提高访问全局内存的效率，r 和 h 必须是 16 的倍数。同时，一个 block 需要的共享内存是8()wrrh+字节。每个 block 需要的共享内存的大小必须小于一个 SM 上的共享存储器的大小，从而得到8()16384wrrh+，同时所有线程需要的寄存器数目要小于一个 SM 上的寄存器数目。每个双精度的数需要占用 2 个 32 位的寄存器， w 和 h 必须满足不等式216384wh 。为了让 w 和 h 足够大，我们取16r =，这样有128wh+。根据矩阵分块乘法，1/1/wh+越小越好，而 h

15、必须是 16 的倍数，在算法中取32w =，64h =。于是有2048wh =。而当前一个block 中的线程数目最多是 512，于是可以取 T=512，每个线程计算 4 个数据。使用 cuda 编译器编译以后，每个线程需要使用 24 个寄存器，一个 block 总共需要12288 个寄存器，需要 12728 个共享存储器。 图 1 矩阵乘法的普通实现 3.2 优化实现 在 3.1 节的算法中，矩阵 A 和 B 的子块都是保存在共享存储器中，而 C 的子块是保存在寄存器中。主要是根据共享存储器的限制来得到 w 和 h 的值，寄存器的使用数目远远没有达到处理器的上限，处理器的寄存器资源没有被充分利用。为了充分利用处理器上的寄存器资源，我们改变数据映射方式。仍然假设每个 block 计算wh的一个子矩阵Cs，使用一个线程来计算Cs中的一列，也就是要计算 w 个元素。注意到一个线程只会访问矩阵 B 的一列， 也就是要计算 w 个元素。注意到一个线程只会访问矩阵 B 的一列，不同的线程访问 B 的元素是不相交的。这样 B 的元素可以放在线程中的寄存器中，而不用存放在共享存储器中。在共享存储器中分配一个wr的数组As存放矩阵 A的子块。这样得到的算法如图 2 所示。 图 2 矩阵乘法的优化实现 下面讨论算法的参数。主要考虑的是两个条件，一个是充分利用片上的存储资源，包括寄存器和共享存储器