mpi矩阵乘法:C=αAB+βC
一、Cannon卡农算法基本介绍
1、二维矩阵串行乘法
两个n维方阵的乘法A×B=C的串行计算公式为:
![](/i/l/?n=22&i=blog/2173981/202208/2173981-20220828201544142-1409089956.png)
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位。
二、对Cannon卡农算法初步探索(主从模式)
1、先实现一下各个进程的数据从主进程发送,(矩阵可以不是方阵、进程总数不是非要开平方)。#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实现
#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实现
#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 、笛卡儿虚拟拓扑
#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头文件
/***** 处理器逻辑阵列坐标到rank的转换 *****/ int get_index(int row,int col,int N){ return ((row+N)%N)*N+(col+N)%N; }
/********** 输出函数 **********/ 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"); }
标签:Matrix,int,double,MPI,Cannon,line,卡农,block From: https://www.cnblogs.com/babyclass/p/16633535.html