MPI.rar_Gauss全主元消去_MPI_全主元消去

#include "stdio.h" #include "stdlib.h" #include "mpi.h" #include "math.h" #define a(x,y) a[x*M+y] #define b(x) b[x] //A为M*M的系数矩阵 #define A(x,y) A[x*M+y] #define B(x) B[x] #define floatsize sizeof(float) #define intsize sizeof(int) int M; int N; int m; float *A; float *B; int my_rank; int p; int l; MPI_Status status; void fatal(char *message) { printf("%s\n",message); exit(1); } void Environment_Finalize(float *a,float *b,float *x,float *f) { free(a); free(b); free(x); free(f); } int main(int argc,char **argv) { int i,j,t,k,my_rank,group_size; int i1,i2; int v,w; float temp; int tem; float *sum; float *f; float lmax; float *a; float *b; float *x; int *shift; FILE *fdA,*fdB; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&group_size); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); p=group_size; if(my_rank==0) { fdA=fopen("dataIn.txt","r"); fscanf(fdA,"%d %d",&M,&N); if(M!=N-1) { printf("the input is wrong\n"); exit(1); } A=(float *)malloc(floatsize*M*M); B=(float *)malloc(floatsize*M); for(i=0;i<M;i++) { for(j=0;j<M;j++) fscanf(fdA,"%f",A+i*M+j); fscanf(fdA,"%f",B+i); } fclose(fdA); } //0号处理器将M广播给所有处理器 MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); m=M/p; if (M%p!=0) m++; //各处理器为主行元素建立发送和接收缓冲区(M+1) f=(float*)malloc(sizeof(float)*(M+1)); //分配至各处理器的子矩阵大小为m*M a=(float*)malloc(sizeof(float)*m*M); //分配至各处理器的子向量大小为m b=(float*)malloc(sizeof(float)*m); sum=(float*)malloc(sizeof(float)*m); x=(float*)malloc(sizeof(float)*M); shift=(int*)malloc(sizeof(int)*M); if (a==NULL||b==NULL||f==NULL||sum==NULL||x==NULL||shift==NULL) fatal("allocate error\n"); for(i=0;i<M;i++) shift[i]=i; //0号处理器采用行交叉划分将矩阵A划分为大小为m*M的p块子矩阵，将B划分为大小为m的p块子向量，依次发送给1-p-1号处理器 if(my_rank==0) { for(i=0;i<m;i++) for(j=0;j<M;j++) a(i,j)=A(i*p,j); for(i=0;i<m;i++) b(i)=B(i*p); for(i=0;i<M;i++) if((i%p)!=0) { i1=i%p; i2=i/p+1; MPI_Send(&A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD); MPI_Send(&B(i),1,MPI_FLOAT,i1,i2,MPI_COMM_WORLD); } } else { for(i=0;i<m;i++) { MPI_Recv(&a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status); MPI_Recv(&b(i),1,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status); } } //消去 for(i=0;i<m;i++) for(j=0;j<p;j++) { //j号处理器负责广播主行元素 if(my_rank==j) { //主元素在原系数矩阵A中的行号和列号为v v=i*p+j; lmax=a(i,v); l=v; //在同行的元素中寻找最大元并确定最大元并确定最大元所在的列 for(k=v+1;k<M;k++) if(fabs(a(i,k))>lmax) { lmax=a(i,k); l=k; } //列交换 if(l!=v) { for(t=0;t<m;t++) { temp=a(t,v); a(t,v)=(t,l); a(t,l)=temp; } tem=shift[v]; shift[v]=shift[l]; shift[l]=tem; } //归一化 for(k=v+1;k<M;k++) a(i,k)=a(i,k)/a(i,v); b(i)=b(i)/a(i,v); a(i,v)=l; for(k=v+1;k<M;k++) f[k]=a(i,k); f[M]=b(i); //发送归一化的主行 MPI_Bcast(&f[0],M+1,MPI_FLOAT,my_rank,MPI_COMM_WORLD); //发送主行中主元素所在的列号 MPI_Bcast(&l,l,MPI_INT,my_rank,MPI_COMM_WORLD); } else { v=i*p+j; //接收归一化后的主行 MPI_Bcast(&f[0],M+1,MPI_FLOAT,j,MPI_COMM_WORLD); //接收主行中主元素所在的列号 MPI_Bcast(&l,l,MPI_INT,j,MPI_COMM_WORLD); //列交换 if(l!=v) { for(t=0;t<m;t++) { temp=a(t,v); a(t,v)=a(t,l); a(t,l)=temp; } tem=shift[v]; shift[v]=shift[l]; shift[l]=tem; } } //编号小于j的处理器(包括j)利用主行对其第i+1,..,m-1行数据做行变换 if(my_rank<=j) for(k=i+1;k<m;k++) { for(w=v+1;w<M;w++) a(k,w)=a(k,w)-f[w]*a(k,v); b(k)=b(k)-f[M]*a(k,v); } //编号大于j的处理器利用主行对其第i,..,m-1行数据做行变换 if(my_rank>j) for(k=i;k<m;k++) { for(w=v+1;w<M;w++) a(k,w)=a(k,w)-f[w]*a(k,v); b(k)=b(k)-f[M]*a(k,v); } } for(i=0;i<m;i++) sum[i]=0.0; //回代 for(i=m-1;i>=0;i--) for(j=p-1;j>=0;j--) if(my_rank==j) { x[i*p+j]=(b(i)-sum[i])/a(i,i*p+j); MPI_Bcast(&x[i*p+j],1,MPI_FLOAT,my_rank,MPI_COMM_WORLD); for(k=0;k<i;k++) { sum[k]=sum[k]+a(k,i*p+j)*x[i*p+j]; } } else { MPI_Bcast(&x[i*p+j],1,MPI_FLOAT,j,MPI_COMM_WORLD); if(my_rank>j) for(k=0;k<i;k++) sum[k]=sum[k]+a(k,i*p+j)*x[i*p+j]; if(my_rank<j) for(k=0;k<=i;k++) sum[k]=sum[k]+a(k,i*p+j)*x[i*p+j]; } if(my_rank!=0) for(i=0;i<m;i++) MPI_Send(&x[i*p+my_rank],1,MPI_FLOAT,0,i,MPI_COMM_WORLD); else //0号处理器从其余各处理器中接收子解向量x for(i=1;i<p;i++) for(j=0;j<m;j++) MPI_Recv(&x[j*p+i],1,MPI_FLOAT,i,j,MPI_COMM_WORLD,&status); if(my_rank==0) { printf("Input of file \"dataIn.txt\"\n"); printf("%d\t%d\n",M,N); for(i=0;i<M;i++) { for(j=0;j<M;j++) printf("%f\t",A(i,j)); printf("%f\n",B(i)); } printf("\nOutput of solution\n"); for(k=0;k<M;k++) { for(i=0;i<M;i++) { if(shift[i]==k) printf("x[%d]=%f\n",k,x[i]); } } } MPI_Finalize(); Environment_Finalize(a,b,x,f); return 0; } /* #include "stdio.h" #include "stdlib.h" #include "mpi.h" #include "math.h" #define a(x,y) a[x*M+y] #define b(x) b[x] #define A(x,y) A[x*M+y] #define B(x) B[x] #define floatsize sizeof(float) #define intsize sizeof(int) int M; int N; int m; float *A; float *B; double starttime; double time1; double time2; int my_rank; int p; int l; MPI_Status status; void fatal(char *message) { printf("%s\n",message); exit(1); } void Environment_Finalize(float *a,float *b,float *x,float *f) { free(a); free(b); free(x); free(f); } int main(int argc, char **argv) { int i,j,t,k,my_rank,group_size; int i1,i2; int v,w; float temp; int tem; float *sum; float *f; float lmax; float *a; float *b; float *x; int *shift; FILE *fdA,*fdB; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&group_size); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); p=group_size; if (my_rank==0) { starttime=MPI_Wtime(); fdA=fopen("dataIn.txt","r"); fscanf(fdA,"%d %d", &M, &N); if (M != N-1) { printf("the input is wrong\n"); exit(1); } A=(float *)malloc(floatsize*M*M); B=(float *)malloc(floatsize*M); for(i = 0; i < M; i++) { for(j = 0; j < M; j++) { fscanf(fdA,"%f", A+i*M+j); } fscanf(fdA,"%f", B+i); } fclose(fdA); } //0号处理器将M广播给所有处理器 MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); m=M/p; if (M%p!=0) m++; //各处理器为主行元素建立发送和接收缓冲区(M+1) f=(float*)malloc(sizeof(float)*(M+1)); //分配至各处理器的子矩阵大小为m*M a=(float*)malloc(sizeof(float)*m*M); //分配至各处理器的子向量大小为m b=(float*)malloc(sizeof(float)*m); sum=(