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位。
二、对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