PROCEDURE 1
INPUT: A , B . A and B are both in processor 0.
OUTPUT: C = A × B in processor 0.
BEGIN:#数据散发阶段，进程0把矩阵A和B分发给各个进程
if myid == 0:for i in range(1,p):send(Ai, i), send(Bi, i)barrier()
else:recv(Amyid, 0), recv(Bmyid, 0)barrier()#并行计算各个分块矩阵的阶段
left = (myid - 1)%p  , right = (myid + 1)%p
for i in range(0, p):k = (i + myid)%pC[myid][k] = Amyid * Bkif i != p-1 :send(Bk,left), recv(Bk+1,right)#聚集阶段，每个节点把自己的计算结果发送给进程0
if myid == 0:for i in range(1, p):j = (i - 1)%pstart_new_thread(recv, (C[j],i))
else:j = (myid - 1)%psend(C[j],0)
return C
END.

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>#define M 1000
#define N 1000
#define K 1000/* 职责：读取数据到矩阵中 */
void read_to_mat(int mat[][],int m,int n){int i = 0,j = 0;char filename[80];FILE *fp = NULL;fscanf(stdin,"%s",filename);fp = fopen(filename,"r");for(i = 0;i < m;i++){for(j = 0;j < n;j++){if(EOF == fscanf(fp,"%d",&mat[i][j])){printf("[ERROR] - filename: %s",filename);exit(1);}}}fclose(fp);return;
}void write_to_file(int mat[][],int m,int n){char filename[] = "result.txt"FILE *fp = fopen(filename,"w");for(int i = 0;i < m;i++){if(i > 0)fprintf(fp,"\n");for(int j = 0;j < n;j++){fprintf(fp,"%d",mat[i][j]);if(j != n-1)fprintf(fp,"\t");}}fclose(fp);return;
}void scatter(int mat[][],int src,int size,int displs[],int sendcnts[]){int myid = 0;MPI_Comm_rank(MPI_COMM_WORLD,&myid);if(myid == src){for(int dest = 0;dest < size;dest++){if(dest != src){for(int j = displs[dest];j < displs[dest] + sendcnts[dest];j++){MPI_Bsend(mat[j],K,MPI_INT,dest,0,MPI_COMM_WORLD);}}}}else{for(int j = displs[myid];j < displs[myid] + sendcnts[myid];j++){MPI_Recv(mat[j],K,MPI_INT,src,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE);}}return;
}typedef struct{int displs[];int sendcnts[];int mat[][];int src;
} RECV_ARG;void receive(RECV_ARG *arg){for(int j = arg->displs[arg->src];j < arg->displs[arg->src] + arg->sendcnts[arg->src];j++){MPI_Recv(arg->mat[j],K,MPI_INT,arg->src,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE);}return;
}void gather(int mat[][],int dest,int size,int displs[],int sendcnts[]){int myid = 0;pthread_t tid = 0;RECV_ARG * arg = (RECV_ARG *)malloc(sizeof(RECV_ARG));MPI_Comm_rank(MPI_COMM_WORLD,&myid);if(myid == dest){for(int src = 0;src < size;src++){if(src != dest){arg->displs = displs;arg->sendcnts = sendcnts;arg->mat = mat;arg->src = src;pthread_create(&tid,NULL,receive,(RECV_ARG *)&arg);}}}else{for(int j = displs[myid];j < displs[myid] + sendcnts[myid];j++){MPI_Send(mat[j],K,MPI_INT,src,0,MPI_COMM_WORLD);}}free(arg);return;
}int main(int argc, char *argv[]){int myid = 0, numprocs = 0,A[M][K],B[N][K],C[M][N],/* 变量B存储了真实矩阵的转置 */sendcnts_row[],sendcnts_col[],displs_row[],displs_col[],*buf = NULL;pack_size = 0,buf_size = 0,left = 0, right = 0;/* 起始阶段，把两个矩阵分块 */MPI_Init(&argc,&argv);MPI_Comm_rank(MPI_COMM_WORLD,&myid);MPI_Comm_size(MPI_COMM_WORLD,&numprocs);if(numprocs < 2){printf("[ERROR] - numprocs = %d", numprocs);exit(2);}sendcnts_row = (int *)malloc(sizeof(int)*numprocs);sendcnts_col = (int *)malloc(sizeof(int)*numprocs);    for(int i = 0;i < numprocs;i++){if(i == 0){displs_row[i] = 0;displs_col[i] = 0;sendcnts_row[i] = (M/numprocs) + M%numprocs;sendcnts_col[i] = (N/numprocs) + N%numprocs;} else {displs_row[i] = i * (M/numprocs) + M%numprocs;displs_col[i] = i * (N/numprocs) + N%numprocs;sendcnts_row[i] = M/numprocs;sendcnts_col[i] = N/numprocs;}}MPI_Pack_size(M*K+K*N,MPI_INT,&pack_size);buf_size = (M + N) * MPI_BSEND_OVERHEAD + pack_size;buf = (int*)malloc(sizeof(int)*buf_size);/* 数据散发阶段，形成通信器的初始状态 */MPI_Buffer_attach(buf,buf_size);if(myid == 0){read_to_mat(A,M,K);    read_to_mat(B,N,K);/* 把真实矩阵映射为转置矩阵 */}scatter(A,0,numprocs,displs_row,sendcnts_row);scatter(B,0,numprocs,displs_col,ndcnts_col);MPI_Buffer_detach(&buf,&buf_size);/* 计算并循环通信阶段，i的值引用了通信器的每个状态 */MPI_Buffer_attach(buf,buf_size);left = (myid - 1) % numprocs;right = (myid + 1) % numprocs;for(int i = 0;i < numprocs;i++){int p = (i + myid) % numprocs;int left_p = (p - 1) % numprocs;for(int n = displs_col[p];n < displs_col[p] + sendcnts_col[p];n++){for(int m = displs_row[myid];m < displs_row[myid] + sendcnts_row[myid];m++){int sum = 0;for(int k = 0;k < K;k++){sum += A[m][k] + B[n][k];}C[m][n] = sum;}if(i != numprocs - 1)MPI_Bsend(B[n],K,MPI_INT,right,0,MPI_COMM_WORLD);/* 一边计算一边发送 */}if(i != numprocs - 1){for(int col = displs_col[left_p];col < displs_col[left_p] + sendcnts_col[left_p];col++){MPI_Recv(B[col],K,MPI_INT,left,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE);}}}MPI_Buffer_detach(&buf,&buf_size);/* 聚集阶段 */MPI_Buffer_attach(buf,buf_size);gather(C,0,numprocs,displs_row,sendcnts_row);MPI_Buffer_detach(&buf,&buf_size);if(myid ==0)write_to_file(C,M,N);free(buf);MPI_Finalize();return 0;
}/* 实际的代码和伪代码有一些区别，
* 实际代码用B存储实际矩阵的转置，
* 实际代码采用一边计算一边传输的方式，
* 循环中消息也是从小序号传递到大序号，
* 从而使得每个状态中，状态的号码等于最大的分块所在的进程号 */