矩阵乘法的MPI并行程序

深圳大学  并行计算  实验六

一、实验目的

1. 学会编写矩阵乘法的MPI并行程序；

2. 掌握6个基本的MPI语句；

3. 对并行程序进行性能分析。

二、实验环境

1. 硬件环境：6×32核CPU、400GB内存的分布内存并行计算平台；

2. 软件环境：Ubuntu Linux、gcc、g++（g++ -O3 -fopenmp -o a.out a.cpp）；

3. 远程登录：本地PowerShell中执行ssh bxjs@hpc.szu.edu.cn；

4. 传输文件：本地PowerShell中执行scp c:\a.cpp bxjs@hpc.szu.edu.cn:/home/bxjs/或ftp://hpc.szu.edu.cn。

三、实验内容

1. 用MPI编写两个n阶方阵a和b的并行相乘程序，结果存放在方阵c中。初始时，a和b都存储在进程0中，其他进程没有数据。结束时，c也存储在进程0中。为了验证结果正确性，将并行计算结果和串行计算结果相比较。

首先，我们考虑用为每一个进程程序，都构建三个长度为n*n的double类型的一维数组a、b、c，用来存储矩阵数据，用a[i*n+j]来表示矩阵第i行第j列的数据，并且考虑将并行计算机的结果存储在c数组所代码的矩阵中，并最终规整到进程编号为0的进程的c数组中。

然后，我考虑使用MPI库中，MPI_Init函数来初始化 MPI 的内部状态和资源，使得 MPI 可以管理并控制进程间的通信；MPI_Comm_rank函数来获取当前进程在通信器中的编号为rank，以此来区分不同进程的不同执行逻辑；MPI_Comm_size函数来获取通信器中包含的进程总数，方便对任务进行划分和规划。

接着，我们需要考虑对数组a、b进行随机的初始化，同时，要保证数据在不同进程间的一致性，所以先只对进程编号为0的进程的a、b进程随机初始化。对此，我采用了cstdlib库中rand()函数获取一个0到RAND_MAX之间的随机整数，再将随机数乘以1.0变成double类型，再除以RAND_MAX，最终获得一个0到1之间的随机的浮点数，然后加上一个数组元素的值，这样就可以获得两个随机的有序的double类型的一维数组。同时，需要对每一个进程程序的c数组进行初始化为0。

在完成初始化后，我考虑对矩阵乘法的计算以a矩阵的行数为计算总量，进行分组，同时，定义两个数组shuzu_a和shuzu_c来存储后续进行矩阵计算时，分割到的需要进行计算中使用的a矩阵的数据到shuzu_a中，并将计算结果存储到shuzu_c中。再用MPI_Wtime()函数，获取当前进程本地的壁钟时间，赋值给double类型的变量t0，作为并行域计算的起始时间。

    然后，我考虑到矩阵a和矩阵b的乘法中，在对矩阵a按行进行计算量的分割后，矩阵b的所有数据在后续计算中基本都会被不同程度的使用，所以需要考虑讲数组b广播给其他所有的进程。这里采用了使用MPI库中MPI_Bcast函数，在进程编号为0的进程中讲数组b进行广播。

在广播完b数组后，需要考虑讲对a矩阵按行分割后的计算中，所需要那几行的数据发送给对应的进程，所以这里采用了在进程编号为0的进程中使用MPI_Recv函数将a数组的部分数据发送给对应的进程，而在其他编号的进程中使用MPI_Send函数接收进程编号为0的进程发送过来的a数组的数据。

在完成并行计算前的数据交换后，通过一个三层for循环，对分割到的矩阵乘法计算的部分任务进行计算，并将结果存储到shuzu_c数组中。在完成并行计算后，仿照前面发送a数组部分数据的思路，将各个进程的并行计算的结果shuzu_c，发送给进程编号为0的进程。同时，进程编号为0的进程接收各进程发送过来的shuzu_c数据存储到对应的c矩阵所在行中。

接着，再次使用MPI_Wtime()函数，获取当前进程本地的壁钟时间，赋值给double类型的变量t1，作为并行域计算的结束时间。继而，通过t1-t0就可以获得该并行矩阵乘法的运行时间。然后，调用MPI_Finalize函数是清理MPI环境并结束MPI会话。

最后，由进程编号为0的进程进行一个三层for循环的串行矩阵a乘以矩阵b的计算，并将计算结果和并行计算c数组的结果进行比较，判断结果是否一致，来判断并行运行的正确性。具体最后实现的代码，在下列“四、代码描述”中可见。

2. 测试并行程序在不同进程数下的执行时间和加速比（与进程数=1时的执行时间相比）。其中，n固定为1000，进程数分别取1、2、4、8、16、32、64时，为减少误差，每项实验进行5次，取平均值作为实验结果。

如下图1所示，在当前目录下创建machinefile文件，并写下如下内容，用来作为后续程序运行命令中的-f参数所指定的文件，指示运行 MPI 任务的主机列表。

图1  machinefile文件的创建

然后，后续以“mpiexec -n 16 -f machinefile ./a.out”为例子，通过修改命令行中的-n的数值大小来控制程序运行的进程数分别为1、2、4、8、16、32、64，并进行5次取平均值作为该进程数下的运行结果。具体的实验结果，在下列“五、实验结果和分析”中可见。

四、代码描述

如下所示，为该实验的完整实现代码， 在main函数中，主要通过在计算前先提取发送数据和回收计算后的结果来实现的多个进程并行实现矩阵a和矩阵b的矩阵乘法的操作。

其中，开始通过MPI_Bcast函数广播b数组，假定采用的是二叉树广播，根进程向一半的进程发送数据，接着这些进程再向其他进程发送，时间复杂度为 。发送a数组的时间开销与a数组的大小直接相关，时间复杂度为 。

然后，根据分割的计算量进行并行地计算，时间复杂度为 。 最后，由进程编号为0的进程接收并行计算的结果，时间复杂度c数组的大小直接相关，时间复杂度为 ​​​​​​​ 。

    所以，总共完成矩阵a和矩阵b的并行矩阵乘法，需要的时间是。

#include "mpi.h"
#include <bits/stdc++.h>
using namespace std;
int main(int argc, char* argv[])
{
	//申请内存，获取进程编号、进程数量 
	const int n = 1000;
	double*a=new double[n*n];
	double*b=new double[n*n];
	double*c=new double[n*n];
	int rank,size;
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD,&size);
	//初始化 
	if(rank==0){
		srand(0);
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++){
			a[i*n+j]=rand()*1.0/RAND_MAX;
			b[i*n+j]=rand()*1.0/RAND_MAX;
			c[i*n+j]=0;
			}
	}else{
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++)
			c[i*n+j]=0;
	}
	int kuai_size=n/size,extra_size=n%size;
	MPI_Status status;
	double t0 = MPI_Wtime();
	//广播b数组 
	MPI_Bcast(b,n*n,MPI_DOUBLE,0,MPI_COMM_WORLD);
	//开局存 
	double* shuzu_a;
	double* shuzu_c;
	if (rank == 0) {
      	shuzu_a=new double[(kuai_size+extra_size)*n];
            shuzu_c=new double[(kuai_size+extra_size)*n];
	}else {
            shuzu_a=new double[kuai_size*n];
            shuzu_c=new double[kuai_size*n];
	}
	//发送部分的a 
	if(rank==0){
		for(int i=1;i<size;i++)
	MPI_Send(a+(kuai_size*i+extra_size)*n,kuai_size*n,MPI_DOUBLE,i,110,MPI_COMM_WORLD);
		for(int i=0;i<(kuai_size+extra_size)*n;i++)
			shuzu_a[i]=a[i]; 
	}else{
	MPI_Recv(shuzu_a,kuai_size*n,MPI_DOUBLE,0,110,MPI_COMM_WORLD,&status); 
	}
	//并行作业
	for(int i=0;i<(rank==0?kuai_size+extra_size:kuai_size);i++){
		for(int j=0;j<n;j++){
			shuzu_c[i*n+j]=0;
			for(int k=0;k<n;k++)
			shuzu_c[i*n+j]+=shuzu_a[i*n+k]*b[k*n+j];
		}
	}
	//发送结果shuzu_c
	if(rank==0){
		for(int i=0;i<(kuai_size+extra_size)*n;i++)
		c[i]=shuzu_c[i]; 
		for(int i=1;i<size;i++)
	MPI_Recv(c+(kuai_size*i+extra_size)*n,kuai_size*n,MPI_DOUBLE,i,130,MPI_COMM_WORLD,&status);
	}else{
	MPI_Send(shuzu_c,kuai_size*n,MPI_DOUBLE,0,130,MPI_COMM_WORLD);
	}
	double t1 = MPI_Wtime();
	MPI_Finalize();
	//串行验证 
	if (rank==0){
		cout<<"time is "<<t1-t0<<" seconds"<<endl;
		int flag=1;
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++){
			double s=0;
			for(int k=0;k<n;k++) s+=a[i*n+k]*b[k*n+j];
				if(s!=c[i*n+j]){
	cout<<"s="<<s<<"c["<<i*n+j<<"]="<<c[i*n+j]<<endl;
					flag=0;
				} 
			}
		if(flag) cout<<"Results are equal!"<<endl;
		else cout<<"Results are not equal!"<<endl;
	}
	delete[] shuzu_a;
	delete[] shuzu_c;
	delete[] a;
	delete[] b;
	delete[] c;
	return 0;
}

五、实验结果和分析

通过在文件目录下输入“ftp://hpc.szu.edu.cn”，使用ftp协议，将本地的a.ccp的代码文件传输到“/2021150233”的目录下。然后，打开shell窗口，使用“ssh bxjs@hpc.szu.edu.cn”连接A0服务器。

然后，通过cd命令打开“/202115023”文件夹，接着，通过“mpic++ -O3 -o a.out a.cpp”命令行将a.ccp文件编译为a.out可执行文件。然后通过“mpiexec -n 1 -f machinefile ./a.out”命令行运行a.out文件，并通过调节命令行中的-n参数来调整程序运行的进程数的设置，最终进程数分别取1、2、4、8、16、32、64时，执行的结果如下图2所示，并将执行的结果记录进表1中，同时计算出加速比（与线程数=1时的执行时间相比）。

 

图2  a.ou文件的执行结果

表1 并行程序在不同线程数下的执行时间（秒）和加速比（n=1000）

线程数 执行时间/s	1	2	4	8	16	32	64
第1次	5.61343	2.57548	1.57893	1.08825	0.77584	0.68420	0.64622
第2次	5.17211	2.71941	1.54579	1.05442	0.88674	0.72676	0.63017
第3次	5.02107	2.86850	1.52255	0.97224	0.80058	0.72967	0.68519
第4次	5.21602	2.80546	1.50722	1.04337	0.79589	0.67116	0.64936
第5次	4.82169	2.72237	1.51086	0.98791	0.73610	0.69187	0.64596
平均值	5.05772	2.73824	1.53307	1.02924	0.79903	0.70073	0.65138
加速比	1.00000	1.84707	3.29908	4.914043	6.32983	7.21777	7.76462

根据上述表1的平均值和加速比数据，绘制出下面的折线图：

图3  平均值和加速比的折线图

根据对上述，表1和图3的分析可知，随着线程数的增加，n=1000的矩阵乘法运算，所需要的时间在逐渐下降，加速比逐渐上升。即，通过增加并行进程数，分割任务进行并行的矩阵乘法可以实现对矩阵乘法运算的优化。

同时，观测到加速比不是简单随线程数的增加而线性增加的，两者不是简单的正比关系，加速比随着进程数的增加而增加，但是增长越来越缓慢了。这与前面代码描述部分中分析的“总共完成矩阵a和矩阵b的并行矩阵乘法，需要的时间是 ​​​​​​​”有关系，执行时间与线程数并不是简单的反比关系。

除此之外，在对矩阵乘法的计算量进行并行分割时，很多时候可能是并不能实现较为均匀的分割的，所以这也在很大程序上影响着程序的运行时间。并且随着进程数的增加，在采用MPI实现进程间的数据交互时，数据交互的代价可能也随之上升，导致加速效果也在一定程度上被削弱，最终导致两者呈现出这种关系。

六、实验结论

在准备实验的初期，我对如何使用MPI来进行并行编程的使用方法存在或多或少的疑惑。虽然，我们之前已经用opnemp对矩阵乘法进行并行化的优化，但是鉴于初次接触MPI，我甚至不知道从何下手，因为我缺乏足够的MPI并行编程经验和对MPI相关运行机理等的理解。

然而，我并没有放弃，而是积极地寻找解决办法。我开始查阅相关资料，阅读MPI并行编程相关的文献和教程，深入了解了MPI基础的运行机理和编程方法。通过学习和实践，我逐渐掌握了如何用MPI将矩阵乘法的计算任务分解成多个并行子任务，并利用多个服务器进程资源进行加速。

在上手写代码实现的过程中，我遇到了许多挑战和问题。例如，如何合理地传输矩阵乘法所需的数据，如何高效对计算任务数据进行划分等等。但是，每当我遇到困难时，我都会认真思考，并寻找解决办法。通过不断地尝试和调试，我逐步解决了这些问题，最终成功地用MPI实现了矩阵乘法的并行程序。

在本次实验中，我们实现了使用 MPI (Message Passing Interface) 并行编程模型对两个 n 阶方阵的矩阵乘法进行了并行计算，并进行了性能分析。利用了 MPI 的 6 个基本函数（MPI_Init、MPI_Comm_rank、MPI_Comm_size、MPI_Send、MPI_Recv、MPI_Finalize）和MPI_Bcast函数进行进程间通信和数据传输，完成了并行计算任务。

实验数据表明，总体上执行时间与线程数呈现递减趋势，但由于通信开销的存在，执行时间和加速比并不是简单的反比关系。随着进程数增加，程序运行的总时间减少，但并行效率提高幅度逐渐减小，甚至在高进程数时可能出现不理想的加速效果。

总的来说，本实验通过使用MPI设计并实现矩阵乘法的并行算法，不仅掌握了MPI并行编程的基本方法和技术，同时通过性能分析对并行程序进行评估，深入理解了并行计算的原理和优化策略，为进一步深入学习并行计算奠定了基础。