MPI 跟OPENMP 混合编程 实现矩阵LU分解
MPI 和OPENMP 混合编程 实现矩阵LU分解
LU分解
将系数矩阵A转变成等价两个矩阵L和U的乘积 ,其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LU,且分解唯一。其中L是单位下三角矩阵,U是上三角矩阵。
方法:
使用openMP和MPI混合编程现实
代码如下:
#include "stdio.h" #include "stdlib.h" #include "mpi.h" #include "omp.h" /***************MPI openMP 混合实现LU分解*******************/ /************ Yingfeng Chen ******************************/ /********************************************************/ #define a(x,y) a[x*M+y] /*A为M*M矩阵*/ #define A(x,y) A[x*M+y] #define l(x,y) l[x*M+y] #define u(x,y) u[x*M+y] #define floatsize sizeof(float) #define intsize sizeof(int) int M,N; int m; float *A; int my_rank; int p; MPI_Status status; void fatal(char *message) { printf("%s\n",message); exit(1); } void Environment_Finalize(float *a,float *f) { free(a); free(f); } int main(int argc, char **argv) { int i,j,k,my_rank,group_size; int i1,i2; int v,w; float *a,*f,*l,*u; FILE *fdA; 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) { puts("The input is error!"); exit(0); } A=(float *)malloc(floatsize*M*M); for(i = 0; i < M; i ++) for(j = 0; j < M; j ++) fscanf(fdA, "%f", A+i*M+j); fclose(fdA); } /*0号进程将M广播给所有进程*/ MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); m=M/p; if (M%p!=0) m++; /*分配至各进程的子矩阵大小为m*M*/ a=(float*)malloc(floatsize*m*M); /*各进程为主行元素建立发送和接收缓冲区*/ f=(float*)malloc(floatsize*M); /*0号进程为l和u矩阵分配内存,以分离出经过变换后的A矩阵中的l和u矩阵*/ if (my_rank==0) { l=(float*)malloc(floatsize*M*M); u=(float*)malloc(floatsize*M*M); } /*0号进程采用行交叉划分将矩阵A划分为大小m*M的p块子矩阵,依次发送给1至p-1号进程*/ if (a==NULL) fatal("allocate error\n"); 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++) if ((i%p)!=0) { i1=i%p; i2=i/p+1; MPI_Send(&A(i,0),M,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); } for(i=0;i<m;i++) for(j=0;j<p;j++) { /*j号进程负责广播主行元素*/ if (my_rank==j) { v=i*p+j; for (k=v;k<M;k++) f[k]=a(i,k); MPI_Bcast(f,M,MPI_FLOAT,my_rank,MPI_COMM_WORLD); } else { v=i*p+j; MPI_Bcast(f,M,MPI_FLOAT,j,MPI_COMM_WORLD); } /*编号小于my_rank的进程(包括my_rank本身)利用主行对其第i+1,…,m-1行数据做行变换*/ /*********MPI 并行优化 ********/ if (my_rank<=j) { #pragma omp parallel shared(a,f,k,v,m) private(k,w) { #pragma omp for for(k=i+1;k<m;k++) { a(k,v)=a(k,v)/f[v]; } for(k=i+1;k<m;k++) { #pragma omp for for(w=v+1;w<M;w++) a(k,w)=a(k,w)-f[w]*a(k,v); } } } /*编号大于my_rank的进程利用主行对其第i,…,m-1行数据做行变换*/ if (my_rank>j) { #pragma omp parallel shared(a,f,k,v,m) private(k,w) { #pragma omp for for(k=i;k<m;k++) { a(k,v)=a(k,v)/f[v]; } for(k=i;k<m;k++) { #pragma omp for for(w=v+1;w<M;w++) a(k,w)=a(k,w)-f[w]*a(k,v); } } } } /*0号进程从其余各进程中接收子矩阵a,得到经过变换的矩阵A*/ if (my_rank==0) { for(i=0;i<m;i++) for(j=0;j<M;j++) A(i*p,j)=a(i,j); } if (my_rank!=0) { for(i=0;i<m;i++) MPI_Send(&a(i,0),M,MPI_FLOAT,0,i,MPI_COMM_WORLD); } else { for(i=1;i<p;i++) for(j=0;j<m;j++) { MPI_Recv(&a(j,0),M,MPI_FLOAT,i,j,MPI_COMM_WORLD,&status); for(k=0;k<M;k++) A((j*p+i),k)=a(j,k); } } if (my_rank==0) { omp_set_num_threads(M); #pragma omp parallel shared(l,u,A,M) private(i,j) { #pragma omp for for(i=0;i<M;i++) { #pragma omp for for(j=0;j<M;j++) { if(i>j) { l(i,j)=A(i,j); u(i,j)=0.0; } else if(i<j) { u(i,j)=A(i,j); l(i,j)=0.0; } else { u(i,j)=0.0; l(i,j)=1.0; } } } } for(i=0;i<M;i++) for(j=0;j<M;j++) { if (i>j) l(i,j)=A(i,j); else u(i,j)=A(i,j); } printf("Input of file \"dataIn.txt\"\n"); printf("%d\t %d\n",M, N); for(i=0;i<M;i++) { for(j=0;j<N;j++) printf("%f\t",A(i,j)); printf("\n"); } printf("\nOutput of LU operation\n"); printf("Matrix L:\n"); for(i=0;i<M;i++) { for(j=0;j<M;j++) printf("%f\t",l(i,j)); printf("\n"); } printf("Matrix U:\n"); for(i=0;i<M;i++) { for(j=0;j<M;j++) printf("%f\t",u(i,j)); printf("\n"); } } MPI_Finalize(); Environment_Finalize(a,f); return(0); }编译命令:mpicc -cc=omcc -o 执行文件 源文件
执行命令:mpirun -np
节点数 ./执行文件