并行计算矩阵分块乘法

目录
一、题目及要求 (1)
1、题目 (1)
2、要求 (1)
二、设计算法、算法原理 (1)
三、算法描述、设计流程 (2)
3.1算法描述 (2)
3.2设计流程 (4)
四、源程序代码及运行结果 (6)
1、超立方 (6)
1.1超立方的源程序代码 (6)
1.2运行结果 (11)
2、网孔连接 (11)
2.1源程序代码 (11)
2.2运行结果 (18)
3、在数学软件中的计算结果 (19)
五、算法分析、优缺点 (19)
1、简单的并行分块乘法的过程为 (19)
2、使用Cannon算法时的算法分析 (20)
3、算法的优缺点 (21)
六、总结 (22)
参考文献 (23)
一、题目及要求
1、题目
简单并行分块乘法：（1）情形1: 超立方连接；（2）情形2：二维环绕网孔连接
已知,177511195310135411274329,7563895712314
2120143321
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛----=⎪⎪⎪⎪⎪⎭⎫
⎝
⎛----=B A 求B A C ⨯=。

2、要求
（1）题目分析、查阅与题目相关的资料； （2）设计算法；
（3）算法实现、代码编写； （4）结果分析和性能分析与改进； （5）论文撰写、答辩；
二、设计算法、算法原理
要考虑的计算问题是C=AB,其中A 与B 分别是n n ⨯矩阵。

①A 、B 和C 分成p p p ⨯=的方块阵ij A ,ij B 和ij C ,大小均为
p
n
p n ⨯ ，p 个处理器编号为1
,1, (1)
0,....,0,0---p p p p
p p , ij P 存放ij A ,ij B 和ij C 。

②通讯:每行处理器进行A 矩阵块的多到多播送(得到ik A , k=0~1-p ) 每列处理器进行B 矩阵块的多到多播送(得到kj B , k=0~ 1-p )
③乘-加运算: ij P 做kj p k ik
ij B A
C ∑-==
1
三、算法描述、设计流程
3.1算法描述
超立方情形下矩阵的简单并行分块算法 输入：待选路的信包在源处理器中 输出：将原处理器中的信包送至其目的地 Begin
（1） for i=1 to n do
11--⊗=i i i d s r
endfor
(2) S V i ==,1 (3) while n i ≤do
(3.1)if 1=i r then 从当前节点V 选路到节点为V ⊗1 （3.2）1+=i i endwhile End
二维网孔情形下矩阵的简单并行分块算法 输入：待选路的信包处于源处理器中 输出：将各信包送至各自的目的地中 Begin
（1） 沿x 维将信包向左或向右选路至目的地的处理器所在的列 （2） 沿y 维将信包向上或向下选路至目的地的处理器所在的行 分块乘法算法
//输入: n n A ⨯,n n B ⨯ ; 子快大小均为
p
n p
n ⨯
输出: n n C ⨯n
Begin
(1)for i=0 to 1-p do for all par-do ij p if i>k then ij A ←()m od ,1j i A +
endif
if j>k then
ij B ← B (i+1)mod , j endif endfor endfor
for i=0 to 1-p do for all ij p par-do ij C =ij A +ij B endfor Endfor End
3.2设计流程
以下是二维网孔与超立方连接设计流程。

如图3-1 二维网孔 步骤：
(1)先进行行播送; (2)再同时进行列播送;
图3-1 二维网孔示意图
4
4 3
超立方
步骤：依次从低维到高维播送, d-立方, d=0,1,2,3,4…; 算法流程如图所示：
图3-2 算法流程
四、源程序代码及运行结果
1、超立方
1.1超立方的源程序代码
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#define intsize sizeof(int)
#define floatsize sizeof(float)
#define charsize sizeof(char)
#define A(x,y) A[x*K+y]
#define B(x,y) B[x*N+y]
#define C(x,y) C[x*N+y]
#define a(x,y) a[x*K+y]
#define b(x,y) b[x*n+y]
#define buffer(x,y) buffer[x*n+y]
#define c(l,x,y) c[x*N+y+l*n]
float *a,*b,*c,*buffer;
int s;
float *A,*B,*C;
int M,N,K,P ;
int m,n;
int myid;
int p;
FILE *dataFile;
MPI_Status status;
double time1;
double starttime,endtime;
void readData()
{
int i,j;
starttime = MPI_Wtime();
dataFile=fopen("yin.txt","r");
fscanf(dataFile,"%d%d", &M, &K); A=(float *)malloc(floatsize*M*K); for(i = 0; i < M; i++) {
for(j = 0; j < K; j++)
{
fscanf(dataFile,"%f", A+i*K+j);
}
}
fscanf(dataFile,"%d%d", &P, &N); if (K!=P) {
printf("the input is wrong\n");
exit(1);
}
B=(float *)malloc(floatsize*K*N); for(i = 0; i < K; i++) {
for(j = 0; j < N; j++)
{
fscanf(dataFile,"%f", B+i*N+j);
}
}
fclose(dataFile);
printf("Input of file \"yin.txt\"\n");
printf("%d\t %d\n",M, K); for(i=0;i<M;i++) {
for(j=0;j<K;j++) printf("%f\t",A(i,j));
printf("\n");
}
printf("%d\t %d\n",K, N); for(i=0;i<K;i++) {
for(j=0;j<N;j++) printf("%f\t",B(i,j));
printf("\n");
}
C=(float *)malloc(floatsize*M*N); }
int gcd(int M,int N,int group_size)
{
int i;
for(i=M; i>0; i--)
{
if((M%i==0)&&(N%i==0)&&(i<=group_size))
return i;
}
return 1;
}
void printResult()
{
int i,j;
printf("\nOutput of Matrix C = AB\n");
for(i=0;i<M;i++)
{
for(j=0;j<N;j++) printf("%f\t",C(i,j));
printf("\n");
}
endtime=MPI_Wtime();
printf("\n");
printf("Whole running time = %f seconds\n",endtime-starttime); printf("Distribute data time = %f seconds\n",time1-starttime); printf("Parallel compute time = %f seconds\n",endtime-time1);
}
int main(int argc, char **argv)
{
int i,j,k,l,group_size,mp1,mm1;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
p=group_size;
if(myid==0)
{
readData();
}
if (myid==0)
for(i=1;i<p;i++)
{
MPI_Send(&M,1,MPI_INT,i,i,MPI_COMM_WORLD);
MPI_Send(&K,1,MPI_INT,i,i,MPI_COMM_WORLD);
MPI_Send(&N,1,MPI_INT,i,i,MPI_COMM_WORLD);
}
else
{
MPI_Recv(&M,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
MPI_Recv(&K,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
MPI_Recv(&N,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
}
p=gcd(M,N,group_size);
m=M/p;
n=N/p;
if(myid<p)
{
a=(float *)malloc(floatsize*m*K);
b=(float *)malloc(floatsize*K*n);
c=(float *)malloc(floatsize*m*N);
if (myid%2!=0)
buffer=(float *)malloc(K*n*floatsize);
if (a==NULL||b==NULL||c==NULL)
printf("Allocate space for a,b or c fail!");
if (myid==0)
{
for (i=0;i<m;i++)
for (j=0;j<K;j++)
a(i,j)=A(i,j);
for (i=0;i<K;i++)
for (j=0;j<n;j++)
b(i,j)=B(i,j);
}
if (myid==0)
{
for (i=1;i<p;i++)
{
MPI_Send(&A(m*i,0),K*m,MPI_FLOAT,i,i,MPI_COMM_WORLD); for (j=0;j<K;j++)
MPI_Send(&B(j,n*i),n,MPI_FLOAT,i,i,MPI_COMM_WORLD);
}
free(A);
free(B);
}
else
{
MPI_Recv(a,K*m,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status); for (j=0;j<K;j++)
MPI_Recv(&b(j,0),n,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
}
if (myid==0)
time1=MPI_Wtime();
for (i=0;i<p;i++)
{
l=(i+myid)%p;
for (k=0;k<m;k++)
for (j=0;j<n;j++)
for (c(l,k,j)=0,s=0;s<K;s++)
c(l,k,j)+=a(k,s)*b(s,j);
mm1=(p+myid-1)%p;
mp1=(myid+1)%p;
if (i!=p-1)
{
if(myid%2==0)
{
MPI_Send(b,K*n,MPI_FLOAT,mm1,mm1,MPI_COMM_WORLD);
MPI_Recv(b,K*n,MPI_FLOAT,mp1,myid,MPI_COMM_WORLD,&status);
}
else
{
for(k=0;k<K;k++)
for(j=0;j<n;j++)
buffer(k,j)=b(k,j);
MPI_Recv(b,K*n,MPI_FLOAT,mp1,myid,MPI_COMM_WORLD,&status);
MPI_Send(buffer,K*n,MPI_FLOAT,mm1,mm1,MPI_COMM_WORLD);
}
}
}
if (myid==0)
for(i=0;i<m;i++)
for(j=0;j<N;j++)
C(i,j)=*(c+i*N+j);
if (myid!=0)
MPI_Send(c,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD);
else
{
for(k=1;k<p;k++)
{
MPI_Recv(c,m*N,MPI_FLOAT,k,k,MPI_COMM_WORLD,&status); for(i=0;i<m;i++)
for(j=0;j<N;j++)
C((k*m+i),j)=*(c+i*N+j);
}
}
if(myid==0)
printResult();
}
MPI_Finalize();
if(myid<p)
{
free(a);
free(b);
free(c);
if(myid==0)
free(C);
if(myid%2!=0)
free(buffer);
}
return (0);
}
1.2运行结果
图4.1 4个处理器的运行结果
2、网孔连接
2.1源程序代码
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include <time.h>
#include <stdio.h>
#include <math.h>
/* 全局变量声明 */
float **A, **B, **C; /* 总矩阵,C = A * B */
float *a, *b, *c, *tmp_a, *tmp_b; /* a、b、c表分块，tmp_a、tmp_b表缓冲区 */
int dg, dl, dl2,p, sp; /* dg:总矩阵维数;dl:矩阵块维数;dl2=dl*dl;p:处理器个数;sp＝sqrt(p) */
int my_rank, my_row, my_col; /* my_rank:处理器ID;(my_row,my_col):处理器逻辑阵列坐标 */
MPI_Status status;
/*
*函数名: get_index
*功能：处理器逻辑阵列坐标至rank号的转换
*输入：坐标、逻辑阵列维数
*输出：rank号
*/
int get_index(int row, int col, int sp)
{
return ((row+sp)%sp)*sp + (col+sp)%sp;
}
/*
*函数名：random_A_B
*功能：随机生成矩阵A和B
*/
void random_A_B()
{
int i,j;
float m;
//srand((unsigned int)time(NULL)); /*设随机数种子*
/*随机生成A,B,并初始化C*/
for(i=0; i<dg ; i++)
for(j=0; j<dg ; j++)
{
scanf("%f",&m);
A[i][j] = m;
C[i][j] = 0.0;
m=0;
}
for(i=0; i<dg ; i++)
for(j=0; j<dg ; j++)
{
scanf("%f",&m);
B[i][j] = m;
m=0;
}
}
/* 函数名：scatter_A_B
* 功能：rank为0的处理器向其他处理器发送A、B矩阵的相关块
*/
void scatter_A_B()
{
int i,j,k,l;
int p_imin,p_imax,p_jmin,p_jmax;
for(k=0; k<p; k++)
{
/*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/
p_jmin = (k % sp ) * dl;
p_jmax = (k % sp + 1) * dl-1;
p_imin = (k - (k % sp))/sp * dl;
p_imax = ((k - (k % sp))/sp +1) *dl -1;
l = 0;
/*rank=0的处理器将A,B中的相应块拷至tmp_a,tmp_b，准备向其他处理器发送*/
for(i=p_imin; i<=p_imax; i++)
{
for(j=p_jmin; j<=p_jmax; j++)
{
tmp_a[l] = A[i][j];
tmp_b[l] = B[i][j];
l++;
}
}
/*rank=0的处理器直接将自己对应的矩阵块从tmp_a,tmp_b拷至a,b*/ if(k==0)
{
memcpy(a, tmp_a, dl2 * sizeof(float));
memcpy(b, tmp_b, dl2 * sizeof(float));
} else /*rank=0的处理器向其他处理器发送tmp_a,tmp_b中相关的矩阵块*/
{
MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD);
}
}
}
/*
*函数名:init_alignment
*功能:矩阵A和B初始对准
*/
void init_alignment()
{
MPI_Sendrecv(a, dl2, MPI_FLOAT, get_index(my_row,my_col-my_row,sp), 1,
tmp_a, dl2, MPI_FLOAT, get_index(my_row,my_col+my_row,sp), 1, MPI_COMM_WORLD, &status);
memcpy(a, tmp_a, dl2 * sizeof(float) );
/*将B中坐标为(i,j)的分块B(i,j)向上循环移动j步*/
MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1,
tmp_b, dl2, MPI_FLOAT, get_index(my_row+my_col,my_col,sp), 1, MPI_COMM_WORLD, &status);
memcpy(b, tmp_b, dl2 * sizeof(float) );
}
/*
*函数名：main_shift
*功能：分块矩阵左移和上移，并计算分块c
*/
void main_shift()
{
int i,j,k,l;
for(l=0; l<sp; l++)
{
/*矩阵块相乘，c+=a*b */
for(i=0; i<dl; i++)
for(j=0; j<dl; j++)
for(k=0; k<dl; k++)
c[i*dl+j] += a[i*dl+k]*b[k*dl+j];
/* 将分块a左移1位 */
MPI_Send(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD);
MPI_Recv(a , dl2, MPI_FLOAT, get_index(my_row, my_col+1, sp), 1, MPI_COMM_WORLD, &status);
/* 将分块b上移1位 */
MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD);
MPI_Recv(b , dl2, MPI_FLOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status);
}
}
/*
*函数名：collect_c
*功能：rank为0的处理器从其余处理器收集分块矩阵c
*/
void collect_C()
{
int i,j,i2,j2,k;
int p_imin,p_imax,p_jmin,p_jmax; /* 分块矩阵在总矩阵中顶点边界值 */
/* 将rank为0的处理器中分块矩阵c结果赋给总矩阵C对应位置 */
for (i=0;i<dl;i++)
for(j=0;j<dl;j++)
C[i][j]=c[i*dl+j];
for (k=1;k<p;k++)
{
/*将rank为0的处理器从其他处理器接收相应的分块c*/
MPI_Recv(c, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status);
p_jmin = (k % sp ) *dl;
p_jmax = (k % sp + 1) *dl-1;
p_imin = (k - (k % sp))/sp *dl;
p_imax = ((k - (k % sp))/sp +1) *dl -1;
i2=0;
/*将接收到的c拷至C中的相应位置,从而构造出C*/
for(i=p_imin; i<=p_imax; i++)
{
j2=0;
for(j=p_jmin; j<=p_jmax; j++)
{
C[i][j]=c[i2*dl+j2];
j2++;
}
i2++;
}
}
}
/*函数名：print
*功能：打印矩阵
*输入：指向矩阵指针的指针，字符串
*/
void print(float **m,char *str)
{
int i,j;
printf("%s",str);
/*打印矩阵m*/
for(i=0;i<dg;i++)
{
for(j=0;j<dg;j++)
printf("%15.0f ",m[i][j]);
printf("\n");
}
printf("\n");
}
/*
*函数名：main
*功能：主过程，Cannon算法，矩阵相乘
*输入：argc为命令行参数个数，argv为每个命令行参数组成的字符串数组 */
int main(int argc, char *argv[])
{
int i;
MPI_Init(&argc, &argv); /* 启动MPI计算 */
MPI_Comm_size(MPI_COMM_WORLD, &p); /* 确定处理器个数 */
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); /* 确定各自的处理器标识符*/
sp = sqrt(p);
/* 确保处理器个数是完全平方数，否则打印错误信息，程序退出 */
if (sp*sp != p)
{
if (my_rank == 0)
printf("Number of processors is not a quadratic number!\n");
MPI_Finalize();
exit(1);
}
if (argc != 2)
{
if (my_rank == 0)
printf("usage: mpirun -np ProcNum cannon MatrixDimension\n"); MPI_Finalize();
exit(1);
}
dg = atoi(argv[1]); /* 总矩阵维数 */
dl = dg / sp; /* 计算分块矩阵维数 */
dl2 = dl * dl;
/* 计算处理器在逻辑阵列中的坐标 */
my_col = my_rank % sp ;
my_row = (my_rank-my_col) / sp ;
/* 为a、b、c分配空间 */
a = (float *)malloc( dl2 * sizeof(float) );
b = (float *)malloc( dl2 * sizeof(float) );
c = (float *)malloc( dl2 * sizeof(float) );
/* 初始化c */
for(i=0; i<dl2 ; i++)
c[i] = 0.0;
/* 为tmp_a、tmp_b分配空间 */
tmp_a = (float *)malloc( dl2 * sizeof(float) );
tmp_b = (float *)malloc( dl2 * sizeof(float) );
if (my_rank == 0)
{
/* rank为0的处理器为A、B、C分配空间 */
A = (float **)malloc( dg * sizeof(float*) );
B = (float **)malloc( dg * sizeof(float*) );
C = (float **)malloc( dg * sizeof(float*) );
for(i=0; i<dg; i++)
{
A[i] = (float *)malloc( dg * sizeof(float) );
B[i] = (float *)malloc( dg * sizeof(float) );
C[i] = (float *)malloc( dg * sizeof(float) );
}
random_A_B(); /* rank为0的处理器随机化生成A、B矩阵 */
scatter_A_B(); /* rank为0的处理器向其他处理器发送A、B矩阵的相关块 */
}
else /* rank不为0的处理器接收来自rank为0的处理器的相应矩阵分块 */
{
MPI_Recv(a, dl2, MPI_FLOAT, 0 , 1, MPI_COMM_WORLD, &status); MPI_Recv(b, dl2, MPI_FLOAT, 0 , 2, MPI_COMM_WORLD, &status);
}
init_alignment(); /* A、B矩阵的初始对准 */
main_shift(); /* 分块矩阵左移、上移, cannon算法的主过程 */
if(my_rank == 0)
{
collect_C(); /* rank为0的处理器从其余处理器收集分块矩阵c */
print(A,"random matrix A : \n"); /* 打印矩阵A */
print(B,"random matrix B : \n"); /* 打印矩阵B */
print(C,"Matrix C = A * B : \n"); /* 打印矩阵C */
} else
{
MPI_Send(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); }
MPI_Barrier(MPI_COMM_WORLD); /* 同步所有处理器 */
MPI_Finalize(); /* 结束MPI计算 */
return 0;
}
2.2运行结果
图4.2 4个处理器的运行结果
3、在数学软件中的计算结果
图4.3 在MATLAB 中的运行结果
五、算法分析、优缺点
1、简单的并行分块乘法的过程为
（1）分块：将: A n ×n 与 B n ×n 分成p 块A i,j 和B i,j (0≤i,j ≤1-p ),每块大小为)/()/(p n p n ⨯，并将它们分配给p p ⨯个处理器 （1,11,00,0,...,,...,---p p p P P P ）。

开始时P i,j 存放在块A i,j 与块B i,j ，然后计算块C i,j 。

（2）通信：为了计算块C i,j ，总结需要所有块A i,j 与块B i,j (0≤k ≤1-p ),为
此在每行处理器之间要进行A 矩阵块的多到多播送，已得到A i,k ；而在每列的处
理器之间要进行B 矩阵块的多到多播送，已得到B k,j 。

（3）乘-加计算:一旦P i,j 收到A i,0, A i,1 (i)
p -1与B 0,j B 0,j ......B p -1,j,就
同时进行乘-加运算。

当并行分块乘法执行在p 个处理器的超立方上时，通信时间为p n t p t p t p n t w s w s 2
222)1((2+=-+），计算时间为p n p
n p q n q /)()(333=⨯=⨯。

所以算法的并行时间为p
n t p t p n T w s p 2
32log ++= 当并行分块乘法执行在⨯p p 的二维网孔上时，其通信时间为
p T =p
n t p t p n w s 2
322++ 2、使用Cannon 算法时的算法分析
当算法执行在p 个处理器的超立方上时，若使用CT 选路法，算法第一步的循环移位时间为)log (22p t p n t t h w s ++，算法第三步的单步移位时间为p p
n t t w s )(22
+，当p 充分大时后者为主项。

算法第三步共执行p 次⨯)/(p n )/(p n 子块乘法，其时间为
p n p n p /)/(33=⨯。

所以在超立方上Cannon 乘法的并行运行时间为
p
n t t p p n T w s p 2
322++= 当算法执行在⨯p p 的二维网孔上时，算法在第一步的循环移位时间为
p p n t t w s )(22+，算法在第三步的单步移位时间为p p
n t t w s )(22
+，算法在第三步的运算时间为p n /3。

所以在二维网孔上Cannon 乘法的并行运行时间为 p
n t p t p n T w s p 2
344++= 使用Cannon 算法的时间分析为
)/())/(()1()()(33321p n O p n p O O p O T T T n T p =⋅++=++=
3、算法的优缺点 简单并行分块乘法的主要缺点是存储要求过大：通信结束时每个处理器拥有p 2块，
每块大小为p n /2，所以每个处理器要求）（p n /O 2存储器，p 个处理器共要求）（p n /O 2，它是串行算法存储要求的
p 倍。

当p=n 2时，t(n)=O(n), c(n)=O(n 3)
通过此次课程设计，使我更加懂得运用所学的知识，在设计过程中虽然遇到了一些问题，但经过一次又一次的思考，一遍又一遍的检查终于找出了原因所在，也暴露了前期我在这方面的知识欠缺和经验不足。

实践出真知，通过亲手动手制作，使我们掌握的知识不再是纸上谈兵。

我们本次做的题目是矩阵的简单并行分块乘法，使用超立方连接和二维网孔连接，刚拿到题目感觉很难，然后通过查找资料，询问同学及老师终于有了头绪。

然后在课设过程中按部就班进行算法设计，终于实现了。

在课程设计过程中，我不断发现错误，不断改正，不断领悟，不断获取。

最终的检测调试环节，本身就是在践行“过而能改，善莫大焉”的行知观。

这次课程设计终于顺利完成了，在设计过程中遇到了很多问题，最后在老师的指导下，终于游逆而解。

在今后社会的发展和学习实践过程中，一定要不懈努力，不能遇到问题就想要退缩。

当然在这个过程中我求教了很多人，老师和同学的帮助才使本次课设顺利完成。

[1] 陈国良, 吴俊敏, 章锋等, 并行计算机体系结构[M], 北京: 高等教育出版
社, 2002.
[2] 陈国良, 并行算法的设计与分析[M]第三版, 北京: 高等教育出版社,
2009.
[3] 陈国良等, 并行算法实践[M], 北京: 高等教育出版社, 2003.
[4] 陈国良, 并行计算----结构·算法·编程[M], 第三版北京：高等教育出版
社, 2011.。