📂mpi
🔖MPI学习笔记
2022-09-05 08:59阅读: 2079评论: 0推荐: 0

MPI学习笔记（四）：矩阵相乘的Cannon卡农算法

mpi矩阵乘法：C=αAB+βC

一、Cannon卡农算法基本介绍

1、二维矩阵串行乘法

两个n维方阵的乘法A×B=C的串行计算公式为：

下图是用图示来表示的这种计算规则：

2、二维块划分的思想

并行化二维矩阵乘法的一种思想是二维块划分方法，将p个进程（p为完全平方数）排列成sqrt(p)×sqrt(p)二维网格，然后将矩阵A、B、C都分成sqrt(p)×sqrt(p)块，均匀分布在网格上，矩阵第（i，j）个子块分别记为Aij、Bij、Cij，储在二维进程网格上坐标为（i，j）的进程Pij上。计算Cij时要将Aij（第i行进程上的所有A的子块）和Bij（第j列进程上的所有B的子块）都收集到进程Pij上，再计算Cij，公式可以表示为：

下图是用图示来表示的这种计算规则：

然而每一个进程都重复收集Aik和Bkj会造成空间的大量浪费，时间效率也比较低，于是有了更高效的Cannon算法，其表达式为：

 3、Cannon算法基本思想

每一个进程只存储A、B、C矩阵的一个子块，本地相对应的A、B子块相乘后将结果累加到本地C子块上，然后再与其他进程交换A、B子块，继续相乘将结果累加到本地C子块上。但是要保证对应的A、B子块相乘，不能错位，我们注意到，在最开始，Pij上的A为所在行的第j个，B为所在列的第i个，A和B子块并不对应，因此将一行上的A子块循环向左移动i格，一列上的B子块循环向上移动j格，这样A和B子块就能对应了，以后每执行一次计算，每行所有A子块循环向左移动1格，每列上的B子块循环向上移动1格，A、B子块相乘的结果累加在本地的C子块上。

4、Cannon算法原理

将矩阵A和B分成p个方块Aij和Bij(0<=i,j<=sqrt(p)-1),每块大小为(n/sqrt(p))*(n/sqrt(p)),并将他们分配给sqrt(p)*sqrt(p)个处理器(P00,P01,...,P(sqrt(p)-1,sqrt(p)-1))。开始时处理器Pij存放有块Aij和Bij，并负责计算Cij，然后算法开始执行：
① 将块Aij(0<=i,j<=sqrt(p)-1)向左循环移动i步；
     将块Bij(0<=i,j<=sqrt(p)-1)向上循环移动j步；
② Pij执行乘 - 加运算；
    然后，将块Aij(0<=i,j<=sqrt(p)-1)向左循环移动1步；
    将块Bij(0<=i,j<=sqrt(p)-1)向上循环移动1步；
③ 重复第②步，在Pij中共执行sqrt(p)次块Aij和Bij的循环单位移步。
5、算法举例
下图示例了在16个处理器上，用Cannon算法执行A(4x4)*B(4x4)的过程。其中(a)和(b)对应于上述算法的第①步；(c)、(d)、(e)和(f)对应于上述算法的第②和第③步。注意在算法第①步时，A矩阵的第0行不移位，第1行循环左移1位，第2行循环左移2位。第3行循环左移3位；类似的，B矩阵的第0列不移位，第1列循环上移1位，第2列循环上移2位。第3列循环上移3位。

其实就是在计算Cij

二、对Cannon卡农算法初步探索（主从模式）

 1、先实现一下各个进程的数据从主进程发送，（矩阵可以不是方阵、进程总数不是非要开平方）。
?
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include <math.h>
#include "dgemm_1.h"
 
int main(int argc, char **argv)
{
   int M=4,N=4,K=4; // 矩阵维度
   int my_rank,comm_sz;
   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);
 
   int a=(int)sqrt(comm_sz); // A B行列分多少块
   if(comm_sz!=a*a || a>M || a>N || a>K){
      if(my_rank==0)
         printf("error process:%d\n",comm_sz);
      MPI_Finalize();
      return 0;
   }
 
   int saveM=M,saveN=N,saveK=K; // 为了A B能均分成块
   if(M%a!=0){
      M-=M%a;
      M+=a;
   }
   if(N%a!=0){
      N-=N%a;
      N+=a;
   }
   if(K%a!=0){
      K-=K%a;
      K+=a;
   }
   int each_M=M/a,each_N=N/a,each_K=K/a; // 矩阵A B每块分多少行列数据
 
   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)); // 保存数据计算结果
 
      // 给矩阵A B,C赋值
      init_Matrix(Matrix_A,Matrix_B,Matrix_C,M,N,K,saveM,saveN,saveK);
      // 输出A,B,C
      //Matrix_print2(Matrix_A,saveM,saveN,N);
      //Matrix_print2(Matrix_B,saveN,saveK,K);
      //Matrix_print2(Matrix_C,saveM,saveK,K);
      printf("a=%d each_M=%d each_N=%d each_K=%d\n",a,each_M,each_N,each_K);
 
      start=MPI_Wtime();
      // 主进程计算第1块
      for(int i=0;i<each_M;i++){
         for(int j=0;j<each_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];
         }
      }
 
      // 向其它进程发送块数据
      for(int i=1;i<comm_sz;i++){
         int beginRow=(i/a)*each_M; // 每个块的行列起始位置(坐标/偏移量)
         int beginCol=(i%a)*each_K;
         for(int j=0;j<each_M;j++)
            MPI_Send(Matrix_C+(beginRow+j)*K+beginCol,each_K,MPI_DOUBLE,i,j+each_M+each_N,MPI_COMM_WORLD);
         // 发送A B每块数据
         for(int k=0;k<a;k++){
            int begin_part=k*each_N; // 移动A的列 B的行 即A列不同程度的左移，B行不同程度的上移
            for(int j=0;j<each_M;j++)
               MPI_Send(Matrix_A+(beginRow+j)*N+begin_part,each_N,MPI_DOUBLE,i,j,MPI_COMM_WORLD);
            for(int p=0;p<each_N;p++)
               MPI_Send(Matrix_B+(begin_part+p)*K+beginCol,each_K,MPI_DOUBLE,i,p+each_M,MPI_COMM_WORLD);
         }
      }
      // 接收从进程的计算结果
      for(int i=1;i<comm_sz;i++){
         int beginRow=(i/a)*each_M;
         int endRow=beginRow+each_M;
         int beginCol=(i%a)*each_K;
         for(int j=beginRow;j<endRow;j++)
            MPI_Recv(result_Matrix+j*K+beginCol,each_K,MPI_DOUBLE,i,j-beginRow+2*each_M+each_N,MPI_COMM_WORLD,&status);
      }
 
      Matrix_print2(result_Matrix,saveM,saveK,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 *buffer_A,*buffer_B,*buffer_C;
      buffer_A=(double*)malloc(each_M*each_N*sizeof(double)); // A的均分行的数据
      buffer_B=(double*)malloc(each_N*each_K*sizeof(double)); // B的均分列的数据
      buffer_C=(double*)malloc(each_M*each_K*sizeof(double)); // C的均分行的数据
 
      // 接收C块数据
      for(int j=0;j<each_M;j++)
         MPI_Recv(buffer_C+j*each_K,each_K,MPI_DOUBLE,0,j+each_M+each_N,MPI_COMM_WORLD,&status);
 
      for(int k=0;k<a;k++){ // 把每块数据求和
         //接收A B块数据
         for(int j=0;j<each_M;j++)
            MPI_Recv(buffer_A+j*each_N,each_N,MPI_DOUBLE,0,j,MPI_COMM_WORLD,&status);
         for(int p=0;p<each_N;p++)
            MPI_Recv(buffer_B+p*each_K,each_K,MPI_DOUBLE,0,p+each_M,MPI_COMM_WORLD,&status);
 
         //计算乘积结果,并将结果发送给主进程
         for(int i=0;i<each_M;i++){
            for(int j=0;j<each_K;j++){
               double temp=0;
               for(int p=0;p<each_N;p++){
                  temp+=buffer_A[i*each_N+p]*buffer_B[p*each_K+j];
               }
               if(k==0)
                  buffer_C[i*each_K+j]=alpha*temp+beta*buffer_C[i*each_K+j];
               else
                  buffer_C[i*each_K+j]+=alpha*temp;
            }
         }
         //matMulti(buffer_A,buffer_B,buffer_C,each_row,N,each_col,alpha,beta);
      }
      // 将结果发送给主进程
      for(int j=0;j<each_M;j++){
         MPI_Send(buffer_C+j*each_K,each_K,MPI_DOUBLE,0,j+2*each_M+each_N,MPI_COMM_WORLD);
      }
 
      free(buffer_A);
      free(buffer_B);
      free(buffer_C);
   }
   MPI_Finalize();
   return 0;
}

三、对等模式的Cannon卡农算法

 1、使用MPI_Isend实现

?
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include <math.h>
#include "dgemm_1.h"
/****** MPI_Isend ******/
int main(int argc, char **argv)
{
   int N=10000; // 矩阵维度
   int my_rank,comm_sz,moveRow,moveCol; // 每个块移动位置
   double start, stop; //计时时间
   double alpha=2,beta=2; // 系数C=aA*B+bC
   MPI_Status status;
   MPI_Request reqs;
 
   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
 
   int block=(int)sqrt(comm_sz); // A B行列分多少块
   if(comm_sz!=block*block || block>N ){
      if(my_rank==0)
         printf("error process:%d\n",comm_sz);
      MPI_Finalize();
      return 0;
   }
 
   int saveN=N; // 为了A B能均分成块
   if(N%block!=0){
      N-=N%block;
      N+=block;
   }
   int line=N/block; // 矩阵A B每块分多少行列数据
   double *buffer_A,*buffer_B,*buffer_C,*ans;
   buffer_A=(double*)malloc(line*line*sizeof(double)); // A的块数据
   buffer_B=(double*)malloc(line*line*sizeof(double)); // B的块数据
   buffer_C=(double*)malloc(line*line*sizeof(double)); // C的块数据
   ans=(double*)malloc(line*line*sizeof(double)); // 临时存储每块的结果数据
 
   if(my_rank==0){
      double *Matrix_A,*Matrix_B,*Matrix_C;
      Matrix_A=(double*)malloc(N*N*sizeof(double));
      Matrix_B=(double*)malloc(N*N*sizeof(double));
      Matrix_C=(double*)malloc(N*N*sizeof(double));
 
      // 给矩阵A B,C赋值
      for(int i=0;i<N;i++){
         for(int j=0;j<N;j++){
            if(i<saveN && j<saveN){
               Matrix_A[i*N+j]=j-i;
               Matrix_B[i*N+j]=i+j;
               Matrix_C[i*N+j]=i*j;
            }
            else{
               Matrix_A[i*N+j]=0;
               Matrix_B[i*N+j]=0;
               Matrix_C[i*N+j]=0;
            }
         }
      }
      //init_Matrix(Matrix_A,Matrix_B,Matrix_C,N,N,N,saveN,saveN,saveN);  
      // 输出A,B,C
      //Matrix_print2(Matrix_A,saveN,saveN,N);
      //Matrix_print2(Matrix_B,saveN,saveN,N);
      //Matrix_print2(Matrix_C,saveN,saveN,N);
      printf("block=%d line=%d\n",block,line);
 
      start=MPI_Wtime();
      // 将每块数据分配给每个进程
      for(int p=0;p<comm_sz;p++){
         int beginRow=(p/block)*line; // 每个块的行列起始位置(坐标/偏移量)
         int beginCol=(p%block)*line;
         if(p==0){
            for(int i=0;i<line;i++)
               for(int j=0;j<line;j++){
                  buffer_A[i*line+j]=Matrix_A[i*N+j];
                  buffer_B[i*line+j]=Matrix_B[i*N+j];
                  buffer_C[i*line+j]=Matrix_C[i*N+j];
               }
         }
         else{
            if((p-p/block)<(p/block)*block)
               moveRow=p-p/block+block;
            else
               moveRow=p-p/block;
            if((p-(p%block)*block)<0)
               moveCol=p-(p%block)*block+comm_sz;
            else
               moveCol=p-(p%block)*block;
            for(int i=0;i<line;i++){
               MPI_Send(Matrix_A+(beginRow+i)*N+beginCol,line,MPI_DOUBLE,moveRow,i,MPI_COMM_WORLD);
               MPI_Send(Matrix_B+(beginRow+i)*N+beginCol,line,MPI_DOUBLE,moveCol,i+line,MPI_COMM_WORLD);
               MPI_Send(Matrix_C+(beginRow+i)*N+beginCol,line,MPI_DOUBLE,p,i+2*line,MPI_COMM_WORLD);
            }
         }
      }
      free(Matrix_A);
      free(Matrix_B);
      free(Matrix_C);
   }
   else{
      // 接收A B C块数据
      for(int i=0;i<line;i++){
         MPI_Recv(buffer_A+i*line,line,MPI_DOUBLE,0,i,MPI_COMM_WORLD,&status);
         MPI_Recv(buffer_B+i*line,line,MPI_DOUBLE,0,i+line,MPI_COMM_WORLD,&status);
         MPI_Recv(buffer_C+i*line,line,MPI_DOUBLE,0,i+2*line,MPI_COMM_WORLD,&status);
      }
   }
   //计算第一次乘积结果
   for(int i=0;i<line;i++){
      for(int j=0;j<line;j++){
         double temp=0;
         for(int p=0;p<line;p++){
            temp+=buffer_A[i*line+p]*buffer_B[p*line+j];
         }
         ans[i*line+j]=alpha*temp+beta*buffer_C[i*line+j];
      }
   }
   // 移动块数据
   if(my_rank==(my_rank/block)*block)
      moveRow=my_rank-1+block;
   else
      moveRow=my_rank-1;
   if((my_rank-block)<0)
      moveCol=my_rank-block+comm_sz;
   else
      moveCol=my_rank-block;
   MPI_Isend(buffer_A,line*line,MPI_DOUBLE,moveRow,3*line,MPI_COMM_WORLD,&reqs);
   MPI_Isend(buffer_B,line*line,MPI_DOUBLE,moveCol,4*line,MPI_COMM_WORLD,&reqs);
   //MPI_Wait(&reqs,&status); //阻塞发送矩阵到结束
   for(int k=1;k<block;k++){ // 把每块数据求和
      //接收A B块数据
      if((my_rank+1)==((my_rank/block)+1)*block)
         moveRow=my_rank+1-block;
      else
         moveRow=my_rank+1;
      if((my_rank+block)>=comm_sz)
         moveCol=my_rank+block-comm_sz;
      else
         moveCol=my_rank+block;
      MPI_Recv(buffer_A,line*line,MPI_DOUBLE,moveRow,3*line,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_B,line*line,MPI_DOUBLE,moveCol,4*line,MPI_COMM_WORLD,&status);
      //计算乘积结果,并将结果发送给主进程
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++){
            double temp=0;
            for(int p=0;p<line;p++){
               temp+=buffer_A[i*line+p]*buffer_B[p*line+j];
            }
            ans[i*line+j]+=alpha*temp;
         }
      }
      //matMulti(buffer_A,buffer_B,buffer_C,each_row,N,each_col,alpha,beta);
      // 移动块数据
      if((my_rank-1)<(my_rank/block)*block)
         moveRow=my_rank-1+block;
      else
         moveRow=my_rank-1;
      if((my_rank-block)<0)
         moveCol=my_rank-block+comm_sz;
      else
         moveCol=my_rank-block;
      MPI_Isend(buffer_A,line*line,MPI_DOUBLE,moveRow,3*line,MPI_COMM_WORLD,&reqs);
      MPI_Isend(buffer_B,line*line,MPI_DOUBLE,moveCol,4*line,MPI_COMM_WORLD,&reqs);
      //MPI_Wait(&reqs,&status); //阻塞发送矩阵到结束
   }
   // 将结果发送给主进程
   if(my_rank!=0){
      for(int i=0;i<line;i++)
         MPI_Isend(ans+i*line,line,MPI_DOUBLE,0,i+5*line,MPI_COMM_WORLD,&reqs);
      //MPI_Wait(&reqs,&status); //阻塞发送矩阵到结束
   }
   if(my_rank==0){
      // 接收从进程的计算结果
      double *result_Matrix;
      result_Matrix=(double*)malloc(N*N*sizeof(double)); // 保存数据计算结果
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++)
            result_Matrix[i*N+j]=ans[i*line+j];
      }
      for(int p=1;p<comm_sz;p++){
         int beginRow=(p/block)*line;
         int endRow=beginRow+line;
         int beginCol=(p%block)*line;
         for(int i=beginRow;i<endRow;i++)
            MPI_Recv(result_Matrix+i*N+beginCol,line,MPI_DOUBLE,p,i-beginRow+5*line,MPI_COMM_WORLD,&status);
      }
 
      //Matrix_print2(result_Matrix,saveN,saveN,N);
      stop=MPI_Wtime();
      printf("rank:%d time:%lfs\n",my_rank,stop-start);
      free(result_Matrix);
   }
 
   free(buffer_A);
   free(buffer_B);
   free(buffer_C);
   free(ans);
 
   MPI_Finalize();
   return 0;
}

2、使用MPI_Sendrecv_replace memcpy实现

?
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "dgemm_1.h"
/*** MPI_Sendrecv_replace memcpy:string.h ***/
int main(int argc, char **argv)
{
   int N=10000; // 矩阵维度
   int my_rank,comm_sz;
   double start, stop; //计时时间
   double alpha=2,beta=2; // 系数C=aA*B+bC
   MPI_Status status;
   MPI_Request reqs;
 
   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
 
   int block=(int)sqrt(comm_sz); // A B行列分多少块
   if(comm_sz!=block*block || block>N ){
      if(my_rank==0)
         printf("error process:%d\n",comm_sz);
      MPI_Finalize();
      return 0;
   }
 
   int saveN=N; // 为了A B能均分成块
   if(N%block!=0){
      N-=N%block;
      N+=block;
   }
   int line=N/block; // 矩阵A B每块分多少行列数据
   int my_Row=my_rank/block; // 每个块当前位置(i,j)
   int my_Col=my_rank%block; // 每个块当前位置(i,j)
 
   double *buffer_A,*buffer_B,*buffer_C,*ans;
   buffer_A=(double*)malloc(line*line*sizeof(double)); // A的块数据
   buffer_B=(double*)malloc(line*line*sizeof(double)); // B的块数据
   buffer_C=(double*)malloc(line*line*sizeof(double)); // C的块数据
   ans=(double*)malloc(line*line*sizeof(double)); // 临时存储每块的结果数据
 
   if(my_rank==0){
      double *Matrix_A,*Matrix_B,*Matrix_C;
      Matrix_A=(double*)malloc(N*N*sizeof(double));
      Matrix_B=(double*)malloc(N*N*sizeof(double));
      Matrix_C=(double*)malloc(N*N*sizeof(double));
 
      // 给矩阵A B,C赋值
      for(int i=0;i<N;i++){
         for(int j=0;j<N;j++){
            if(i<saveN && j<saveN){
               Matrix_A[i*N+j]=j-i;
               Matrix_B[i*N+j]=i+j;
               Matrix_C[i*N+j]=i*j;
            }
            else{
               Matrix_A[i*N+j]=0;
               Matrix_B[i*N+j]=0;
               Matrix_C[i*N+j]=0;
            }
         }
      }
      //init_Matrix(Matrix_A,Matrix_B,Matrix_C,N,N,N,saveN,saveN,saveN);  
      // 输出A,B,C
      //Matrix_print2(Matrix_A,saveN,saveN,N);
      //Matrix_print2(Matrix_B,saveN,saveN,N);
      //Matrix_print2(Matrix_C,saveN,saveN,N);
      printf("block=%d line=%d\n",block,line);
 
      start=MPI_Wtime();
      // 将每块数据分配给每个进程
      for(int p=1;p<comm_sz;p++){
         int beginRow=(p/block)*line; // 每个块的行列起始位置(坐标/偏移量)
         int beginCol=(p%block)*line;
         for(int i=0;i<line;i++){
            memcpy(buffer_A+i*line,Matrix_A+(beginRow+i)*N+beginCol,line*sizeof(double));  // memory copy内存复制
            memcpy(buffer_B+i*line,Matrix_B+(beginRow+i)*N+beginCol,line*sizeof(double));
            memcpy(buffer_C+i*line,Matrix_C+(beginRow+i)*N+beginCol,line*sizeof(double));
            //for(int j=0;j<line;j++){
               //buffer_A[i*line+j]=Matrix_A[(beginRow+i)*N+beginCol+j];
               //buffer_B[i*line+j]=Matrix_B[(beginRow+i)*N+beginCol+j];
               //buffer_C[i*line+j]=Matrix_C[(beginRow+i)*N+beginCol+j];
            //}
         }
         MPI_Send(buffer_A,line*line,MPI_DOUBLE,p,0,MPI_COMM_WORLD);
         MPI_Send(buffer_B,line*line,MPI_DOUBLE,p,1,MPI_COMM_WORLD);
         MPI_Send(buffer_C,line*line,MPI_DOUBLE,p,2,MPI_COMM_WORLD);
      }
      // id为0的处理器直接拷贝过去
      for(int i=0;i<line;i++){
         memcpy(buffer_A+i*line,Matrix_A+i*N,line*sizeof(double));  // memory copy内存复制
         memcpy(buffer_B+i*line,Matrix_B+i*N,line*sizeof(double));
         memcpy(buffer_C+i*line,Matrix_C+i*N,line*sizeof(double));
         //for(int j=0;j<line;j++){
            //buffer_A[i*line+j]=Matrix_A[i*N+j];
            //buffer_B[i*line+j]=Matrix_B[i*N+j];
            //buffer_C[i*line+j]=Matrix_C[i*N+j];
         //}
      }
 
      free(Matrix_A);
      free(Matrix_B);
      free(Matrix_C);
   }
   else{
      // 接收A B C块数据
      MPI_Recv(buffer_A,line*line,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_B,line*line,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_C,line*line,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_C,line*line,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&status);
   }
   /*
     将A中坐标为(i, j)的分块循环左移i步
     将B中坐标为(i, j)的分块循环上移j步
   */
   //MPI_Sendrecv(buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col-my_Row,block),3,buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col+my_Row,block),3,MPI_COMM_WORLD,&status);
   //MPI_Sendrecv(buffer_B,line*line,MPI_DOUBLE,get_index(my_Row-my_Col,my_Col,block),4,buffer_B,line*line,MPI_DOUBLE,get_index(my_Row+my_Col,my_Col,block),4,MPI_COMM_WORLD,&status);
   MPI_Sendrecv_replace(buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col-my_Row,block),3,get_index(my_Row,my_Col+my_Row,block),3,MPI_COMM_WORLD,&status);
   MPI_Sendrecv_replace(buffer_B,line*line,MPI_DOUBLE,get_index(my_Row-my_Col,my_Col,block),4,get_index(my_Row+my_Col,my_Col,block),4,MPI_COMM_WORLD,&status);
   for(int k=0;k<block;k++){
      // 把每块数据求和
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++){
            double temp=0;
            for(int p=0;p<line;p++){
               temp+=buffer_A[i*line+p]*buffer_B[p*line+j];
            }
            if(k==0)
               ans[i*line+j]=alpha*temp+beta*buffer_C[i*line+j];
            else
               ans[i*line+j]+=alpha*temp;
         }
      }
      // 矩阵A每行左移一位，矩阵B每行上移一位
      //MPI_Sendrecv(buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col-1,block),5,buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col+1,block),5,MPI_COMM_WORLD,&status);
      //MPI_Sendrecv(buffer_B,line*line,MPI_DOUBLE,get_index(my_Row-1,my_Col,block),6,buffer_B,line*line,MPI_DOUBLE,get_index(my_Row+1,my_Col,block),6,MPI_COMM_WORLD,&status);
      MPI_Sendrecv_replace(buffer_A,line*line,MPI_DOUBLE,get_index(my_Row,my_Col-1,block),5,get_index(my_Row,my_Col+1,block),5,MPI_COMM_WORLD,&status);
      MPI_Sendrecv_replace(buffer_B,line*line,MPI_DOUBLE,get_index(my_Row-1,my_Col,block),6,get_index(my_Row+1,my_Col,block),6,MPI_COMM_WORLD,&status);
   }
   // 将结果发送给主进程
   if(my_rank!=0){
      MPI_Send(ans,line*line,MPI_DOUBLE,0,7,MPI_COMM_WORLD);
   }
   else{
      // 接收从进程的计算结果
      double *result_Matrix;
      result_Matrix=(double*)malloc(N*N*sizeof(double)); // 保存数据计算结果
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++)
            result_Matrix[i*N+j]=ans[i*line+j];
      }
      for(int p=1;p<comm_sz;p++){
         int beginRow=(p/block)*line;
         int beginCol=(p%block)*line;
         MPI_Recv(ans,line*line,MPI_DOUBLE,p,7,MPI_COMM_WORLD,&status);
         for(int i=0;i<line;i++){
            memcpy(result_Matrix+(beginRow+i)*N+beginCol,ans+i*line,line*sizeof(double));  // memory copy内存复制
            //for(int j=0;j<line;j++)
               //result_Matrix[(beginRow+i)*N+beginCol+j]=ans[i*line+j];
         }
      }
 
      //Matrix_print2(result_Matrix,saveN,saveN,N);
      stop=MPI_Wtime();
      printf("rank:%d time:%lfs\n",my_rank,stop-start);
      free(result_Matrix);
   }
 
   free(buffer_A);
   free(buffer_B);
   free(buffer_C);
   free(ans);
 
   MPI_Finalize();
   return 0;
}

3、使用MPI_Sendrecv_replace、向量派生数据:MPI_Type_vector 、笛卡儿虚拟拓扑

?
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "dgemm_1.h"
/*** MPI_Sendrecv_replace 向量派生数据:MPI_Type_vector 笛卡儿虚拟拓扑 ***/
int main(int argc, char **argv)
{
   int N=10000; // 矩阵维度
   int my_rank,comm_sz;
   double start, stop; //计时时间
   double alpha=2,beta=2; // 系数C=aA*B+bC
   MPI_Status status;
   MPI_Request reqs;
 
   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
 
   int block=(int)sqrt(comm_sz); // A B行列分多少块
   if(comm_sz!=block*block || block>N ){
      if(my_rank==0)
         printf("error process:%d\n",comm_sz);
      MPI_Finalize();
      return 0;
   }
 
   int saveN=N; // 为了A B能均分成块
   if(N%block!=0){
      N-=N%block;
      N+=block;
   }
   int line=N/block; // 矩阵A B每块分多少行列数据
   int my_Row=my_rank/block; // 每个块当前位置(i,j)
   int my_Col=my_rank%block; // 每个块当前位置(i,j)
 
   double *buffer_A,*buffer_B,*buffer_C,*ans;
   buffer_A=(double*)malloc(line*line*sizeof(double)); // A的块数据
   buffer_B=(double*)malloc(line*line*sizeof(double)); // B的块数据
   buffer_C=(double*)malloc(line*line*sizeof(double)); // C的块数据
   ans=(double*)malloc(line*line*sizeof(double)); // 临时存储每块的结果数据
 
   if(my_rank==0){
      double *Matrix_A,*Matrix_B,*Matrix_C;
      Matrix_A=(double*)malloc(N*N*sizeof(double));
      Matrix_B=(double*)malloc(N*N*sizeof(double));
      Matrix_C=(double*)malloc(N*N*sizeof(double));
 
      // 给矩阵A B,C赋值
      for(int i=0;i<N;i++){
         for(int j=0;j<N;j++){
            if(i<saveN && j<saveN){
               Matrix_A[i*N+j]=j-i;
               Matrix_B[i*N+j]=i+j;
               Matrix_C[i*N+j]=i*j;
            }
            else{
               Matrix_A[i*N+j]=0;
               Matrix_B[i*N+j]=0;
               Matrix_C[i*N+j]=0;
            }
         }
      }
      //init_Matrix(Matrix_A,Matrix_B,Matrix_C,N,N,N,saveN,saveN,saveN);  
      // 输出A,B,C
      //Matrix_print2(Matrix_A,saveN,saveN,N);
      //Matrix_print2(Matrix_B,saveN,saveN,N);
      //Matrix_print2(Matrix_C,saveN,saveN,N);
      printf("block=%d line=%d\n",block,line);
 
      MPI_Datatype columntype; // 定义新列向量类型
      MPI_Type_vector(line,line,N,MPI_DOUBLE,&columntype);// 块数 块长度 跨度
      MPI_Type_commit(&columntype);
      start=MPI_Wtime();
      // 将每块数据分配给每个进程
      for(int i=0,k=0;i<block;i++,k+=block)
         for (int j=0;j<block;j++){
            if(i==0&&j==0)
               ++j;
            MPI_Send(Matrix_A+(i*line)*N+j*line,1,columntype,k+j,0,MPI_COMM_WORLD);
            MPI_Send(Matrix_B+(i*line)*N+j*line,1,columntype,k+j,1,MPI_COMM_WORLD);
            MPI_Send(Matrix_C+(i*line)*N+j*line,1,columntype,k+j,2,MPI_COMM_WORLD);
         }
 
      // id为0的处理器直接拷贝过去
      for(int i=0;i<line;i++){
         memcpy(buffer_A+i*line,Matrix_A+i*N,line*sizeof(double));  // memory copy内存复制
         memcpy(buffer_B+i*line,Matrix_B+i*N,line*sizeof(double));
         memcpy(buffer_C+i*line,Matrix_C+i*N,line*sizeof(double));
      }
 
      free(Matrix_A);
      free(Matrix_B);
      free(Matrix_C);
   }
   else{
      // 接收A B C块数据
      MPI_Recv(buffer_A,line*line,MPI_DOUBLE,0,0,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_B,line*line,MPI_DOUBLE,0,1,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_C,line*line,MPI_DOUBLE,0,2,MPI_COMM_WORLD,&status);
   }
   // 上下左右进程号 每一维的进程数 一维上网格的周期性 标识数是否可以重排序
   int nbrs[4],dims[2]={block,block},periods[2]={1,1},reorder=0,coords[2];
   MPI_Comm cartcomm;
   // 定义虚拟进程拓扑 它是一个2维 dims:block*block的网格
   MPI_Cart_create(MPI_COMM_WORLD,2,dims,periods,reorder,&cartcomm);
 
   // 将进程rank的顺序编号转换为笛卡儿坐标coords 2:维数
   MPI_Cart_coords(cartcomm,my_rank,2,coords);
   // 得到当前进程上下j 左右i 的进程标识
   MPI_Cart_shift(cartcomm,0,coords[1],&nbrs[0],&nbrs[1]);
   MPI_Cart_shift(cartcomm,1,coords[0],&nbrs[2],&nbrs[3]);
 
   /*
     将A中坐标为(i, j)的分块循环左移i步
     将B中坐标为(i, j)的分块循环上移j步
   */
   MPI_Sendrecv_replace(buffer_A,line*line,MPI_DOUBLE,nbrs[2],3,nbrs[3],3,MPI_COMM_WORLD,&status);
   MPI_Sendrecv_replace(buffer_B,line*line,MPI_DOUBLE,nbrs[0],4,nbrs[1],4,MPI_COMM_WORLD,&status);
 
   // 得到当前进程上下左右的进程标识
   MPI_Cart_shift(cartcomm,0,1,&nbrs[0],&nbrs[1]);
   MPI_Cart_shift(cartcomm,1,1,&nbrs[2],&nbrs[3]);
   for(int k=0;k<block;k++){
      // 把每块数据求和
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++){
            double temp=0;
            for(int p=0;p<line;p++){
               temp+=buffer_A[i*line+p]*buffer_B[p*line+j];
            }
            if(k==0)
               ans[i*line+j]=alpha*temp+beta*buffer_C[i*line+j];
            else
               ans[i*line+j]+=alpha*temp;
         }
      }
      // 矩阵A每行左移一位，矩阵B每行上移一位
      MPI_Sendrecv_replace(buffer_A,line*line,MPI_DOUBLE,nbrs[2],5,nbrs[3],5,MPI_COMM_WORLD,&status);
      MPI_Sendrecv_replace(buffer_B,line*line,MPI_DOUBLE,nbrs[0],6,nbrs[1],6,MPI_COMM_WORLD,&status);
   }
   // 将结果发送给主进程
   if(my_rank!=0){
      MPI_Send(ans,line*line,MPI_DOUBLE,0,7,MPI_COMM_WORLD);
   }
   else{
      // 接收从进程的计算结果
      double *result_Matrix;
      result_Matrix=(double*)malloc(N*N*sizeof(double)); // 保存数据计算结果
      for(int i=0;i<line;i++){
         for(int j=0;j<line;j++)
            result_Matrix[i*N+j]=ans[i*line+j];
      }
      for(int p=1;p<comm_sz;p++){
         int beginRow=(p/block)*line;
         int beginCol=(p%block)*line;
         MPI_Recv(ans,line*line,MPI_DOUBLE,p,7,MPI_COMM_WORLD,&status);
         for(int i=0;i<line;i++){
            memcpy(result_Matrix+(beginRow+i)*N+beginCol,ans+i*line,line*sizeof(double));  // memory copy内存复制
            //for(int j=0;j<line;j++)
               //result_Matrix[(beginRow+i)*N+beginCol+j]=ans[i*line+j];
         }
      }
 
      //Matrix_print2(result_Matrix,saveN,saveN,N);
      stop=MPI_Wtime();
      printf("rank:%d time:%lfs\n",my_rank,stop-start);
      free(result_Matrix);
   }
 
   free(buffer_A);
   free(buffer_B);
   free(buffer_C);
   free(ans);
 
   MPI_Finalize();
   return 0;
}　　

 四、附录

1、dgemm_1.h头文件

?
1
2
3
4
/***** 处理器逻辑阵列坐标到rank的转换 *****/
int get_index(int row,int col,int N){
   return ((row+N)%N)*N+(col+N)%N;
}
?
1
2
3
4
5
6
7
8
9
10
/********** 输出函数 **********/
void Matrix_print2(double *A,int M,int saveN,int N)
{
   for(int i=0;i<M;i++){
      for(int j=0;j<saveN;j++)
         printf("%.1lf ",A[i*N+j]);
      printf("\n");
   }
   printf("\n");
}

　　

 

 

 

 

 

 

 

 

本文作者：惊小呆

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

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

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