高斯列主元消元法mpi实现

列主元法mpi实现

列主元消元法在求解线性方程组时很好的解决了因为计算机本身字长的限制而产生的问题。本文中使用mpi将矩阵分行处理，并将通信量最小化。这是课程作业，学弟学妹们参考参考就好

直接放代码

#include<iostream>#include"mpi.h"#include<fstream>#include<algorithm>#include<cmath>using namespace std;const double eps=1e-10;const int N=1000;int main(){    double **a=new double*[N];    double *answer=new double[N];    int *local=new int[N+1];    double max;    int i,j,myid,numprocs,madex;    for(i=0;i<N;i++)        a[i]=new double[N+1];//malloc        /*local 数组用于对应主元的位置和原来的位置的映射关系*/    for(i=0;i<N+1;i++)        local[i]=i;    ifstream fin("matrix.dat");//fscanf    for(i=0;i<N;i++)        for(j=0;j<N+1;j++)            fin>>a[i][j];    MPI_Init(0,0);    MPI_Comm_rank(MPI_COMM_WORLD,&myid);    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);    int numtasks=N/numprocs;    for(j=0;j<N;j++){        if(myid==j%numprocs){            max=madex=0;            for(int k=0;k<N;k++)                if(local[k]>=j&&abs(a[j][local[k]])>=max){                    max=abs(a[j][local[k]]);                    madex=k;                }        }        MPI_Bcast(&madex,1,MPI_INT,j%numprocs,MPI_COMM_WORLD);        MPI_Bcast(a[j],N+1,MPI_DOUBLE,j%numprocs,MPI_COMM_WORLD);        swap(local[j],local[madex]);//交换两个部分的值        for(i=myid;i<N;i+=numprocs){            if(i!=j){                double temp=a[i][local[j]]/(a[j][local[j]]+eps);                for(int k=0;k<N+1;k++){                    a[i][k]-=a[j][k]*(temp);                }            }        }    }    for(j=0;j<N;j++){        if(myid==j%numprocs){            answer[local[j]]=a[j][N]/(a[j][local[j]]+eps);        }        MPI_Bcast(answer+local[j],1,MPI_DOUBLE,j%numprocs,MPI_COMM_WORLD);    }    if(myid==0){        for(j=0;j<N;j++){            cout<<"X["<<j<<"]="<<answer[j]<<endl;        }    }    MPI_Finalize();}