MPI并行计算学习笔记6——行主元高斯消去法

最新推荐文章于 2024-04-01 01:26:26 发布
原创 最新推荐文章于 2024-04-01 01:26:26 发布 · 3.6k 阅读 · 17 ·
本内容遵循CC 4.0 BY-SA版权协议
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

1.运行环境：VS2017 + MPI

2.参考资料：MPI与OpenMP并行程序设计 （Michael J.Quinn 著； 程文光 武永卫译 ） 第十二章线性方程组的求解（P243）

3.算法说明：直接求解线性方程组，为保证数值稳定性，大部分情况下采用列主元高斯消元法，但理论上列主元高斯消元法的并行性没有行主元高斯消元法的可扩展性好（书上有详细介绍）；

         描述： 将矩阵进行 行交叉分解，即第0行分给0号进程，第1行分给1号进程,.....,第p-1行分给p-1号进程 ，第p行分给0号进程，依此循环分配； 求得第k行的主元（枢纽），并将该行所有元素及枢纽位置发送给每个进程，每个进程按照该行元素及枢纽位置将各进程中枢纽列的k + 1行至dim - 1行全部消为0；其中k从0，1，2，....，到dim - 2

4.MPI通信实现方法：

    方法一：进行第k次迭代时，将所得的枢纽及枢纽所在行的元素利用MPI群集通信函数广播至每个进程；

    方法二：采用非阻塞通信，各进程组成环，0号进程将消息发送给1号，1号发送给2号，......，p - 1号发送给0号；这一过程中使用                    MPI_Isend非阻塞发送，可实现异步计算；

    注：书上说到方法二是可扩展性相当好的方法，但在PC机上实现的效果非常差；个人猜测与硬件有关，会在以后的文章中探讨；也欢迎有经验的前辈来交流！！！

5.源代码（头文件Matrix.h与第三篇文章中的类似）

#include<iostream>
#include<cmath>
#include<ctime>
#include"Matrix.h"

#include"mpi.h"

using namespace std;

int myid, numprocs, masterNode;//进程标识号与进程总数，设为全局变量可在任意个函数中访问
int partRow,up,down;//各进程分配到的行数，相邻的上下两个进程
int dim = 1500;//问题规模
MPI_Comm partComm;//不能整除的部分在该通信域做一次scatter/gather即可
double start, finish;

void seqGaussSolver(Matrix A, Matrix b, Matrix &x)//串行 行主元高斯消去法
{
	int N = A.row, picked = 0;
	int *loc = new int[N] {0};
	for (int i = 0; i < N; ++i)
		loc[i] = i;
	for (int k = 0; k < N; ++k)
	{
		//找行主元
		double max = 0;
		for (int j = k; j < N; ++j)
		{
			if (abs(A(k, loc[j])) > max)
			{
				max = abs(A(k, loc[j]));
				picked = j;
			}
		}
		//picked为当前行k的主元，但当前行记录的主元为loc[k]
		int tmp = loc[k];
		loc[k] = loc[picked];
		loc[picked] = tmp;

		//以枢纽为基准开始消元
		for (int i = k + 1; i < N; ++i)
		{
			double t = A(i, loc[k]) / A(k, loc[k]);
			for (int j = k; j < N; ++j)
				A(i, loc[j]) -= t * A(k, loc[j]);
			b(i, 0) -= t * b(k, 0);
		}
	}
	for (int i = N - 1; i >= 0; --i)
	{
		x(loc[i], 0) = b(i, 0) / A(i, loc[i]);
		for (int j = 0; j < i; ++j)
			b(j, 0) -= x(loc[i], 0)*A(j, loc[i]);
	}
	//cout << x << endl;
	delete[]loc;

}
void parallelInit()
{
	partRow = dim / numprocs;
	int remainder = dim % numprocs;
	if (myid < remainder)++partRow;
	
        //进程成环
	up = myid - 1, down = myid + 1;
	if (myid == 0)up = numprocs - 1;
	if (myid == numprocs - 1)down = 0;

        //创建新的通信域，方便不能整除时的播撒数据
	MPI_Group worldGroup;
	MPI_Comm_group(MPI_COMM_WORLD, &worldGroup);
	MPI_Group partGroup;

	int *groupRank = new int[remainder];
	for (int i = 0; i < remainder; ++i)
		groupRank[i] = i;
	
	MPI_Group_incl(worldGroup, remainder, groupRank, &partGroup);
	MPI_Comm_create(MPI_COMM_WORLD, partGroup, &partComm);
	delete[]groupRank;
}
void distributeTask(const Matrix &lAb,Matrix &Ab)
{
	
	int counts = dim / numprocs;
	for (int i = 0; i < counts; ++i)
		MPI_Scatter(&lAb(i*numprocs, 0), dim + 1, MPI_DOUBLE, &Ab(i, 0), dim + 1, MPI_DOUBLE, masterNode, MPI_COMM_WORLD);
	if (counts != partRow)
		MPI_Scatter(&lAb(counts*numprocs, 0), dim + 1, MPI_DOUBLE, &Ab(partRow - 1, 0), dim + 1, MPI_DOUBLE, masterNode, partComm);

	
}
void gatherResult(const Matrix &partResult,Matrix &result)
{
	int counts = dim / numprocs;
	for (int i = 0; i < counts; ++i)
		MPI_Gather(&partResult(i, 0), 1, MPI_DOUBLE, &result(i*numprocs, 0), 1, MPI_DOUBLE, masterNode, MPI_COMM_WORLD);
	if (counts != partRow)
		MPI_Gather(&partResult(partRow - 1, 0), 1, MPI_DOUBLE, &result(counts*numprocs, 0), 1, MPI_DOUBLE, masterNode, partComm);
}
void parallelGaussSolver(const Matrix &_Ab, Matrix &partResult, Matrix &result)
{
	Matrix Ab = _Ab;
	MPI_Request request;
	MPI_Status status;
	int *loc = new int[dim + 1]{ 0 };
	for (int i = 0; i < dim + 1; ++i)
		loc[i] = i;
	int picked = 0;

	Matrix sta(1, dim + 1);//定义基准行元素

	double s1, s2, dur = 0.0;
	double rec1, rec2, dur2 = 0;

	//群集通信函数实现并行
	/*int p = 0;
	for (int k = 0; k < dim; ++k)
	{
		int source = k % numprocs;
		if (myid == source)
		{
			double max = 0.0;
			for (int j = k; j < dim; ++j)
			{
				if (abs(Ab(p, loc[j])) > max)
				{
					max = abs(Ab(p, loc[j]));
					picked = j;
				}
			}
			
			for (int j = 0; j < dim + 1; ++j)
				sta(0, j) = Ab(p, j);
			p++;
		}
		
		if (myid == masterNode)
			rec1 = MPI_Wtime();
		MPI_Bcast(&sta(0, 0), dim + 1, MPI_DOUBLE, source, MPI_COMM_WORLD);
		MPI_Bcast(&picked, 1, MPI_DOUBLE, source, MPI_COMM_WORLD);
		if (myid == masterNode)
		{
			rec2 = MPI_Wtime();
			dur2 += rec2 - rec1;
		}
		int tmp = loc[k];
		loc[k] = loc[picked];
		loc[picked] = tmp;
		if (myid == masterNode)
			s1 = MPI_Wtime();
		for (int i = p; i < partRow; ++i)
		{
			double t = Ab(i, loc[k]) / sta(0, loc[k]);
			for (int j = k; j < dim + 1; ++j)
				Ab(i, loc[j]) -= t * sta(0, loc[j]);
		}
		if (myid == masterNode)
		{
			s2 = MPI_Wtime();
			dur += s2 - s1;
		}
	}*/
	
	//非阻塞通信实现异步计算
	int p = 0;
	for (int k = 0; k < dim; ++k)
	{
		int source = k % numprocs;
		if (myid == source)
		{
			double max = 0.0;
			for (int j = k; j < dim; ++j)
			{
				if (abs(Ab(p, loc[j])) > max)
				{
					max = abs(Ab(p, loc[j]));
					picked = j;
				}
			}
			MPI_Isend(&Ab(p, 0), dim + 1, MPI_DOUBLE, down, picked, MPI_COMM_WORLD, &request);

			for (int j = 0; j < dim + 1; ++j)
				sta(0, j) = Ab(p, j);

			p++;
		}
		else
		{
			MPI_Irecv(&sta(0, 0), dim + 1, MPI_DOUBLE, up, MPI_ANY_TAG, MPI_COMM_WORLD, &request);
			if (myid == masterNode)
				rec1 = MPI_Wtime();
			MPI_Wait(&request, &status);

			if (myid == masterNode)
			{
				rec2 = MPI_Wtime();
				dur2 += rec2 - rec1;
			}
			picked = status.MPI_TAG;
			if (down != source)
				MPI_Isend(&sta(0, 0), dim + 1, MPI_DOUBLE, down, picked, MPI_COMM_WORLD, &request);
		}

		if (myid == masterNode)
			s1 = MPI_Wtime();
		int tmp = loc[k];
		loc[k] = loc[picked];
		loc[picked] = tmp;


		for (int i = p; i < partRow; ++i)
		{

			double t = Ab(i, loc[k]) / sta(0, loc[k]);

			for (int j = k; j < dim + 1; ++j)
				Ab(i, loc[j]) -= t * sta(0, loc[j]);
		}

		if (myid == masterNode)
		{
			s2 = MPI_Wtime();
			dur += s2 - s1;
		}


	}
	if (myid == masterNode)
	{
		cout << "计算时间：" << dur << endl;
		cout << "通信时间："<<dur2 << endl;
	}

	//并行回代
	int rp = partRow - 1;
	double val = 0;
	for (int k = dim - 1; k >= 0; --k)
	{
		int source = k % numprocs;
		if (myid == source)
		{
			partResult(rp, 0) = Ab(rp, dim) / Ab(rp, loc[k]);
			val = partResult(rp--, 0);
		}
		MPI_Bcast(&val, 1, MPI_DOUBLE, source, MPI_COMM_WORLD);
		for (int j = rp; j >= 0; --j)
			Ab(j, dim) -= val * Ab(j, loc[k]);
	}

	gatherResult(partResult, result);


	//解向量的元素顺序重排
	if (myid == masterNode)
	{
		Matrix rtmp = result;
		for (int i = 0; i < dim; ++i)
			result(loc[i], 0) = rtmp(i, 0);
		//cout << result << endl;
	}
	delete[]loc;
}

int main(int argc, char* argv[])
{
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);//获得进程数
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);//获得当前进程标识号0,1,2,3,....,numprocs - 1
	
	masterNode = 0;

	Matrix lA, lb,lAb,result;//分别是系数矩阵，等号右端常数向量，增广矩阵，解向量
	if (myid == masterNode)
	{
		/*cout << "请输入问题规模，系数矩阵及右端常数向量:" << endl;
		cin >> dim;*/
		lA = Matrix(dim, dim);
		lb = Matrix(dim, 1);
		
		result = Matrix(dim, 1);
		lA.MatrixInit("data_A"), lb.MatrixInit("data_b");//矩阵文件输入
		//lA.ranCreate(),lb.ranCreate(); //自动生成矩阵
		

		lAb = lA.merge(lb, COLUMN);//A，b合并得增广矩阵
		if (numprocs == 1)
		{
			start = MPI_Wtime();
			seqGaussSolver(lA, lb, result);
			finish = MPI_Wtime();
			cout << finish - start << endl;
			MPI_Finalize();
			//system("pause");
			return 0;
		}
	}
	MPI_Bcast(&dim, 1, MPI_INT, masterNode, MPI_COMM_WORLD);//将问题规模广播至各个进程
	parallelInit();//并行初始化工作
	
	Matrix Ab(partRow, dim + 1),partResult(partRow,1);
	distributeTask(lAb,Ab);//将增广矩阵lAb交叉散射到每个进程的Ab中
	start = MPI_Wtime();
	parallelGaussSolver(Ab, partResult, result);
	finish = MPI_Wtime();
	if (myid == masterNode)
	{
		
		cout << finish - start << endl;
		//cout << result << endl;
	}

	
	MPI_Finalize();
	return 0;
}

6.结果

一.使用群集通信的结果：符合预期，计算时间下降，通信略有增加

 二.使用非阻塞通信的结果：仅在进程调用为2的时候符合预期，其余并行部分通信与计算基本1比1；

      将问题规模增至4000，情况依旧。

      

                   

标签
#并行计算 #c++
@Abel
关注
三种线性方程组求解方法的串行算法，和并行算法。包含高斯消元法，LU分解法，追赶法。
06-17
使用CUDA实现的线性方程组的求解并行化。
高斯消去法的MPI实现
06-14
实现高斯消去法解线性方程组的MPI编程，与SSE（或AVX）编程结合，并与Pthread、OpenMP（结合SSE或AVX）版本对比。
pthread & OpenMP 优化高斯消元算法
qq_62617798的博客
07-24 1587
本次实验结合Arm平台NEON指令集和X86平台的AVX指令集来对高斯消元算法进行pthread的多线程优化，并探究不同优化方法（信号量或barrier同步）、编程策略的优化幅度；此外，本次实验还将在x86平台上利用OpenMP进行高斯消元算法的多线程优化，并探究不同任务划分方案的影响。
数值分析—行主元消去法解线性方程组—FORTRAN程序
qq_43630847的博客
11-12 1065
数值分析—行主元消去法解线性方程组—FORTRAN程序 program main implicit none real8,dimension(:
基于高斯消去法解线性方程组（MPI）
11-22
基于高斯消去法解线性方程组（MPI），高斯消去法把Ax=b归约为上三角方程组Tx=c，这样利用回带算法求解x。第i次迭代时，选取i列的最大元素作为主元，主元所在的行称为枢轴行（枢轴行的行数会被标记），枢轴行与第i行进行交换，算法利用枢轴行和第i+1到n-1行各行的倍数将第i列上所有的非零元归约成零。最终将nxn的稠密矩阵化成上三角形，再用回带的方法算出每一个元素的值。
计算机数值方法：高斯消元法解线性方程组
王蒟蒻
04-15 2409
题目 用高斯消元法解下列线性方程组（要求按三位小数计算） * 高斯消元法简介
英特尔oneAPI—高斯消元算法并行化
zouxianyu的博客
05-15 1936
本文基于高斯消元算法串行算法，使用英特尔oneAPI实现该算法的并行化，并结合该算法介绍了oneAPI中一些基本组件的使用方法。在文章的最后设计实验对比了串行算法与并行算法之间的性能差异。
并行计算：MPI总结
热门推荐
JacksonKim的博客
08-03 5万+
一、MPI简介 1.什么是MPI Massage Passing Interface:是消息传递函数库的标准规范，由MPI论坛开发。 一种新的库描述，不是一种语言。共有上百个函数调用接口，提供与C和Fortran语言的绑定 MPI是一种标准或规范的代表，而不是特指某一个对它的具体实现 MPI是一种消息传递编程模型，并成为这种编程模型的代表和事实上的标准 2.MPI的特点 MPI有以下的特点： 消息传递式并行程序设计 指用户必须通过显式地发送和接收消息来实现处理机间的数据交换。 在这种并行编程中，每个
并行计算MPI(一）：MPI入门
sj_lxd的专栏
03-06 1万+
这个学期老板可能让我做并行计算方面的东西并推荐MPI给我学习。私下找了一些资料，这个系列作为学习MPI的一个记录吧，所用的主要有两本书：《高性能计算并行编程技术》，都志辉，清华大学出版社；还有一本《MPI并行程序设计实例》，清华大学出版社。代码有自己写的，书上的，网上的。   1.MPI介绍   说到底，MPI就是一个库，一个抽象了的库，这样在你编写多进程，多线程程序的时候不用关心
我的并行计算之路（二）MPI点对点通信MPI_Send和MPI_Recv
度凡心的博客
08-31 1万+
学习了MPI的点对点技术后，来博客记录一下。先上完整地代码： 　　#include<bits/stdc++.h> #include<mpi.h>using namespace std; int comm_sz=0; int my_rank=0; int arra[]={1,2,3,4,5},arrb[]={10,9,8,7,6},arrc[]={0,0,0,0,0};bool get_input(
我的并行计算之路（一）Ubuntu 16.04下的MPI安装
度凡心的博客
08-30 2万+
现在大四，准备保研了。目前联系好了一个做并行计算的老师，可能以后就走这条路了。人工智能的三大基石之一就是计算大量数据的能力，目前这个技术就是并行计算，无论是用GPU还是大型的集群。然后我现在算是开始学习并行计算，并且准备将学习的经历记录下来。一是供我以后复习用，做个备份，再是想给学习并行计算的朋友一个参考（如果觉得有参考价值的话，摊手）。 　　 下面我就从MPI的安装开始谈起。MPI(Messag
MPI并行计算 - mpi4py
qq_33934427的博客
04-12 2371
MPI并行计算 - mpi4py 博客参考： MPI，OpenMPI 与深度学习 MPI文档（MPI-1.0，2.2，3.0） MPI tutorial introduction (中文版) MPI在Deep Learning的主流时代背景下除了传统计算领域外对DL的应用前景如何，MPI与NCCL的区别在哪？？？ 【科普】超算到底是咋算的？ 代码参考： Python并行编程cookbook（含mpi4py，Lock，Celery） mpi4py 简介 Python使用MPI实现分布式计算（mpi4py
高斯消元法并行MPI解线性方程
08-06
在MPI编译环境下在C源代码的 基础上编写的并行程序：高斯消元法求解线性方程组。
并行计算——MPI
zyw2002的博客
11-10 6398
文章目录MPI 环境搭建下载并安装MPICH2在VS中配置MPI环境CentOS 安装OpenMPI3.1.0MPI 原理以及编程MPI 简介MPI的概念MPI的 通信机制MPI 点对点的通信类型点对点的通信模型点对点通信的常用函数MPI聚合通信课堂编程实例Hello WorldNonblocking send/receive MPI 环境搭建 下载并安装MPICH2 下载并打开“mpich2-1.4.1p1-win-ia64.msi”。 安装64位的MPICH2 如果安装的时候，给出如下的错误提
并行计算（MPI + OpenMP）
weixin_44885334的博客
12-21 6471
Message Passing Interface（MPI）：一种基于信息传递的并行编程技术，定义了一组具有可移植性的编程接口标准（并非一种语言或者接口）。支持点对点通信和广播。MPI 的目标是高性能、大规模性、可移植性，在今天仍为高性能计算的主要模型。OpenMPI 函数库，微软 MPI 文档程序结构（C语言版） 编译： 运行： 数据类型 MPI 数据类型 C语言数据类型 MPI_INT int MPI_FLOAT float MP
MPI安装+CentOs6.5多机环境下MPI并行编程+MPI矩阵并行计算（超详细）
m0_49893559的博客
05-11 6660
MPI 创建多进程，输出进程号和进程数，运行多进程并行例子程序，编程实现大规模向量矩阵并行计算
并行计算OpenMP与MPI环境搭建
mc_故事与你的博客
04-01 4252
目前公开文献可查到比较早提出的定义是Gottlieb等人在1989年《高性能计算》一书中给出的。Gottlieb等认为：并行计算是一种计算类型，其中许多计算或计算过程的执行是同时执行的。并行计算的目的是加快计算速度。我国学者张林波在其编写的教材《并行计算导论》中对并行计算进行了描述性定义：并行计算（parallel computing）是指，在并行机上，将一个应用分解成多个子任务，分配给不同的处理器，各个处理器之间相互协同，并行地执行子任务，从而达到加速求解速度，或者求解应用问题规模的目的。
Linux上搭建并行计算MPI环境
receiling的博客
04-19 2439
最近因为并行计算课程的需求，需要搭建一个多机的MPI运行环境，故将搭建过程分享在这篇博文里。 具体的系统是Ubuntu 18.04，在四台机器上搭建环境。 记四台机器为A，B，C，D，四台机器的IP为IP-A，IP-B，IP-C，IP-D，还有一个要求就是在四台机器上需要同一个用户名，以xiaomin为例。 接下来就是正式开始搭建了： Step1: 创建用户 用ssh以root用户登录四...
MPI并行计算学习笔记2——圆周率pi的并行求解
askjff的博客
05-04 4507
调式环境：VS2017 + MPI 圆周率PI的数值计算公式如下： #include<iostream> #include"mpi.h" #include<ctime> #include<cmath> using namespace std; const int N = 1000000; double start,finish;//计算的起止...
MPI并行计算第一课（2小时入门）
保持写作习惯，完成知识沉淀
01-21 5623
MPI 全名叫 Message Passing Interface，即信息传递接口，作用是可以通过 MPI 可以在不同进程间传递消息，从而可以并行地处理任务，即进行。
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值