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

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 - 1masterNode = 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，情况依旧。

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

tornado学习笔记day04-执行顺序
响应输出 -> write 原型 self.write()函数源码中是这样定义的 def write(self, chunk: Union[str, bytes, dict]) -> N ...
【温故知新】CSS学习笔记（行高简介）
行高简介 [行高的测量] 真正的行高有如下图所示的四条线组成(顶线.中线.基线.底线). 顶线:以文字的最高点为准: 中线:以文字的最高点和最低点中间取一条线为中线: 基线:以文字中最短的底线为基准: ...
【MySQL 】学习笔记千行总结
/* Windows服务 */ -- 启动MySQLnet start mysql -- 创建Windows服务sc create mysql binPath= mysqld_bin_path(注意: ...
云计算学习笔记006---运行hadoop的例子程序：统计字符--wordcount例子程序
04-运行wordcount例子程序下面可以看下hadoop的例子程序: hadoop-0.20.2-examples.jar 注意这里用到的wordcount.txt中的内容为: hello ha ...
python学习笔记15-执行环境
可调用对象代码对象一.可调用对象 Python 有4 种可调用对象:函数,方法,类,以及一些类的实例.记住这些对象的任何引用或者别名都是可调用的. 1.函数 python 有3 种不同类型函数对象 ...
[Shell学习笔记] 命令行下的高级网络工具cURL命令
原文: http://www.1987.name/365.html Linux curl命令是一个利用URL规则在命令行下工作的文件传输工具.它支持文件的上传和下载,所以是综合传输工具,但按传统,习惯 ...
[前端项目学习笔记] 200行代码网站首页轮播实现（html,css,js）
目录 1.设置基本布局 2.添加轮播按钮 3.轮播代码初步实现 4.给按钮添加点击事件实现轮播 5.添加圆点轮播 6.将列表项替换为图片,并给图片加上超链接 7.最终效果 1.设置基本布局整体布局我 ...
jmh学习笔记-State共享对象
系列文章目录 jmh学习笔记-源代码编译与bench mode jmh学习笔记-State共享对象 jmh学习笔记-State共享对象前后置方法 jmh学习笔记-代码清除 jmh学习笔记-常量折叠 j ...
[学习笔记]高斯消元求解两种特殊问题(带状矩阵/主元法)
本文章是[学习笔记]概率与期望进阶的一部分由于时间问题我写的比较简略,所以我把大佬的总结链接贴上来了(应该没什么吧qwq). 1 概述最常见的当然是随机游走问题了- • fu=∑pu,v∗(fv+ ...

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

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

最新文章

热门文章