MPI列主元高斯消元法的应用
20180914重新编辑并对原算法进行修正
想不到我也会有朝一日当助教,又重新写了一次mpi的高斯消元法,发现思维已经没有以前灵活了,感慨良多。
列主元高斯消元法
高斯消元法思想非常简单,学过线性代数的基本都对此比较熟悉了。高斯消元法是求解Ax=b的一种方法,列主元高斯消元法是对高斯消元法的一种拓展,克服了由于机器字长限制,将小主元误差放大的后果,串行基本步骤如下:
- 初始化映射数组,i=0
- 将A与b列成增广矩阵
- 通过映射数组开始寻找第i列主元,并记录主元到映射数组
- 通过主元行消去映射数组中未定义的行,让在映射数组中未定义的行的第i个元素为0
- i++并重复步骤3 4 直到i=N-1,此时已经化成上三角矩阵了
- 将上三角矩阵从下往上消去(前向替换算法)获得c
串行程序中,最花时间的就是4步,复杂度为O(n^2),我们尝试将这部分并行掉
mpi并行列主元消元法思路
考虑到消元部分前面的处理数据较多(上三角上面的部分),后面处理的数据较少(上三角尖的部分),为了能让每个核处理的数据量能大致相似,即使列主元数据也具有可能有序性,我们采用交叉分配的方式分配内存,整个流程如下:
- 给0号进程输入增广矩阵
- 将增广矩阵按交叉分配的方式(如总共有4个进程,1号进程就负责1,5,9,13。。。。行)先将矩阵重排列,让同一个进程需要的内存先排列成连续的内存,然后再讲这个矩阵分发出去给每个进程,保证每个进程只用到自己应该用的部分
- 寻找主元行,这里注意,主元行的寻找需要通过MPI的规约函数,规约到当前列的最大值后,与进程中该列的值对比可以得到最大的主元行,然后将主元行的行号广播出去。
- 广播主元行并消元,消元法同串行
- 获得“上三角”矩阵后再收集回0号进程,在0号进程中进行前向替换运算,并在0进程中输出结果
这种方式,时间和空间都有比较不错的效果
下面是源代码
/**yhf at YUQUAN**/
#include"mpi.h"
#include<iostream>
#include<string>
#include<vector>
#include<cmath>
#include<cstdlib>
#include<algorithm>
using namespace std;
const int N = 2000;
int myid, numproc;
void copyMemory(const double *M, double *A) {//复制矩阵使得同个进程中需要的内存是连续的
for (int i = 0; i < N*(N + 1); i++)
A[((i / (N + 1)) % numproc)*( (N * (N + 1)) / numproc) + (N+1) * ((i / (N + 1)) / numproc) + i % (N + 1)] = M[i];
}
int main(int argc,char* argv[]) {
MPI_Init(&argc,&argv);
double t1, t2;
t1 = MPI_Wtime();
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
MPI_Comm_size(MPI_COMM_WORLD,&numproc);
double *M, *c, *M_trans,*A,*tempRow;
int *map = new int[N]; int *rmap = new int[N];
for (int i = 0; i < N; i++) {
map[i] = -1; rmap[i] = -1;
}
if (myid == 0) {//内存分配
M = new double[N*(N+1)];
M_trans = new double[N*(N + 1)];
for (int i = 0; i < N*(N + 1); i++) {
M[i] = rand() % 10;
}
copyMemory(M, M_trans);
A = new double[N*(N + 1) / numproc + 1];
c = new double[N];
}
else {
M = new double[1];
M_trans = new double[1];
A = new double[N*(N + 1) /numproc+1];
c = new double[N];
}
tempRow = new double[N + 1];
MPI_Scatter(M_trans, N*(N+1) / numproc,MPI_DOUBLE,A,N*(N+1)/numproc,MPI_DOUBLE,0,MPI_COMM_WORLD);//分发矩阵
delete M_trans;
double local_max_value, global_max_value;
int local_max_id, max_proc_id, max_pro_id_1, global_max_id;
for (int row = 0; row < N-1; row++) {
local_max_value = -1e10; local_max_id = -1; max_pro_id_1 = -1, global_max_id=-1;
for (int i = 0; i < N / numproc; i++) {
if (rmap[i * numproc + myid] < 0) {
if (A[i*(N + 1) + row] > local_max_value) {
local_max_value = A[i*(N + 1) + row];
local_max_id = i;
}
}
}//找到对应进程的主元
MPI_Allreduce(&local_max_value, &global_max_value, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);//将每个进程的主元找到一个最大值
if (global_max_value > local_max_value-0.00001&&global_max_value < local_max_value+0.00001) {//对比最大值找到主元
max_pro_id_1 = myid;
for (int i = 0; i < N + 1; i++) {
tempRow[i] = A[local_max_id*(N + 1) + i];
}
global_max_id= local_max_id * numproc + myid;
}
MPI_Allreduce(&max_pro_id_1,&max_proc_id,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);//广播主元所在的进程便于下面的操作
MPI_Bcast(&max_proc_id,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(tempRow,N+1,MPI_DOUBLE,max_proc_id,MPI_COMM_WORLD);//广播主元行
MPI_Bcast(&global_max_id,1,MPI_INT,max_proc_id,MPI_COMM_WORLD);//广播主元id
map[row] = global_max_id; rmap[global_max_id] = row;
for (int i = 0; i < N/numproc; i++) {//求解过程
if (rmap[i * numproc + myid] < 0) {
double temp = A[i*(N + 1)+row] / tempRow[row];
for (int j = row; j < N + 1; j++) {
A[i*(N + 1) + j] -= tempRow[j] * temp;
}
}
}
}
MPI_Gather(A,N*(N+1)/numproc,MPI_DOUBLE,M,N*(N+1)/numproc,MPI_DOUBLE,0,MPI_COMM_WORLD);//求解完三角阵后规约起来
if (myid == 0) {
for (int i = 0; i < N; i++) {
if (rmap[i] == -1) {
map[N - 1] = i;
}
}
printf("\n");
}
if (myid == 0) { //输出
for (int i = N-1; i >=0; i--) {
int index = map[i] % numproc*(N / numproc) + map[i] / numproc;
for (int j = N-1; j >i; j--) {
M[index*(N+1)+N] -= c[j] * M[index*(N + 1) + N - j];
}
c[i] = M[index*(N + 1) + N] / M[index*(N+1)+i];
}
for (int i = 0; i < N; i++) {
printf("C[%d]=%.2f\n",i, c[i]);
}
}
t2 = MPI_Wtime();
printf("it take %.2lfs\n", t2-t1);
MPI_Finalize();
}
经测试,在N=2000左右的时候,4核的加速比可以达到3.5,这说明该程序还是很好的执行了并行的模块
不足
- 需要在主进程额外地开辟一块内存
- 只可用于计算核数能整除行数的情况