高斯约旦消去算法
具体的算法过程参考及函数参数说明可参考
选主元
从第一列开始遍历,选择每列绝对值最大的元素为主元,标记主元所在行以保证选择过的行不应再被选择(这只是一种选择主元的方法,适用于mpi程序)。
初始化代码:
//用到的变量名和数组
int rank, p;
const int n = 1000;//矩阵为n*n矩阵
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]; // 存储系数矩阵和常数项
double(*z)[n + 1] = new double[m][n + 1]; // 备份的系数矩阵,用于计算误差
int* pivot = new int[m]; // 记录每个进程处理的行的主元所在的列
double x[n]; // 存储最终的解
int which[n]; // 记录每个列的主元所在的进程号
double picked_row[n + 1];
// 选中的主元行,将存储当前遍历的主元行所在进程发送过来的主元行,用来消去使得选定主元所在列的其他元素都为0;
// 初始化
for (int i = 0; i < m; i++)
pivot[i] = -1;//-1表示该行还未被选定主元
srand(rank * 100);
for (int i = 0; i < m; i++)
for (int j = 0; j < n + 1; j++)
a[i][j] = z[i][j] = rand() * 1.0 / RAND_MAX;//随机生成a,z
选择主元(外层对每列的遍历省略):
//k为外层循环计数器,表示第k列
double local_max = 0.0;
int local_row = -1;
for (int i = 0; i < m; i++) {
if (pivot[i] == -1 && fabs(a[i][k]) > local_max) {//该行需未被选定主元
local_max = fabs(a[i][k]);
local_row = i;
}
}
struct {
double value;
int rank;
} local_data, global_data;
local_data.value = local_max;
local_data.rank = rank;
// 记录每个进程找到的局部最大值
MPI_Allreduce(&local_data, &global_data, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
//Mpi的规约函数,将从每个进程关于第k列所给出的本地最大值,返回全局的最大值,存到每个进程的global_data中,此时global_data中存储的是主元行的最大值和所在进程号。
主元行需执行归一化,消去:
将主元行同时除以主元,得到归一化后的主元行,其他各行需减去这个主元行的倍数,最后主元所在这一列只有主元不为0。
归一化:
if (global_data.rank == rank) {
// 主元所在的进程将主元行广播给所有进程
double f = a[local_row][k];
for (int j = 0; j < n + 1; j++)
a[local_row][j]/=f;
//归一化
MPI_Bcast(a[local_row], n + 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
//广播至各进程picked_row。
pivot[local_row] = k;
//表示该行以及选定了主元在第k列
for (int j = 0; j < n + 1; j++)
picked_row[j] = a[local_row][j];
//对于本进程的其他行来说也要执行消去,本进程的picked_row也需赋值。
which[k] = rank;//记录第k列的主元,也就是最终解x[k]在哪个进程里。
MPI_Bcast(&which[k], 1, MPI_INT, rank, MPI_COMM_WORLD);//广播
}
else {
MPI_Bcast(picked_row, n + 1, MPI_DOUBLE, global_data.rank, MPI_COMM_WORLD);
//接收
MPI_Bcast(&which[k], 1, MPI_INT, global_data.rank, MPI_COMM_WORLD);
//接收
}
归一化后消去:
// 消去k列的其他元素
for (int i = 0; i < m; i++) {
//不是主元行的行执行消去
if (!(i==local_row&&rank== global_data.rank)) {
double factor = a[i][k];
//该行主元所在列的元素,将决定消去多少倍的归一化后的主元行
for (int j = k; j < n + 1; j++) {
//前k列执行过消去后,传过来的picked_row(0~k-1),肯定为0,所以可以从k开始。
a[i][j] -= factor * picked_row[j];
}
}
}
求解
历经k次循环后可想而知,矩阵A的形式将会变成做数次行变换后的单位矩阵,x的解也就是现在的a[0~n-1][n]列,但是位置上却不是x[i]=a[i][n]的对应关系,这时which数组将帮我们定位x[k]所在的进程。
代码:
//这里对列的遍历已经结束。此k非彼k。
for (int k = 0; k <n; k++) {
if(rank == which[k]) {//定位进程
for (int i = 0; i < m; i++) {
if (a[i][k] == 1)//定位解所在的行
x[k] = a[i][n];
}
MPI_Bcast(&x[k], 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);//广播
}
else
MPI_Bcast(&x[k], 1, MPI_DOUBLE, which[k], MPI_COMM_WORLD);//其余进程接收。
}
完整代码
#include "mpi.h"
#include <iostream>
#include <cmath>
#include <algorithm>
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]; // 存储系数矩阵和常数项
double(*z)[n + 1] = new double[m][n + 1]; // 备份的系数矩阵,用于计算误差
int* pivot = new int[m]; // 记录每个进程处理的行的主元所在的列
double* max_coefficient = new double[p]; // 每个进程的最大系数
double x[n]; // 存储最终的解
int which[n]; // 记录每个列的主元所在的行
double picked_row[n + 1]; // 选中的主元行
// 初始化
for (int i = 0; i < m; i++)
pivot[i] = -1;
srand(rank * 100);
for (int i = 0; i < m; i++)
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++) {
// 寻找第k列的主元
double local_max = 0.0;
int local_row = -1;
for (int i = 0; i < m; i++) {
if (pivot[i] == -1 && fabs(a[i][k]) > local_max) {
local_max = fabs(a[i][k]);
local_row = i;
}
}
struct {
double value;
int rank;
} local_data, global_data;
local_data.value = local_max;
local_data.rank = rank;
// 记录每个进程找到的局部最大值
//将每个进程关于第k列的最大系数的最大值传入global_data
MPI_Allreduce(&local_data, &global_data, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
if (global_data.rank == rank) {
// 主元所在的进程将主元行广播给所有进程
double f = a[local_row][k];
for (int j = 0; j < n + 1; j++)
a[local_row][j]/=f;
MPI_Bcast(a[local_row], n + 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
pivot[local_row] = k;
for (int j = 0; j < n + 1; j++)
picked_row[j] = a[local_row][j];
which[k] = rank;
MPI_Bcast(&which[k], 1, MPI_INT, rank, MPI_COMM_WORLD);
}
else {
MPI_Bcast(picked_row, n + 1, MPI_DOUBLE, global_data.rank, MPI_COMM_WORLD);
MPI_Bcast(&which[k], 1, MPI_INT, global_data.rank, MPI_COMM_WORLD);
}
// 消去k列的其他元素
for (int i = 0; i < m; i++) {
if (i==local_row&&rank== global_data.rank) {
}
else {
double factor = a[i][k];
for (int j = k; j < n + 1; j++) {
a[i][j] -= factor * picked_row[j];
}
}
}
}
// 回代求解
for (int k = 0; k <n; k++) {
if(rank == which[k]) {
for (int i = 0; i < m; i++) {
if (a[i][k] == 1)
x[k] = a[i][n];
}
MPI_Bcast(&x[k], 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
}
else
MPI_Bcast(&x[k], 1, MPI_DOUBLE, which[k], 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;
}