高斯列主元消元法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*
;
double *answer=new double
;
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]
/(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();
}