5
MPI学习笔记（二）：矩阵相乘的两种实现方法 - 惊小呆520
source link: https://www.cnblogs.com/babyclass/p/16518726.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
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 900int 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+bCMPI_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 1100int 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+bCMPI_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");} |
本文作者：惊小呆520
本文链接：https://www.cnblogs.com/babyclass/p/16518726.html
版权声明：本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK