高斯-约旦消去法的MPI并行程序

深圳大学  并行计算  实验八

一、实验目的

1. 学会编写高斯-约旦消去法的MPI并行程序；

2. 掌握MPI_Allgather、MPI_Bcast、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编写选主元的高斯-约旦消去法的MPI并行程序，求解n元线性方程组Ax=b。为了验证结果正确性，将求解得到的x带入原方程组，求出最大误差，即向量Ax-b中绝对值最大的分量。

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

然后，通过获取的处理器总数和当前进程在通信器中的编号的信息，来对数据块进行针对性的进行划分，让每个进程负责m行线性方程组的数据。rank小于n%p的进程，负责m=n/p+1行线性方程组的数据的存储和；rank大于n%p的进程，负责m=n/p行线性方程组的数据的存储和。同时，不用刻意去区别不同进程负责的线性方程组的相应的顺序。然后，定义double类型的二维数组a[m][n+1]和z[m][n+1]来存储各自进程负责的线性方程组的数据，其中，z数组主要是为了备份a数组的数据，在后续计算误差时使用。而b参数就存储在a[i][n]和z[i][n]中。

接着，定义double类型的一维数组max_coefficient[p]，主要用于后续主元选取过程中的数据交换的存储。定义double类型的一维数组x[n]，用于存储线性方程组的解向量。定义int类型的一维数组pivot[m]用于指示当前进程负责的m行线性方程组被选择的作为主元的列位置，方便后续计算解向量x。定义int类型的一维数组which[n]来指示每一列中被选作主元的所在行是由哪个进程负责，方便最后对解向量x的广播共享。

紧接着，对pivot数组进行初始化为-1，表示负责的行还未选取出主元。然后，通过rand()*1.0/RAND_MAX生成小于1的随机浮点数，赋值给a数组和z数组，作为当前进程负责的几行线性方程组的相关参数。所有的进程的数据就构成了一个完整的随机线性方程组的所有系数。

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

然后，根据高斯-约旦消去法，构造一个循环n次的for循环，来实现将a数组存储的线性方程组变为每行只有一个主元为非0参数的线性方程组。在第k次循环中，每个进程首先各自通过一个m次，寻找自己负责的行数据在对应的想要寻找主元的k列上的最大值local_max和其所在的行local_row（当然这里已经出现主元的行不进行查找）。然后，通过MPI_Allgather函数实现将各个进程的寻找到的local_max进行数据交互，并各个进程找的local_max，都存储在每个进程的max_coefficient数组中，然后通过distance函数，找出最大local_max是在进程编号为max_index的进程，并将switch[k]赋值为max_index。

此时，进程编号为max_index的进程就需要根据local_row将pivot[local_row]设置为k，来表示当前local_row行已经选取出第k列的数据作为主元了，并对应的local_row行数据拷贝出来存储picked_row数组。然后，使用MPI_Bcast函数对picked_row数组进行广播，让每个进程都获取到当前选取出来主元行的各个位置上的参数。接着，每个进程根据picked_row数组的数据，通过一个m*n的双层for循环将自己负责的线性方程组的行（当然，除了先前选出的主元行picked_row）中的第k个参数进行消元。

就这样循环上述的操作k次，就可以完成a数组存储的线性方程组变为每行只有一个主元为非0参数的线性方程组。然后，每个进程分别通过一个m次for循环将负责的行的选出的主元x[pivot[i]]变量进行求解，赋值为a[i][n]/a[i][pivot[i]]。最后，通过一个n次for循环，根据先前的switch数组的记录，通过使用MPI_Bcast函数，由switch[i]进程将解向量中的x[i]广播给所有进程，最后让每个进程都拥有一个最后完整的解向量x。

接着，再次使用MPI_Wtime()函数，获取当前进程本地的壁钟时间，赋值给double类型的变量t1，作为并行域计算的结束时间。继而，通过t1-t0就可以获得该并行高斯-约旦消去法程序的运行时间。

然后，还有计算线性方程组的最大误差，通过一个n*m的双层for循环将每个进程负责的线性方程组的所有行的最大误差local_dif计算出来，然后通过MPI_Reduce函数，从所有进程求解出来最大误差local_dif中，选取出一个最大误差max_dif，并由进行编号为0的进程输出最大误差的信息和进程运行时间的信息。

最后，使用delete来回收数组的内存，调用MPI_Finalize函数是清理MPI环境并结束MPI会话。而针对于max_dif的输出的值，可以通过判断max_dif为一个较小的浮点数值，来判断并行运行的正确性。具体最后实现的代码，在下列“四、代码描述”中可见。

2. 测试并行程序在不同进程数下的执行时间和加速比（与进程数=1时的执行时间相比）。其中，n固定为1000，进程数分别取1、2、4、8、16、32、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函数中，先获取进程通信的基本信息来确定负责的线程方程组的行数，然后申请相应数组空间作为线性方程组存储和并行计算的需要，并进行随机初始化作为线性方程组对应行的数组。

接着，通过一个n次for循环实现，选取主元，并进行线性方程组的消元的操作，实现将线性方程组矩阵转变为一个每行只有一个主元为非0参数的线性方程组。其中，主要通过“MPI_Allgather(&local_max, 1, MPI_DOUBLE, max_coefficient, 1, MPI_DOUBLE, MPI_COMM_WORLD);”进行进程间通信，从而实现主元的确定；主要通过“MPI_Bcast(& picked_row[k], n+1-k, MPI_DOUBLE, max_index,MPI_COMM_WORLD);”对主元所在行方程系数进行广播。

然后，每个进程计算对应行的主元的值，即为解向量x中的一个维度上的值，最后通过MPI_Bcast函数对各个进程中的部分的解向量x进行广播，实现每个进程上都拥有完整的解向量x。

在线性方程组的求解后，各个进程将解进行带入，求解每个方程的误差，并选取出最大的误差local_dif，然后通过MPI_Reduce函数获取进程间的最大的local_dif到0号进程中，由0号进程对结果进行输出。

在这个程序运行过程中主要时间消耗来源于，每一个进程都平均需要完成O(n×m×n)次的消元运算，O(n)次的对O(n)级别的数据的广播。

#include <iostream>
#include <cmath>
#include <algorithm>
#include <string.h>
#include "mpi.h"
using namespace std;
int main(int argc, char* argv[])
{
	int rank, p;
	const int n = 1000;
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	MPI_Comm_size(MPI_COMM_WORLD, &p);
	int m = rank < n % p ? n / p + 1 : n / p; //该进程处理的矩阵行数
	double(*a)[n + 1] = new double[m][n + 1]; //b[i]存储在a[i][n]中
	double(*z)[n + 1] = new double[m][n + 1]; //a的备份，用于计算最大误差
	double* max_coefficient = new double[p]; // 存储每个进程找到的最大主元值
	int* pivot = new int[m]; // 存储每行的主元列
	int which[n];// 记录主元的行号
	double x[n]; // 解向量
	double picked_row[n + 1]; // 存储被选为主元的行的n+1个数据 
	for (int i = 0; i < m; i++)  pivot[i] = -1;
	srand(rank*100);
	for(int i=0;i<m;i++)// 随机初始化矩阵a和备份矩阵z
		for(int j=0;j<n+1;j++)
			a[i][j]=z[i][j]=rand()*1.0/RAND_MAX;
	// 获取开始时间
	double t0 = MPI_Wtime();
	// 开始高斯-约旦消去法
	for (int k = 0; k < n; k++) {
		// 确定主元 
		double local_max = -RAND_MAX-1;
		int local_row = -1;
		for (int i = 0; i < m; i++) {
			if (pivot[i] == -1 && a[i][k] > local_max) {
			local_max = a[i][k];
			local_row = i;
			}
		}
		MPI_Allgather(&local_max, 1, MPI_DOUBLE, max_coefficient, 1,
					 MPI_DOUBLE, MPI_COMM_WORLD);
		int max_index = distance(max_coefficient, max_element(
				max_coefficient, max_coefficient + p));
		// 同步主元行数据 
		if(max_index==rank){
		memcpy(picked_row, a[local_row], (n+1) * sizeof(double));
			pivot[local_row] = k;
		}
		which[k] = max_index;
		MPI_Bcast(&picked_row[k], n+1-k, MPI_DOUBLE, max_index,
							 MPI_COMM_WORLD);
		// 用主元行数据消元 
		for (int i = 0; i < m; i++) {
			if (pivot[i] !=k) {
			double factor = a[i][k] / picked_row[k];
			a[i][k] = 0;
			for (int j = k+1; j < n + 1; j++) {
				a[i][j] -= factor * picked_row[j];
			}
			}
		}
	}
	// 根据标志计算负责的行的主元值 
	for(int i=0;i<m;i++)
		x[pivot[i]]=a[i][n]/a[i][pivot[i]];
	// 广播各自的结果 
	for (int i = n - 1; i >= 0; i--)
		MPI_Bcast(&x[i], 1, MPI_DOUBLE, which[i], MPI_COMM_WORLD);
	// 获取结束时间
	double t1 = MPI_Wtime();
	// 计算最大误差
	double max_dif;
	double local_dif = 0;
	for (int i = 0; i < m; i++){
		double dif = z[i][n];
		for (int j = 0; j < n; j++)
		dif -= z[i][j] * x[j];
		local_dif = max(local_dif, abs(dif));
	}
	MPI_Reduce(&local_dif, &max_dif, 1, MPI_DOUBLE, MPI_MAX, 0,
							 MPI_COMM_WORLD);
	// 输出结果
	if (rank == 0){
		cout << "max difference is " << max_dif << endl;
		cout << "time is " << t1 - t0 << " seconds" << endl;
	}
	delete[] a;
	delete[] z;
	delete[] pivot;
	delete[] max_coefficient;
	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=1000）

进程数 执行时间/s	1	2	4	8	16	32	64
第1次	0.41976	0.60901	0.79552	1.06614	1.28682	2.08940	3.00667
第2次	0.41622	0.62445	0.81880	1.09330	1.22624	1.93889	2.97225
第3次	0.41905	0.59711	0.81761	1.09367	1.21281	1.94413	2.83950
第4次	0.41613	0.60464	0.78574	1.06966	1.26164	2.03220	2.96655
第5次	0.44109	0.62025	0.79317	1.06371	1.22930	1.98670	3.11017
平均值	0.42245	0.61109	0.80217	1.07730	1.24336	1.99826	2.97903
加速比	1.00000	0.69131	0.52663	0.39214	0.33977	0.21141	0.14181

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

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

根据上述图2的实验的运行结果，可以发现，由于程序在生成线性方程组的系数时，随机种子采用“srand(rank*100);”进行设置，所以导致程序的在多次进程数相同的运行中，其实随机种子在多次运行中并没有改变，导致在进程数相同的5次运行中max_dif是保持不变的，所以该现象是合理且正常的。

然后，对上述的表1和图3进行分析可知，随着进程数的增加，n=1000的线性方程组的求解，所需要的时间在逐渐上升，加速比在逐渐下降。即，在MPI通信模型下，通过增加并行进程数，分割任务进行并行的线性方程组的求解，并没有对线性方程组的求解起到优化，反倒是随着进程数的增加，运行时间越来越长。

这主要是与前面代码描述部分中分析的“在这个程序运行过程中主要时间消耗来源于，每一个进程都平均需要完成O(n×m×n)次的消元运算，O(n)次的对O(n)级别的数据的广播。”有关系。虽然随着进程数的增加，每个进程负责的线性方程组的行数在减少，O(n×m×n)次的消元运算在逐渐减少。但是进程数的增加也导致了进程间的通信越来越复杂，通信的开销越来越高，O(n)次的对O(n)级别的数据的广播的所需时间越来越高，远不是并行减少消元运算的增益所能抵消的。最终导致呈现出这种‘进程数越多，运行时间越长’的关系。

六、实验结论

在准备实验的初期，鉴于上课的讲解和课后的研究，我自觉已经对选主元的高斯-约旦消去法的算法了然于胸了。但是，我在一个下午的实验过程中磕磕碰碰，甚至对于我已经掌握的知识产生了怀疑，怀疑自己的掌握是否不彻底。

在调试了几个小时后，我发现自己的程序仍然跑不出运行max_dif很小的结果，我一步步的调试起来。我甚至开始怀疑起我之前所写的所有代码的正确性了，我平复了心情后，开始逐渐着手解决问题。我首先通过编写编写串行的选主元的高斯-约旦消去法的算法，然后通过一步一步添加并行代码，将最后的串行的选主元的高斯-约旦消去法程序修改为并行的串行的选主元的高斯-约旦消去法程序。

但是，但我不断调试过程中，我发现我可以保持max_dif为很小的值来保证程序逻辑的正确性，但是运行时间却和实验提供的参考的1.2s相差甚远。我发现自己在1000次的对picked_row数组进行广播所需的时间已经达到了10s了，远高于参考的1.2s了。然后，又陷入了自我怀疑中，但是通过多方调试，我可以保证确实是这个通信时间，并不是我代码问题。

然而，在多次调试的过程中，最终，发现服务器A1、A2崩溃了，无法运行程序了。然后，在群里信息的沟通下，我通过在machinefile中删除A1和A2后，成功跑到了1.2s。

除此之外，我遇到了许多挑战和问题。例如，运行时间过长的问题等等。但是，每当我遇到困难时，我都会认真思考，并寻找解决办法。通过不断地尝试和调试，我逐步解决了这些问题，最终成功地实现选主元的高斯-约旦消去法的并行程序。

本次实验使用 MPI并行编程模型实现线性方程组进行并行的选主元的高斯-约旦消去法的算法，并进行了性能分析，我不仅掌握了MPI并行编程的基本方法和技术，同时通过性能分析对并行程序进行评估，深入理解了选主元的高斯-约旦消去法的算法的原理和优化策略，为进一步深入学习并行计算奠定了基础。