高斯-约旦消去法的MPI并行程序C++实现

高斯约旦消去算法

具体的算法过程参考及函数参数说明可参考

高斯约旦算法

MPI_Bcast说明

MPI_Allreduce说明

选主元

从第一列开始遍历，选择每列绝对值最大的元素为主元，标记主元所在行以保证选择过的行不应再被选择（这只是一种选择主元的方法，适用于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;
}