标签：malloc Matrix int double 矩阵 MPI 相乘 sizeof
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_Marix);
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");
}
结束。
标签：malloc,Matrix,int,double,矩阵,MPI,相乘,sizeof 来源： https://www.cnblogs.com/babyclass/p/16518726.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。