基于MPI+CUDA异步模型的并行矩阵乘法

基于MPI+CUDA异步模型的并行矩阵乘法
刘青昆;马名威;阎慰椿
【摘要】Matrix multiplication plays an important role in scientific computing. Different structural models can improve the performance of parallel matrix multiplication. In the existing MPI + CUDA synchronization model, the host-side need enter the waiting state and cannot continue to work until the device completes the task, which obviously wastes time. Concerning this question, a parallel matrix multiplication based on MPI + CUDA asynchronous model was proposed. This model prevented host-side's entering into the waiting state, and used CUDA-stream technology to solve the problem of data bulk over GPU memory. By analyzing the speedup ratio and efficiency of the asynchronous model, the experimental results show that MPI + CUDA parallel programming obviously promotes parallel efficiency and large-scale matrix multiplication' s speed, which exerts the advantages of the distributional memory between the nodes and the share memory in the node. It is an effective and feasible parallel strategy.%矩阵乘法在科学计算领域中起着重要的作用,不同结构模型能够改善并行矩阵乘的性能.现有的MPI+CUDA同步模型中,主机端需要进入等待状态,直到设备端完成任务后才能继续工作,这显然浪费时间.针对上述问题,提出一种基于MPI+ CUDA异步模型的并行矩阵乘法.该模型避免了主机端进入等待状态,并采用CUDA 流技术解决数据量超过GPU内存问题.通过分析异步模型的加速比和效率,实验结果表明,此方法显著提高了并行效率和大型矩阵乘法的运算速度,充分发挥了节点间分布式存储和节点内共享内存的优势,是一种有效可行的并行策略.
【期刊名称】《计算机应用》
【年(卷),期】2011(031)012
【总页数】4页(P3327-3330)
【关键词】矩阵乘法;并行计算;混合编程;消息传递接口;统一计算设备架构
【作者】刘青昆;马名威;阎慰椿
【作者单位】辽宁师范大学计算机与信息技术学院,辽宁大连 116081;辽宁师范大学计算机与信息技术学院,辽宁大连 116081;辽宁师范大学计算机与信息技术学院,辽宁大连 116081
【正文语种】中文
【中图分类】TP311;TP316.4
0 引言
矩阵乘法在科学计算领域起着重要的作用。

由于矩阵乘法计算量大，有效矩阵乘并行算法在实际应用中就显得非常重要。

在国内，刘杰等人于1997年提出了一种B 迁移矩阵乘并行算法[1];吴建平等人于1999年将cannon算法由正方形网格处理机阵列推广到长方形网格处理机阵列[2]。

在国外，文献[3]提出的基于OTIS－Mesh网络结构的并行矩阵乘算法;文献[4]提出六种不同的并行矩阵乘法在非均质多核处理器计算集群的实验研究。

这些成果表明，根据不同的结构模型开发并行矩阵乘法能够有效地改善并行算法的性能。

随着计算机图形处理器(Graphic Processing Unit，GPU)的计算能力越来越强大[5－6]，它已经具备了实现大规模快速计算的编程能力，2007年NVIDIA公司提
出的计算统一设备架构(Computer Unified Device Architecture，CUDA)技术就是这方面的杰出代表[7]。

目前普遍使用的并行程序设计模型是消息传递的消息传
递接口 (Message Passing Interface，MPI)模型，这种消息传递模型一般面向分
布式内存的，不能充分发挥集群节点内共享内存的优势。

本文将消息传递模型
MPI和共享内存模型CUDA结合在一起，开发一套MPI+CUDA异步模型，该异
步模型主机端不需等待设备端完成任务就可以进行数据通信，使CPU与GPU执
行任务期间互不影响，更加充分地发挥了消息传递和共享内存两种编程模型的优势。

1 MPI+CUDA并行模型
MPI[8]是消息传递接口的简称，是一种消息传递编程模型。

由于MPI细粒度并行
会引发大量的通信，从而产生昂贵的通信代价，造成较大开销。

通过开发混合并行程序，MPI只需要建立粗粒度并行，实现节点间通信即可，有效地克服了MPI通
信代价大的缺点。

CUDA是基于GPU进行通用计算的开发平台，适合于大规模的并行数据计算。

在CUDA编程模型中，GPU被视为一个协处理器，能够执行大量线程并行。

在其架
构下面，一个程序分为两部分:host端和device端[9]。

host端是指在CPU上执
行的部分，而device端是在显示芯片上执行的部分。

当处理数据并行计算时，CUDA首先将程序编译成GPU能够执行的程序，而后传送到GPU。

运行在GPU
上的CUDA并行计算函数叫作“kernel(内核函数)”。

在实际运行中，CUDA会
产生大量在GPU上执行的线程，线程间只能通过共享存储器和同步机制进行通信，能充分发挥节点内共享内存的优势，非常适合MPI+CUDA模型中细粒度并行部分。

混合编程模型的优势在于能取单一模型所具备的长处，克服单一模型所存在的缺点与不足。

在编写并行程序时，使用MPI将任务分解成若干子任务，每个节点完成
一个子任务，由于每个节点只能读写本地内存中的数据，所以各子任务之间通过消
息传递的方式来进行通信完成对远程数据的访问。

每个节点在GPU上采用多线程共享存储机制来实现数据共享进行细粒度并行，充分发挥单个节点计算能力，从而提高计算性能和加速比。

混合粒度的并行化可以兼顾性能和可扩展性，MPI解决
多处理器间的粗粒度的并行通信，而GPU采用由硬件管理的轻量级线程并行机制，可以很好地解决处理器内部的通信，CUDA提供大量线程数目可以快速地进行大
规模矩阵乘计算。

在MPI+CUDA并行模型中，MPI并行位于顶层，CUDA位于
底层。

这种模型很好地映射了多处理器计算机组成的集群体系结构[10]。

文献[11]提出了MPI+CUDA并行编程同步的计算模型，如图l所示，处理器CPU 将分配的矩阵块数据传给GPU，CUDA就会产生许多在GPU上执行的线程来实
现块矩阵乘加运算，然后将计算结果从GPU调回CPU上，而后进行进程间的数
据通信。

图1 MPI+CUDA同步模型
很显然，这种调用同步版本的GPU函数在设备完成请求任务前，都不会返回主机线程，而主机端线程将进入让步(yield)、阻滞(block)或者自旋(spin)状态，这就浪费了时间。

如果当GPU进行计算或数据传输时就返回主机线程，那么主机线程不必等待GPU运行完毕就可以继续进行计算，从而节省了处理器等待时间，因此我们提出了一种MPI+CUDA异步模型。

该模型使CPU和GPU可以进行并行工作，有效提高了并行效率。

MPI+CUDA异步模型显示如图2。

异步模型是指GPU进行的操作(内核启动或异步存储器拷贝数据)从主机端启动后，在GPU真正完成这些操作之前，在主机端就可以得到这些函数的返回值，CPU线程就可以继续执行下一步操作。

换句话说，CPU可以在GPU端进行运算或者数据传输的同时进行其他操作，这样更加有效地利用了系统中的计算资源。

2 并行矩阵乘法
2.1 矩阵映射策略
已知A、B和C是n阶矩阵，矩阵乘法[12]可表示为:C=A×B。

这里A和B为进
行乘法运算的初始矩阵，C为存放结果的矩阵，C在计算前为零矩阵。

假设使用
p2个处理器，将其修改成二维网格逻辑列阵p×p形式。

在处理器Pij存储矩阵元素或分块矩阵这一部分，矩阵元素或块矩阵乘法要计算所有匹配的Aij和Bij，然
而只有在主对角线的进程才是匹配的。

因此，矩阵映射策略可以直接将匹配的矩阵元素或者块矩阵映射到二维网格处理器中。

若p=n时，矩阵映射策略可以直接根
据二维网格中处理器和矩阵元素进行一一映射，例如矩阵元素a01和b11被映射到处理器P01上。

同理，结果矩阵也用相同方法存储到处理器中。

若进行大规模
矩阵运算(p≪n)时，显然这种情况不能简单地进行处理器和矩阵元素的一一映射。

这时，可以将矩阵A，B和C分成p2个块矩阵Aij、Bij和Cij(0≤i，j＜p)，每个
矩阵块的大小为(n/p)×(n/p)，然后将它们按照上述映射策略把矩阵块分配给各处
理器。

图2 MPI+CUDA异步模型
2.2 算法实现
在MPI+CUDA异步模型中进行两个n×n矩阵的并行乘法运算时，首先，根据上
述策略对矩阵A、B和C进行映射，存储到二维处理器网格中。

而后操作如下。

步骤1 各处理器启动CUDA并行计算Cij=Aij×Bij。

Pij在GPU上实现乘－加运算的同时，CPU将矩阵块Aij向左循环移动1个位置，矩阵块Bij向上循环移动1个位置。

步骤2 块Aij和Bij的循环单步移位操作完成后，如果CUDA已经将计算结果传回CPU，那么执行步骤1，如果CUDA仍在计算，那么CPU等待CUDA计算结束，传回计算结果后再执行步骤1。

重复执行步骤1p次，每次移动后各处理器启动CUDA执行一次乘—加运算。

在整个算法执行过程，首先，通过映射策略，利用MPI消息传递将矩阵A中的元
素或矩阵块Aik和矩阵B中的元素或矩阵块Bkj分别存储在二维网格处理器Pij中。

而步骤1执行乘加运算这部分是调用kernel完成的。

kernel先开辟大小合适的矩阵块a和b的数据空间后，将Aik和Bkj的数据从CPU调入GPU中，CUDA就
会产生许多在GPU上执行的线程来实现块矩阵a与b相乘，然后将计算结果从GPU调回CPU上。

而CPU将数据传给GPU后，不需要等待CUDA完成计算结果，就可以利用MPI消息传递的发送和接收功能将矩阵A块和B块分别执行循环左移和循环上移单步移位操作。

为了避免在通信过程中发生死锁，我们错开奇数号及偶数号处理器的收发顺序，利用MPI消息传递功能，使偶数号处理器先发送后
接收;而奇数号处理器先将A的矩阵块存于缓冲区buffer中，然后接收编号在其后面的处理器发送的A矩阵块，最后再将缓冲区中原矩阵A的矩阵块发送给编号在
其前面的处理器。

再利用二维网格奇偶行数的收发顺序，使所在二维网格的偶数行的处理器先发送后接收，而所在二维网格的奇数行的处理器先将B矩阵块存于缓
冲区buffer中，然后接收编号在其二维网格下面的处理器发送的B矩阵块，最后
再将缓冲区中原矩阵B的矩阵块发送给编号在其二维网格上面的处理器。

移位操
作结束后，CUDA完成任务返回结果，各处理器再次调入kernel实现Cij=Aij×Bij，总共需执行p次乘加运算和p次块Aij和Bij的循环单步移位。

经过这些步骤求出来的矩阵C元素或矩阵块Cij存储在处理器Pij中，最后经过整理，就得到全部的
结果矩阵C了。

需要注意的是，为了提高CPU与GPU的数据传输宽带，我们使用页锁定主机内
存(cudaHostAlloc)来分配主机内存，使该内存始终停留在物理内存中，不会被操作系统分页从而交换到磁盘上，防止了该内存被破坏或重新定位，有效地提高了主机与设备端之间的数据传输速率。

由于GPU的内存远少于主机内存，当计算数据量超过GPU内存时，整个缓冲区
无法一次性地填充到GPU。

我们提出使用CUDA流思想将输入缓冲区划分更小的
块，将其中一部分复制到GPU，运行内核函数，然后将输出缓冲区中的这部分结
果复制回主机。

CUDA流相当于GPU操作队列，并按指定顺序执行。

例如将核函数启动、内存复制等操作添加到流的顺序作为执行顺序。

每个流可视为GPU上的一个任务，这些任务是可以并行执行的。

在内核函数调用的尖括号中带有一个流参数，这样内核函数调用将是异步的。

从顺序上来说，第一次放入流中的复制操作将在第二次复制操作之前执行，第二次复制操作将在内核函数启动之前完成，而内核函数将在第三次复制操作开始前完成。

CudaMemcpyAsync函数调用来完成内存数组数据的复制。

程序将大规模矩阵乘分成很多小矩阵块相乘。

GPU中使用线程块block，把两个小块矩阵加载到共享
内存上，使得小矩阵相乘不再需要存取外部的内存。

由于GPU上的内存控制器从某个固定倍数地址开始读取，才会有最高的效率。

因此根据其支持的线程数，最好的办法就是将矩阵分解为多个16×16的方阵，每个线程块使用256个线程并行计算。

为了方便线程块的计算，每个block分配16×16个thread，而后建立了(WC/16)×(HC/16)个block(WC和HC分别是结果矩阵C的列数和行数)，这样划分充分利用了GPU的多线程特点，从而能大大提高运算速度。

2.3 算法性能分析
从算法中可以知道，循环体中有一次处理器间的数据通信操作和一次四则运算，设ts为启动时间，处理器平均传输每个数据的时间为tw，一次循环操作传输的数据
大小为n2/p2，则传输数据所需要的时间为ts+twn2/p2;那么p次循环所需时间
是tsp+twn2/p。

又假设一次乘法和加法运算时间为一个单位时间，由于每个处理器要进行p次(n/p)×(n/p)子矩阵的乘法和加法运算，因此每个处理器花在计算中的总时间为p×(n/p)3=n3/p。

故其总执行时间为:
Tp=tsp+twn2/p+n3/p。

在MPI+CUDA异步模型中，CUDA利用GPU多线程计算矩阵乘—加运算的同时，
处理器则利用MPI消息传递进行数据的交换，将通信开销隐藏在计算开销中，大
大提高了并行效率，图3显示了MPI+CUDA模型同步与异步执行方式比较。

显然，不使用异步执行时的程序运行时间为(初始化数据时间+计算时间+通信时间);而使用异步后，时间将为(初始化数据时间+计算时间)。

此时，CPU和GPU间的数据传输时间被有效地隐藏了，程序的性能得到了改善。

图3 MPI+CUDA模型同步与异步执行方式比较
3 数值实验
为了测试这种编程模型的性能，将其对同步模型和MPI模型、CUDA模型进行模
拟数值实验，实验环境是:4台带有支持CUDA的GeForce GTX 470显卡的多核
处理器。

采用大规模矩阵相乘作为测试程序。

机器配置如下:Intel core 2 CPU
**********GHz，1.00 GB 内存，OpenSuse11.2 操作系统，MPICH2函数库，CUDA，通过千兆交换机连接成集群系统。

通过编写MPI程序控制这些处理器间的通信实现MPI模型和MPI+CUDA模型的逻辑连接。

首先，在固定处理器数目的条件下(P=4)选取五组固定矩阵规模，对本文算法和同步算法以及MPI算法、CUDA算法在执行时间方面进行了实验测试，
运行时间如表1所示，根据时间得到不同模型的加速比如表2所示，再由加速比
计算出的效率如表3所示。

而后还在单台PC上对传统O(n3)的串行矩阵乘算法在执行时间方面进行了实验测试，得出实验数据如图4所示。

表1 模型运行时间表 s矩阵阶数模型MPI+CUDA同步 MPI+CUDA异步MPI CUDA 4800 5.43 5.08 406.46 6.03 6400 10.38 8.96 956.52 11.64 8000 16.79 13.91 1915.08 20.89 9600 24.95 20.04 3322.55 33.28 12800 47.63 35.77 8718.88 68.07
表2 加速比实验值矩阵阶数模型MPI+CUDA同步 MPI+CUDA异步MPI CUDA 4800 251.84 281.54 3.52 237.09 6400 360.45 417.92 3.91 321.52 8000
297.70 442.54 3.88 355.60 9600 527.34 656.51 3.96 395.39 12800 1563.34 2081.75 8.54 1093.94
表3 表效率实验值矩阵阶数模型MPI+CUDA同步 MPI+CUDA异步MPI CUDA 4800 62.96 70.38 0.880 59.27 6400 90.11 104.48 0.977 80.38 8000 74.43 110.63 0.970 88.90 9600 131.83 164.12 0.990 98.85 12800 390.83 520.43
2.130 27
3.48
从图4上看，MPI+CUDA模型的执行时间远远小于MPI的并行执行时间，表现
出非常高的执行效率，相对于同步模型而言，MPI+CUDA异步模型执行时间更少，执行效率的提高更加明显。

图4 执行时间的比较
由表1和表2的实验数据可知，随着矩阵阶数的增大，本文算法加速比呈线性增长，效率越来越高，充分证明该模型具有很高的可行性。

4 结语
实现MPI+CUDA异步模型，通过对矩阵乘法的测试，实验证明，MPI+CUDA异步模型具有更高的执行效率，更加充分发挥节点间分布式存储和节点内共享存储的优势。

在多核处理器环境下该模型可以有效提高并行计算性能，是一种高效可行的并行编程策略。

参考文献:
[1] 刘杰，胡庆丰，迟利华.并行矩阵乘的B迁移算法[J].计算机工程与科学，1997，19(3):50 －54.
[2] 吴建平，迟学斌.分布式系统上并行矩阵乘法[J].计算数学，1999，21(1):99 －108.
[3] WANG C F， SAHNI S.Matrix multiplication on the OTIS－mesh optoelectronic computer[J].IEEE Transactions on Computers，2001，
50(7):635－645.
[4] ALONSO P，REDDY R，LASTOVETSKY A.Experimental study of six different implementations of parallel matrix multiplication on heterogeneous computational clusters of multicore processors[C]//18th Euromicro Conference on Parallel， Distributed and Network－based Processing.Washington， DC:IEEE Computer Society，2010:263 －270. [5] 刘伟峰，杨权一，曹邦功，等.基于 GPU的高度并行Marching Cubes改进算法[J].微电子学与计算机，2008，25(9):151 －157.
[6] 钱悦.图形处理器CUDA编程模型的应用研究[J].计算机与数字工程，2008，36(12):177 －180.
[7] 董荦，葛万成，陈康力CUDA并行计算的应用研究[J].信息技术，2010(4):11－13.
[8] GROPP W，LUSK E，DOSS N， et al.A high－performance， portable implementation of the MPI message－passing interface standard[J].Parallel Computing，1996，22(6):789 －828.
[9] NVIDIA Corporation.CUDA Programming Guide 2.3[M].[s.l.]:NVDIA Coporation，2009:11 －13.
[10] 赵永华，迟学斌.基于SMP集群的MPI+OpenMP混合编程模型及有效实现[J].微电子学与计算机，2005，22(10):7 －11.
[11] 许彦芹，陈庆奎.基于SMP集群的MPI+CUDA模型的研究与实现[J].计算机工程与设计，2010，31(15):3408 －3412.
[12] 夏丽华，张丽晓.PC集群环境下大规模矩阵乘法算法的研究[J].计算机工程与设计，2007，28(19):4603 －4605，4615.。