📂mpi
🔖MPI学习笔记
2022-07-25 20:27阅读: 1948评论: 0推荐: 1

MPI学习笔记（二）：矩阵相乘的两种实现方法

mpi矩阵乘法（C=αAB+βC）

        最近领导让把之前安装的软件lapack、blas里的dgemm运算提取出来独立作为一套程序，然后把这段程序改为并行的，并测试一下进程规模扩展到128时的并行效率。
        我发现这个是dgemm.f文件，里面主要是对C=αAB+βC的实现，因此在此总结一下MPI的矩阵乘法使用。
        其主要思想：是把相乘的矩阵按行分解（任务分解），分别分给不同的进程，然后在汇总到一个进程上，在程序上实现则用到了主从模式，人为的把进程分为主进程和从进程，主进程负责对原始矩阵初始化赋值，并把数据均匀分发（为了负载均衡）到从进程上进行相乘运算，主要用到的知识是MPI点对点通信和组通信的机制。

一、使用简单的MPI_Send和MPI_Recv实现

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#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实现

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#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;
}

三、结果分析

下图为上面两种方法的耗时：

1、 执行时间分析：
并行时，随着进程数目的增多，并行计算的时间越来越短；当达到一定的进程数时，执行时间小到最小值；然后再随着进程数的增多，执行时间反而越来越长。
2、加速比分析：
随着进程数的增大，加速比也是逐渐增大到最大值；再随着进程数的增大，加速比逐渐减小。
3、执行效率分析：
随着进程数的增大，程序执行效率不断降低

由于消息传递需要成本，而且不是每个进程都同时开始和结束，所以随着进程数的上升，平均每进程的效率下降

四、头文件functions.h内容

?
1
2
3
4
5
6
7
8
9
10
/********** 输出函数 **********/
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

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#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的程序中，矩阵的行数据并不能平均的分配给每个子进程，当不能矩阵的行数据不能与总进程数整除时，就需要单独计算一下矩阵剩余的行数据。
?
1
2
3
4
if(M%comm_sz!=0){
    M-=M%comm_sz;
    M+=comm_sz;
}
　　为了使矩阵的行数据能平均的分配给每个子进程，因此记comm_sz为进程数，矩阵的维度需要是comm_sz的倍数。我们将矩阵的维度扩展到comm_sz的倍数，多余的部分用0填充，保证正确性。
?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#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;
}

　　

 

 

 

结束。

 

本文作者：惊小呆

本文链接：https://www.cnblogs.com/babyclass/p/16518726.html

版权声明：本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   惊小呆  阅读(1948)  评论(0)  编辑  收藏  举报
点击右上角即可分享
微信分享提示
💬
评论
📌
收藏
💗
关注
👍
推荐
🚀
回顶
➖
收起
  1. 1 404 not found REOL
404 not found - REOL
00:00 / 00:00
An audio error has occurred.