矩阵乘法的MPI并行程序

最新推荐文章于 2024-09-26 09:13:29 发布

绝不收费

最新推荐文章于 2024-09-26 09:13:29 发布

阅读量1k

点赞数 12

分类专栏：并行计算实验文章标签： c++

本文链接：https://blog.csdn.net/weixin_60705585/article/details/140508909

版权

并行计算实验专栏收录该内容

8 篇文章 0 订阅

订阅专栏

深圳大学并行计算实验七

一、实验目的

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

2. 掌握MPI_Sendrecv、MPI_Isend、MPI_Reduce等MPI语句；

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

二、实验环境

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

2. 软件环境：Ubuntu Linux、mpicc、mpic++（mpic++ -O3 -o a.out a.cpp）、mpiexec（mpiexec -n 36 -f machinefile /home/bxjs/a.out）

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的Cannon矩阵乘法程序，结果存放在方阵C中。为了简单起见，假设进程数p为完全平方数，且n可以被整除。a、b、c都采用块棋盘划分存储，每个子块的大小都为。A和B的初始值为：Ai,j=Bi,j=i×n+j，0≤i,j<n。在计算完C=AB后，求C中所有元素的和sum，即sum= ，然后进程0在终端输出sum值。

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

然后，通过获取的处理器总数信息，来对数据块进行针对性的进行划分。先使用math.h库中的sqrt函数，将处理器隐式的划分为q*q的排列，然后，通过m=n/q的计算，将n*n的矩阵划分为了m*m的矩阵，每个处理器就隐式的拥有着A矩阵、B矩阵和C矩阵划分出来的m*m的子块。同时，根据前面获得的进程在通信器中的编号rank，根据rank/q和rank%q可以获知改子块在对应的A矩阵、B矩阵和C矩阵中的位置为 [rank/q，rank%q] 。

针对这些矩阵的子块数据，需要将其具体化，要为其提供数据的表示和存储。所以，考虑用为每一个进程程序，都构建三个长度为m*m的int类型的一维数组a、b、c，用来存储A、B、C矩阵在进程中的对应的块数据，同时，可以用a[i*m+j]来表示子块中第i行第j列的数据，也表示着矩阵中第rank/q*m+i行第rank%q*m+j的数据，并且考虑将并行计算过程中的结果存储在c数组。除此之外，为了满足后续进行进程间数据通信时的数据交换的需要，需要再构建一个长度为m*m的int类型的一维数组buf。

接着，我们需要考虑对数组a、b进行初始化，这里采用的初始化方式为将a[i*m+j]和b[i*m+j]都初始化为(rank/q*m+i)*n+rank%q*m+j对应的值。同时，需要对每一个进程程序的c数组进行初始化为0，方便后续存储矩阵乘法结果的需要。

在完成初始化后，用MPI_Wtime()函数，获取当前进程本地的壁钟时间，赋值给double类型的变量t0，作为并行域计算的起始时间。

然后，根据Cannon矩阵乘法的需要，需要先对一些进程的数据块进行循环的移动，来满足后续并行矩阵乘法中连续的数据交换的需要。“将所有进程的Ai,j(0≤i, j≤-1)分块向左循环移动i步(按行移位)，将所有块Bi,j(0≤i,j≤-1)向上循环移动j步(按列移位)”，具体在代码实现上的体现就是，将a数组数据发送给yi*q+((er-yi+q)%q)进程（这里yi=rank/q, er=rank%q），然后从yi*q+((er+yi)%q)进程接收传输过来的数据来替换a数组；将b数组数据发送给((yi-er+q)%q)*q+er进程，然后从((yi+er)%q)*q+er进程接收传输过来的数据来替换b数组的内容。

同时，这里在涉及数据的发送和接收时，由于数组a、b的大小为m*m，是比较大的，缓冲区不足以放下，所以在发送时需要采用非阻塞式的MPI_Isend函数，接收时可以采用阻塞式的MPI_Recv函数。并且，为了防止数组a或b还没发送完，就已经获取到了来替换的数据了，而发生数据的数据覆盖和出错问题。在接收时，需采用buf数组接收，并且需要使用MPI_Wait函数来保证非阻塞式的MPI_Isend函数结束了，才能将buf数组的数据拷贝到对应的a数组和b数组中。

在完成了并行矩阵乘法前的数据移位后，就可以考虑根据Cannon矩阵乘法进行并行计算了。这里我采用了do-while循环来实现，在循环里面我用一个三层for循环来实现子块a和子块b的矩阵乘法，并将结果放置到子块c中。然后，进行计数判断，最多进行q次这样的子块矩阵乘法的计算。未跳出循环，则“将A的每个块向左循环移动一步，将B的每个块向上循环移动一步”，具体在代码实现上的体现就是，将a数组数据发送给yi*q+((er-1+q)%q)进程，然后从yi*q+((er+1)%q)进程接收传输过来的数据来替换a数组；将b数组数据发送给((yi-1+q)%q)*q+er进程，然后从((yi+1)%q)*q+er进程接收传输过来的数据来替换b数组的内容，与前面的数据交换相似。

在完成q次子块矩阵乘法后，将c数组的值累加到sum变量中。然后调用MPI_Reduce函数“MPI_Reduce(&sum,&rec_sum,1,MPI_LONG_LONG,MPI_SUM,0,MPI_COMM_WORLD);”，将p个进程的sum遍历的值累加到0号进程的rec_sum变量中，再将rec_sum变量赋值给sum。

接着，再次使用MPI_Wtime()函数，获取当前进程本地的壁钟时间，赋值给double类型的变量t1，作为并行域计算的结束时间。继而，通过t1-t0就可以获得该并行矩阵乘法的运行时间。然后，对sum的值和程序运行的时间信息进行输出。

最后，使用delete来回收数组的内存，调用MPI_Finalize函数是清理MPI环境并结束MPI会话。而针对于sum的输出的值，可以与先前计算好的该确定的A、B矩阵的乘法结果进行比较，判断结果是否一致，来判断并行运行的正确性。具体最后实现的代码，在下列“四、代码描述”中可见。

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

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

图1 machinefile文件的创建

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

四、代码描述

如下所示，为该实验的完整实现代码，在main函数中，先申请数组空间并初始化，然后进行并行计算前面的数据交换，实现子块数据的移位，然后通过do-while循环进行q次（其中进程数=q*q）并行的子块矩阵乘法和相应q-1次的数据交换，从而完成矩阵的乘法运算。

其中，数据交换，因为发送和接收缓存空间的大小是小于发送数据的块的大小的，所以通过非阻塞式MIP_Isend发送数据，和阻塞式MIP_Recv接收数据，并且用MIP_Wait保证数据的覆盖是在发送完成之后。同时，程序中的三层for循环，计算进程中子块a和子块b的矩阵乘法中，我针对性的调整了for循环的顺序，来提高for循环过程中的cache的命中率，来提升程序的性能。在进行数组数据的覆盖式，我采用的memcpy来提高数组数据的拷贝效率。

在完成举证乘法后，各个进程对各自的c数组进行求和，然后通过MPI_Reduce函数将结果，累加到0号进程中，由0号进程对结果进行输出。

在这个程序运行过程中，主要的时间开销来源于前面的矩阵乘法的计算。在这个矩阵乘法过程中，每一个进程都平均需要完成q次的a数组子块数据交换和覆盖、q次的b数组子块数据交换和覆盖,以及q次子块矩阵的乘法。

#include "mpi.h"
#include <iostream>
#include <math.h>
#include <string.h>
using namespace std;
int main(int argc, char** argv)
{
	int n=1680; //矩阵大小为n*n
	long long sum=0;
	int rank,p;
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD,&p);
	int q=int(sqrt(double(p))); //处理器数为q*q 
	int m=n/q; //数据块大小为m*m 
	int* a=new int[m*m]; //用一维数组表示子块，a[i*m+j]为子块第i行第j列元素
	int* b=new int[m*m];
	int* c=new int[m*m];
	int* buf=new int[m*m];
	for (int i=0;i<m;i++)
		for (int j=0;j<m;j++){
		a[i*m+j]=b[i*m+j]=(rank/q*m+i)*n+rank%q*m+j;
		c[i*m+j]=0;
		}
	double t0=MPI_Wtime();
	//在以下部分添加并行代码，main函数其他地方保持不变，不可再新建数组。
	int yi=rank/q,er=rank%q,temp=0; //所在位置为 [yi,er] 
	MPI_Status status;
	MPI_Request send_request;
	//提前偏移 
	if(yi>0){
		MPI_Isend(a,m*m,MPI_INT,yi*q+((er-yi+q)%q),231,MPI_COMM_WORLD, &send_request);
		MPI_Recv(buf,m*m,MPI_INT,yi*q+((er+yi)%q),231,MPI_COMM_WORLD, &status);
		MPI_Wait(&send_request, &status);
		memcpy(a,buf,m*m*sizeof(int));
	}
	if(er>0){
		MPI_Isend(b,m*m,MPI_INT,((yi-er+q)%q)*q+er,123,MPI_COMM_WORLD, &send_request);
		MPI_Recv(buf,m*m,MPI_INT,((yi+er)%q)*q+er,123,MPI_COMM_WORLD, &status);
		MPI_Wait(&send_request, &status);
		memcpy(b,buf,m*m*sizeof(int));	
	}
	//并行计算机分块矩阵 和 A、B的循环轮换 
	do{
		//分块矩阵乘法 
		for(int i=0;i<m;i++)
		for(int k=0;k<m;k++)
			for(int j=0;j<m;j++)
			c[i*m+j]+=a[i*m+k]*b[k*m+j];
		if(++temp==q) break;
		//A矩阵的循环轮换 
	MPI_Isend(a,m*m,MPI_INT,yi*q+((er-1+q)%q),321+temp, MPI_COMM_WORLD,&send_request);
		MPI_Recv(buf,m*m,MPI_INT,yi*q+((er+1)%q), 321+temp, MPI_COMM_WORLD,&status);
		MPI_Wait(&send_request,&status);
		memcpy(a,buf,m*m*sizeof(int));
		//B矩阵的循环轮换
		MPI_Isend(b,m*m,MPI_INT,((yi-1+q)%q)*q+er,132+temp, MPI_COMM_WORLD,&send_request);
		MPI_Recv(buf,m*m,MPI_INT,((yi+1)%q)*q+er,132+temp, MPI_COMM_WORLD,&status);
		MPI_Wait(&send_request,&status);
		memcpy(b,buf,m*m*sizeof(int));	
	}while(1);
	//求据矩阵C的和 
	for(int i=0;i<m*m;i++) sum+=c[i];
	long long rec_sum;
	MPI_Reduce(&sum,&rec_sum,1,MPI_LONG_LONG,MPI_SUM,0,MPI_COMM_WORLD); 
	sum=rec_sum;
	double t1 = MPI_Wtime();
	if (rank == 0)
		cout<<"sum is "<<sum<<endl<<"time is "<<t1-t0<<" seconds"<<endl;
	delete[]a;
	delete[]b;
	delete[]c;
	delete[]buf;
	MPI_Finalize();
	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、4、9、16、25、36、49、64时，执行的结果如下图2所示，并将执行的结果记录进表1中，同时计算出加速比（与线程数=1时的执行时间相比）。

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

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

进程数执行时间/s	1	4	9	16	25	36	49	64
第1次	2.89266	1.30866	0.58598	0.44086	0.38822	0.28133	0.30589	0.35962
第2次	2.63143	1.33882	0.58782	0.43107	0.38257	0.24617	0.28209	0.38724
第3次	2.70635	1.24244	0.57660	0.45349	0.36055	0.29789	0.28745	0.37237
第4次	2.20432	1.31150	0.53749	0.43695	0.39060	0.26278	0.30967	0.37083
第5次	2.54433	1.07306	0.62559	0.44890	0.40158	0.24978	0.29668	0.36671
平均值	2.52161	1.24146	0.58188	0.44260	0.38383	0.26416	0.29397	0.37429
加速比	1.00000	2.03117	4.33359	5.69723	6.56969	9.54595	8.57771	6.73709

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

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

根据对上述的表1和图3进行分析可知，随着线程数的增加，n=1680的矩阵乘法运算，所需要的时间在先逐渐下降然后逐渐上升，加速比逐渐上升然后逐渐下降。即，通过增加并行进程数，分割任务进行并行的矩阵乘法可以实现对矩阵乘法运算的优化，但是随着进程数的增加，到达一定的值后，优化的效果逐渐减弱。

这主要是与前面代码描述部分中分析的“矩阵乘法过程中，每一个进程都平均需要完成q次的a数组子块数据交换和覆盖、q次的b数组子块数据交换和覆盖,以及q次子块矩阵的乘法”有关系，随着进程数的增加，子块大小在逐渐减少，每个进行子块矩阵的乘法的计算的时间开销在减小，通信的数据量大小是保持不变的，但是通信次数的增加和通信结点是在增加的。所以，一开始，由于通信结点的数量较少，通信的代价并不会非常大，并行子块矩阵的乘法的计算的时间带来的优化较大。但是随着进程数的增加，通信次数的增加，通信结点增加导致通信网络的复杂度的增加，导致通信的开销越来越大，逐渐覆盖了，子块矩阵的乘法的计算的时间带来的优化，导致了运行时间的回升和加速比的回降，最终导致呈现出这种关系。