高性能计算 第二次作业

张钧宇
18340211
计算机科学与技术
zhangjunyu@nelson-cheung.cn

 

运行方法

进入到code/build文件夹下,输入如下命令。

code/build会生成4个可执行文件,名字形如test_*.exe,如下所示。

运行test_p2p.exetest_set.exe时,在mpiexec加上运行的进程数量即可,例如使用4个进程运行test_p2p.exe

结果如下。

1-2

修改src/mpi_p2p.cppsrc/mpi_set.cpp中的如下语句即可修改矩阵的规模。

其中,矩阵乘法为C=ABA的规模为REAL_M x REAL_KB的规模为REAL_K x REAL_N

硬件数据

运行环境:VMware 16.x虚拟机

操作系统:Ubuntu 18.04

内存:8G

CPU:6个逻辑CPU

通用矩阵乘法

对于一个M×K的矩阵AK×N的矩阵B之间的乘积,其定义为

C=AB

其中

Ci,j=kAi,k×Bk,j1iM, 1jN,1kK

Lab1中已经实现了通用矩阵乘法和相应的优化方法,取其中性能最好的串行版本作为MPI实现方法的比较对象。对于不同规模的矩阵,取定测试矩阵都为N×N的方阵,矩阵的每一个元素从正态分布N(0,12)​中随机生成,则串行版本的性能测试结果如下。

51210242048
0.027s0.209s1.563s

测试代码放置在src/ijk_comp.cpp下。

矩阵乘法库函数

首先在include/matrix.h下声明两个新函数。

gen_rand_mat用于生成一个MxN的矩阵,然后按照行优先的方式放置在mat中。

mat_mul用于计算MxK的矩阵aKxN的矩阵b的乘积,结果放置在c中。

注意到这里的矩阵虽然是二维的,但是都统一使用了行优先的方式存储在一维数组中。这是因为使用二维数组来表示一个二维矩阵虽然比较直观,但是不利于使用MPI进行传输,所以在后面使用MPI来实现矩阵乘法中,矩阵的表示都是一维数组。这里虽然没有涉及到MPI,但为了和后面的实现兼容,也使用了一维数组来表示二维矩阵。

声明了函数后,在src/matrix.cpp中实现这两个函数,如下所示。

函数的实现使用了Lab1中实现的矩阵类Matrix的成员和方法,可以理解为是在Matrix上的一层封装。

实现完成后,将函数的实现编译成动态链接文件libmatrix.so,对应的makefile如下。

编写测试代码test/test.cpp

将测试代码编译成可执行文件test_so

运行可执行文件test_so

可以看到输出了结果矩阵。

1-1

test_so的编译命令中并没有使用到gen_rand_matmat_mul的实现,而是通过-L.-lmatrix参数使用了动态链接库libmatrix.so

注意-lmatrix参数必须放置在源代码文件名后面,否则即使动态链接库libmatrix.so存在,也会报找不到函数实现的错误。

点对点通信的矩阵乘法优化

基本思路

使用点对点方式通信的方式可以理解为P2P​的方式,即在计算过程中,每一个进程是同等的计算节点,计算过程和计算量也相同。

在给出MPI点对点通信的矩阵乘法算法前,首先假设有p个进程执行矩阵乘法优化,矩阵乘法为

C=AB

其中,A的规模为M×KB的规模为K×N

由于每一个进程都是独立和等同的计算节点,借助于矩阵分块乘法法则,可以将矩阵C划分为若干独立的块,每一个进程执行其中一块的计算,最后由进程0将每一个进程的计算结果聚合成最后的矩阵C​。这就是点对点通信实现矩阵乘法的基本思路,详细算法如下。

首先,将矩阵A按行的方式划分为p块。

A=[A0A1Ap1]

将矩阵B按列的方式划分为p块。

B=[B0B1Bp1]

注意到MN不一定整除p,假设有

M=mp+pm, 0pm<pN=np+pn, 0pn<p

此时,只要令A0的大小为(m+pm)×KA1,,Ap1的大小为m×K,即可解决在划分的时候M不一定能够整除p的问题。同理,在划分的时候令B0的大小为K×(n+pn)B1,,Bp1的大小为K×n

矩阵乘法可化为

C=AB=[A0A1Ap1][B0B1Bp1]=[A0B0A0B1A0Bp1A1B0A1B1A1Bp1Ap1B0Ap1B1Ap1Bp1]

此时,只要令进程i计算第i行的矩阵即可,即

[AiB0AiB1AiBp1]

上述矩阵可化为

Ai[B0B1Bp1]=AiB

在计算之前可以将A0,A1,,Ap1​​​分别发送给进程P0,P1,,Pp1​​​。但是,每一个进程的计算都需要用到矩阵B​​​。若每一个进程都存储矩阵B​的副本,则在计算过程中的空间开销较大。为了减小矩阵存储的开销,可以先将B0,B1,,Bp1​分别发送给进程P0,P1,,Pp1​。每一个进程只存储矩阵B​的其中一块。计算完成后,每一个进程i将自己存储的子矩阵Bi发送给序号+1的进程,即进程0发送给进程1,进程1发送给进程2,如此类推。特别地,进程p1发送给进程0​。

如此重复p​​​轮后,便完成了矩阵乘法的计算。

可以将上述算法总结如下。

  1. 进程0读入两个矩阵A,B,将矩阵A按行的方式分成p块,其中A0的规模为(m+pm)×KA1,,Ap1的大小为m×K;将矩阵B按列的方式分成p块,其中B0的规模为K×(n+pn)B1,,Bp1的大小为K×n
  2. 进程0分别将A0,A1,,Ap1​发送给进程P0,P1,,Pp1​,然后将B0,B1,,Bp1​发送给进程P0,P1,,Pp1​。
  3. 假设P0,P1,,Pp1​存储的子矩阵B为Bk0,Bk1,,Bkp1​,进程Pi​计算矩阵乘法AiBki​。
  4. 计算完成后,矩阵Pi将矩阵Bki发送给Pi+1。特别地,Pp1发送给P0
  5. 重复2-4步p次。
  6. P1,P2,,Pp1将计算结果发送给P0P0输出结果。

代码实现

实现代码放置在src/mpi_p2p.cppmain函数中,下面是对main函数各部分的解析。

首先,对MPI进程初始化,获取进程号和进程数。

进程0确定矩阵规模,并采用点对点的方式将矩阵的规模发送给其他进程。

根据前面的算法,每一个进程每次只会计算一个A的子矩阵和B的子矩阵相乘的结果。所以每一个进程需要分配的空间包含了存储A​的子矩阵local_aB的子矩阵local_b和计算结果C的子矩阵local_c的空间。

由于M,N都不一定能整除p​,即

M=mp+pm, 0pm<pN=np+pn, 0pn<p

在最开始划分的时候,B0​的规模可能比B1,B2,,Bp1​的要大。而且,根据前面的算法,每一个进程i​都需要计算AiB0​的值。为了存下的所有子矩阵和对每一个进程屏蔽进程规模的差异,在分配空间的时候,local_alocal_blocal_c实际上是按照存储A,B的子矩阵需要的最大规模来分配的。

然后,进程0生成随机矩阵A,B

对于不同的进程,这里会产生不同的行为。

对于进程0,由于其要存储矩阵A,B和计算结果C,在分配空间的时候需要给abc分配和矩阵A,B,C规模相同大小的空间。然后,进程0生成随机矩阵A,B

generate_matrix利用了矩阵类Matrix的成员函数生成了两个随机矩阵,然后按行优先的方式将二维矩阵存储在一维数组中。

对于其他进程,由于其只需要存储最终结果的部分元素,所以只分配了对应的空间。

由于不能整除的问题导致了B​划分出来的子矩阵规模并不相同,并且每一个进程在计算的时候都会用到B​​的所有子矩阵。为了屏蔽子矩阵的规模差异和减少判断语句采用统一的计算过程,在传输矩阵之前,需要先传送矩阵的规模,然后再传送矩阵的元素。由于MPI保证了同一进程发送出的消息的有序性,矩阵的规模必然会在矩阵的元素发送之前收到。此时,定义一个结构体MatrixType来表示矩阵的行和列。

通过mpi_type_create_struct来将上述结构体定义为MPI函数可以识别和处理的MPI_Datatype

进程0Ai,Bi发送给进程i

发送的方式是先发送矩阵的规模,然后再发送数据。子矩阵的获取函数get_submatrix实现如下。

收到矩阵后,每一个进程执行相同的计算过程,计算子矩阵乘积,然后将矩阵B的子矩阵发送给下一进程,接收上一进程发送来的B的子矩阵,如此重复p轮。p是用于计算的进程数量。

由于MPI_SendMPI_Recv有可能会发生阻塞,从而导致死锁。为了解决这一问题,这里让进程号为奇数的进程先接收然后发送,进程号为偶数的进程先发送再接收。

矩阵的计算函数calculate_matrix实现如下。

实现的基本方式是Lab1中的内存重排实现方式。由于二维矩阵采用类行优先的方式存储在一维数组中。对于一个M×N的矩阵,第i行第j列的元素在一维数据的索引为i×N+j

在每一轮的计算中,每个进程只会计算两个子矩阵的乘积。计算完成后,需要放到进程最终计算结果中,实现函数set_submatrix如下。

注意到矩阵B​​的每一个子矩阵在每一轮的计算中都是有序地传递,即使没有传输子矩阵的下标,每一个进程都可以根据进程号和计算轮数计算出来。

计算完成后,进程0收集其他进程的计算结果生成最终的结果。

点对点通信的实现大致如上所示。

集合通信的矩阵乘法优化

基本思路

和点对点的实现方式不同,集合通信采用的是主从模式。在集合通信的实现方式中,进程0充当主节点,其他进程充当从节点。主节点负责将数据分发到从节点,从节点接收到数据后进行计算,最后主节点收集从节点的计算结果后得到最终的结果。粗略地来看,集合通信的实现方式和点对点的实现方式基本相同。但是,在点对点的实现方式中,每一个进程都是对等的,数据可以在进程中相互流动。但是在集合通信的方式中,才用了主从模式,从而规定数据只能从主节点流向从节点,从节点之间是不会发生数据传输的。

集合通信的实现方式同样利用了矩阵分块乘法的思想。假设有p个进程,对于矩阵乘法

C=AB

其中,矩阵A的规模是M×K,矩阵B的规模是K×N

可以按列将矩阵A​​按行的方式分成p​块

A=[A0A1Ap1]

矩阵乘法可化为

C=AB=[A0A1Ap1]B=[A0BA1BAp1B]

只要令进程i​计算AiB​​​即可将矩阵乘法并行化。

在对矩阵A​进行划分的时候,依然会存在M​不能整除p​​的情况,即

M=mp+pm, 0<pm<p

这里的解决方法和点对点通信的实现方法不同。在集合通信的划分中,要求A0,A1,,Ap1的规模均为m×K。若pm>0,则会存在一个规模为pm×K的矩阵Ap没有归纳到划分中。注意到pm应该是远小于m​的,此时令进程0在最后计算ApB并加入到结果即可。

由此可以得到集合通信的矩阵乘法计算方法。

  1. 进程0读入两个矩阵A,B,将矩阵A按行的方式分成p块,其中,A0,A1,,Ap1的规模为m×K。​
  2. 进程0广播矩阵B​到其他进程,将A​散射(Scatter)给所有的进程。
  3. 进程i计算AiB​的结果。
  4. 进程0回收收集(Gather)其他进程的计算结果。
  5. 若存在M不能整除p的情况,进程0计算ApB并加入到最终结果。

代码实现

实现代码放置在src/mpi_set.cppmain函数中,下面是对main函数各部分的解析。

首先,初始化MPI,得到进程号和进程数量。

进程0读入矩阵A,B的规模,然后使用MPI_Bcast广播给其他进程。

分配存储用于计算的空间。

根据前面的算法,进程i负责计算AiB。因此,每一个进程需要分配的空间包含了存储矩阵Ai的空间local_a,存储矩阵B的空间b和存储计算结果的空间local_c

进程0生成随机矩阵A,B

其中,a用于存储矩阵Ac用于存储最后的结果,即矩阵C。生成矩阵的函数generate_matrix和点对点通信的方法相同。

进程0利用函数MPI_Bcast将矩阵B广播给其他进程。

进程0利用函数MPI_Scatter将矩阵A​散射(Scatter)给其他进程。

进程i计算AiB的结果。

进程0利用函数收集(Gather)其他进程的计算结果。

若存在不能整除的情况,即在对矩阵A进行划分时会多余出一个矩阵Ap没有被纳入到计算过程中,则由进程0负责计算ApB并加入到最终结果中。

集合通信的实现方法大致如上所示。

结果

正确性验证

任何优化都必须以正确性的保证作为首要前提。在得到MPI优化的矩阵乘法结果后,利用前面的矩阵类Matrix实现的乘法得到另一个结果。然后,计算这两个结果的均方误差。注意到浮点数计算会存在一定的误差,只要误差是一个接近于零的值,则表明MPI的实现方法是正确的。

代码实现如下。

性能比较

测试样例采用规模为N×N​的方阵,512N2048​,进程数量从1​到16​递增。得到的结果如下。

点对点通信的实现方法。

 N=512N=1024N=2048
p=20.019s0.090s1.333s
p=40.012s0.086s1.146s
p=80.186s0.257s1.013s
p=160.252s1.035s1.421s

集合通信的实现方法。

 N=512N=1024N=2048
p=20.015s0.142s1.403s
p=40.011s0.145s1.330s
p=80.059s0.331s1.509s
p=160.248s0.765s1.833s

加速比如下所示。

3-1

从整体上来看,当并行的进程数量没有超过系统中的CPU的核数(测试机上的核数为6)时,MPI的实现均能够对原来的串行版本起到加速效果。同时,加速效果随着进程数量的增多而有所增强。但是,当并行的进程数量超过系统中的CPU的核数后,加速比就降至1以下,即矩阵乘法的计算反而变慢。这是因为当并行的进程数量不大于系统中的CPU的核数时,并行的进程可以在CPU的不同的核上运行,同时执行独立的计算任务,从而能够实现加速。但是,当并行的进程数量大于系统中的CPU的核数后,必然会存在若干个进程运行在同一核上。运行在同一核上的进程是并发执行的,这些进程之间是会相互影响的,使得进程在计算过程中不可避免地会由于操作系统的调度而产生阻塞和睡眠,导致了进程的计算时间变长。

在并行的进程数量没有超过系统中的CPU的核数的前提下,分析MPI实现的差别。

对于矩阵规模小于10243​下,集合通信的实现方法的加速效果比点对点的实现方法更佳;但当矩阵规模大于等于