MPI学习笔记(二):矩阵相乘的两种实现方法
mpi矩阵乘法(C=αAB+βC)
最近领导让把之前安装的软件lapack、blas里的dgemm运算提取出来独立作为一套程序,然后把这段程序改为并行的,并测试一下进程规模扩展到128时的并行效率。我发现这个是dgemm.f文件,里面主要是对C=αAB+βC的实现,因此在此总结一下MPI的矩阵乘法使用。
其主要思想:是把相乘的矩阵按行分解(任务分解),分别分给不同的进程,然后在汇总到一个进程上,在程序上实现则用到了主从模式,人为的把进程分为主进程和从进程,主进程负责对原始矩阵初始化赋值,并把数据均匀分发(为了负载均衡)到从进程上进行相乘运算,主要用到的知识是MPI点对点通信和组通信的机制。
一、使用简单的MPI_Send和MPI_Recv实现
#include <stdio.h> #include "mpi.h" #include <stdlib.h> #include "functions.h" #define M 1000 // 矩阵维度 #define N 1100 #define K 900 int main(int argc, char **argv) { int my_rank,comm_sz,line; double start, stop; //计时时间 MPI_Status status; char Processorname[20]; double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C; double alpha=2,beta=2; // 系数C=aA*B+bC MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); line=M/comm_sz; // 每个进程分多少行数据 Matrix_A=(double*)malloc(M*N*sizeof(double)); Matrix_B=(double*)malloc(N*K*sizeof(double)); Matrix_C=(double*)malloc(M*K*sizeof(double)); buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的数据 buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的数据 ans=(double*)malloc(line*K*sizeof(double)); // 临时保存部分数据计算结果 // 给矩阵A B,C赋值 if(my_rank==0){ start=MPI_Wtime(); for(int i=0;i<M;i++){ for(int j=0;j<N;j++) Matrix_A[i*N+j]=i+1; } for(int i=0;i<N;i++){ for(int j=0;j<K;j++) Matrix_B[i*K+j]=j+1; } for(int i=0;i<M;i++){ for(int j=0;j<K;j++) Matrix_C[i*K+j]=1; } // 输出A,B,C /*Matrix_print(Matrix_A,M,N); Matrix_print(Matrix_B,N,K); Matrix_print(Matrix_C,M,K); */ /*将矩阵广播出去*/ for(int i=1;i<comm_sz;i++){ MPI_Send(Matrix_A+(i-1)*line*N,line*N,MPI_DOUBLE,i,66,MPI_COMM_WORLD); MPI_Send(Matrix_C+(i-1)*line*K,line*K,MPI_DOUBLE,i,99,MPI_COMM_WORLD); } MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 接收从进程的计算结果 for(int p=1;p<comm_sz;p++){ MPI_Recv(ans,line*K,MPI_DOUBLE,p,33,MPI_COMM_WORLD,&status); for(int i=0;i<line;i+=comm_sz) for(int j=0;j<K;j++) Matrix_C[((p-1)*line+i)*K+j]=ans[i*K+j]; } // 计算A剩下的行数据 for(int i=(comm_sz-1)*line;i<M;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++) temp+=Matrix_A[i*N+p]*Matrix_B[p*K+j]; Matrix_C[i*K+j]=alpha*temp+beta*Matrix_C[i*K+j]; } } //Matrix_print(Matrix_C,M,K); stop=MPI_Wtime(); printf("rank:%d time:%lfs\n",my_rank,stop-start); free(Matrix_A); free(Matrix_B); free(Matrix_C); free(buffer_A); free(buffer_C); free(ans); } else{ //接收广播的数据 MPI_Recv(buffer_A,line*N,MPI_DOUBLE,0,66,MPI_COMM_WORLD,&status); MPI_Recv(buffer_C,line*K,MPI_DOUBLE,0,99,MPI_COMM_WORLD,&status); MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); //计算乘积结果,并将结果发送给主进程 for(int i=0;i<line;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++){ temp+=buffer_A[i*N+p]*Matrix_B[p*K+j]; } ans[i*line+j]=alpha*temp+beta*buffer_C[i*K+j]; } } MPI_Send(ans,line*K,MPI_DOUBLE,0,33,MPI_COMM_WORLD); } MPI_Finalize(); return 0; }
二、使用较高级的MPI_Scatter和MPI_Gather实现
#include <stdio.h> #include "mpi.h" #include <stdlib.h> #include "functions.h" #define M 1200 // 矩阵维度 #define N 1000 #define K 1100 int main(int argc, char **argv) { int my_rank,comm_sz,line; double start, stop; //计时时间 MPI_Status status; double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C,*result_Matrix; double alpha=2,beta=2; // 系数C=aA*B+bC MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); line=M/comm_sz; // 每个进程分多少行数据 Matrix_A=(double*)malloc(M*N*sizeof(double)); Matrix_B=(double*)malloc(N*K*sizeof(double)); Matrix_C=(double*)malloc(M*K*sizeof(double)); buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的数据 buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的数据 ans=(double*)malloc(line*K*sizeof(double)); // 保存部分数据计算结果 result_Matrix=(double*)malloc(M*K*sizeof(double)); // 保存数据计算结果 // 给矩阵A B,C赋值 if(my_rank==0){ start=MPI_Wtime(); for(int i=0;i<M;i++){ for(int j=0;j<N;j++) Matrix_A[i*N+j]=i+1; for(int p=0;p<K;p++) Matrix_C[i*K+p]=1; } for(int i=0;i<N;i++){ for(int j=0;j<K;j++) Matrix_B[i*K+j]=j+1; } // 输出A,B,C //Matrix_print(Matrix_A,M,N); //Matrix_print(Matrix_B,N,K); //Matrix_print(Matrix_C,M,K); } // 数据分发 MPI_Scatter(Matrix_A,line*N,MPI_DOUBLE,buffer_A,line*N,MPI_DOUBLE,0,MPI_COMM_WORLD); MPI_Scatter(Matrix_C,line*K,MPI_DOUBLE,buffer_C,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 数据广播 MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 计算 结果 for(int i=0;i<line;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++) temp+=buffer_A[i*N+p]*Matrix_B[p*K+j]; ans[i*K+j]=alpha*temp+beta*buffer_C[i*K+j]; } } // 结果聚集 MPI_Gather(ans,line*K,MPI_DOUBLE,result_Matrix,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 计算A剩下的行数据 if(my_rank==0){ int rest=M%comm_sz; if(rest!=0){ for(int i=M-rest-1;i<M;i++) for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++) temp+=Matrix_A[i*N+p]*Matrix_B[p*K+j]; result_Matrix[i*K+j]=alpha*temp+beta*Matrix_C[i*K+j]; } } //Matrix_print(result_Matrix,M,K); stop=MPI_Wtime(); printf("rank:%d time:%lfs\n",my_rank,stop-start); } free(Matrix_A); free(Matrix_B); free(Matrix_C); free(ans); free(buffer_A); free(buffer_C); free(result_Matrix); MPI_Finalize(); return 0; }
三、结果分析
下图为上面两种方法的耗时:并行时,随着进程数目的增多,并行计算的时间越来越短;当达到一定的进程数时,执行时间小到最小值;然后再随着进程数的增多,执行时间反而越来越长。
2、加速比分析:
随着进程数的增大,加速比也是逐渐增大到最大值;再随着进程数的增大,加速比逐渐减小。
3、执行效率分析:
随着进程数的增大,程序执行效率不断降低
由于消息传递需要成本,而且不是每个进程都同时开始和结束,所以随着进程数的上升,平均每进程的效率下降
四、头文件functions.h内容
/********** 输出函数 **********/ void Matrix_print(double *A,int M,int N) { for(int i=0;i<M;i++){ for(int j=0;j<N;j++) printf("%.1f ",A[i*N+j]); printf("\n"); } printf("\n"); }
五、对以上两种程序的不同写法
1、MPI_Send和MPI_Recv
#include <stdio.h> #include "mpi.h" #include <stdlib.h> #include "dgemm_1.h" #define M 1200 // 矩阵维度 #define N 1000 #define K 1300 int main(int argc, char **argv) { int my_rank,comm_sz,line; double start, stop; //计时时间 double alpha=2,beta=2; // 系数C=aA*B+bC MPI_Status status; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); line=M/comm_sz; // 每个进程分多少行数据 // 给矩阵A B,C赋值 if(my_rank==0){ double *Matrix_A,*Matrix_B,*Matrix_C,*result_Matrix; Matrix_A=(double*)malloc(M*N*sizeof(double)); Matrix_B=(double*)malloc(N*K*sizeof(double)); Matrix_C=(double*)malloc(M*K*sizeof(double)); result_Matrix=(double*)malloc(M*K*sizeof(double)); // 保存数据计算结果 start=MPI_Wtime(); for(int i=0;i<M;i++){ for(int j=0;j<N;j++) Matrix_A[i*N+j]=i+1; } for(int i=0;i<N;i++){ for(int j=0;j<K;j++) Matrix_B[i*K+j]=j+1; } for(int i=0;i<M;i++){ for(int j=0;j<K;j++) Matrix_C[i*K+j]=1; } // 输出A,B,C //Matrix_print(Matrix_A,M,N); //Matrix_print(Matrix_B,N,K); //Matrix_print(Matrix_C,M,K); /*将矩阵广播出去*/ for(int i=1;i<comm_sz;i++){ MPI_Send(Matrix_A+(i-1)*line*N,line*N,MPI_DOUBLE,i,66,MPI_COMM_WORLD); MPI_Send(Matrix_C+(i-1)*line*K,line*K,MPI_DOUBLE,i,99,MPI_COMM_WORLD); } MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 接收从进程的计算结果 for(int i=0;i<comm_sz-1;i++){ MPI_Recv(result_Matrix+i*line*K,line*K,MPI_DOUBLE,i+1,33,MPI_COMM_WORLD,&status); } // 计算A剩下的行数据 //matMulti(Matrix_A+(comm_sz-1)*line,Matrix_B,result_Matrix,M,N,K,alpha,beta); for(int i=(comm_sz-1)*line;i<M;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++) temp+=Matrix_A[i*N+p]*Matrix_B[p*K+j]; result_Matrix[i*K+j]=alpha*temp+beta*Matrix_C[i*K+j]; } } //Matrix_print(result_Matrix,M,K); stop=MPI_Wtime(); printf("rank:%d time:%lfs\n",my_rank,stop-start); free(Matrix_A); free(Matrix_B); free(Matrix_C); free(result_Matrix); } else{ double *Matrix_B,*buffer_A,*buffer_C,*ans; Matrix_B=(double*)malloc(N*K*sizeof(double)); buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的数据 buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的数据 //接收广播的数据 MPI_Recv(buffer_A,line*N,MPI_DOUBLE,0,66,MPI_COMM_WORLD,&status); MPI_Recv(buffer_C,line*K,MPI_DOUBLE,0,99,MPI_COMM_WORLD,&status); MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); //计算乘积结果,并将结果发送给主进程 //matMulti(buffer_A,Matrix_B,buffer_C,line,N,K,alpha,beta); for(int i=0;i<line;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++){ temp+=buffer_A[i*N+p]*Matrix_B[p*K+j]; } buffer_C[i*line+j]=alpha*temp+beta*buffer_C[i*K+j]; } } MPI_Send(buffer_C,line*K,MPI_DOUBLE,0,33,MPI_COMM_WORLD); free(Matrix_B); free(buffer_A); free(buffer_C); } MPI_Finalize(); return 0; }
2、MPI_Scatter和MPI_Gather
在上面的MPI_Scatter和MPI_Gather的程序中,矩阵的行数据并不能平均的分配给每个子进程,当不能矩阵的行数据不能与总进程数整除时,就需要单独计算一下矩阵剩余的行数据。if(M%comm_sz!=0){ M-=M%comm_sz; M+=comm_sz; }
#include <stdio.h> #include "mpi.h" #include <stdlib.h> #include "dgemm_1.h" //#define M 3 // 矩阵维度 #define N 2 #define K 4 int main(int argc, char **argv) { int M=3; int my_rank,comm_sz,line; double start, stop; //计时时间 MPI_Status status; double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C,*result_Matrix; double alpha=2,beta=2; // 系数C=aA*B+bC MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); int saveM=M; if(M%comm_sz!=0){ M-=M%comm_sz; M+=comm_sz; } line=M/comm_sz; // 每个进程分多少行数据 Matrix_A=(double*)malloc(M*N*sizeof(double)); Matrix_B=(double*)malloc(N*K*sizeof(double)); Matrix_C=(double*)malloc(M*K*sizeof(double)); buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的数据 buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的数据 ans=(double*)malloc(line*K*sizeof(double)); // 保存部分数据计算结果 result_Matrix=(double*)malloc(M*K*sizeof(double)); // 保存数据计算结果 // 给矩阵A B,C赋值 if(my_rank==0){ for(int i=0;i<M;i++){ for(int j=0;j<N;j++) if(i<saveM) Matrix_A[i*N+j]=i+1; else Matrix_A[i*N+j]=0; for(int p=0;p<K;p++) Matrix_C[i*K+p]=1; } for(int i=0;i<N;i++){ for(int j=0;j<K;j++) Matrix_B[i*K+j]=j+1; } start=MPI_Wtime(); // 输出A,B,C Matrix_print(Matrix_A,saveM,N); Matrix_print(Matrix_B,N,K); Matrix_print(Matrix_C,saveM,K); } // 数据分发 MPI_Scatter(Matrix_A,line*N,MPI_DOUBLE,buffer_A,line*N,MPI_DOUBLE,0,MPI_COMM_WORLD); MPI_Scatter(Matrix_C,line*K,MPI_DOUBLE,buffer_C,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 数据广播 MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD); // 计算结果 for(int i=0;i<line;i++){ for(int j=0;j<K;j++){ double temp=0; for(int p=0;p<N;p++) temp+=buffer_A[i*N+p]*Matrix_B[p*K+j]; ans[i*K+j]=alpha*temp+beta*buffer_C[i*K+j]; } } // 结果聚集 MPI_Gather(ans,line*K,MPI_DOUBLE,result_Matrix,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD); if(my_rank==0){ Matrix_print(result_Matrix,saveM,K); stop=MPI_Wtime(); double E=(double)(5.957/(stop-start))/comm_sz; printf("time:%lfs 并行效率:%.2lf%\n",stop-start,100*E); } free(Matrix_A); free(Matrix_B); free(Matrix_C); free(ans); free(buffer_A); free(buffer_C); free(result_Matrix); MPI_Finalize(); return 0; }
结束。