北航并行计算矩阵相乘作业

矩阵相乘并行实现

1、算法描述：

设有如下矩阵相乘：

C=A×B

其中A，B分别是m×k和k×n矩阵，C是m×n矩阵。若处理器个数为p，且它们的编号依次是0,1，…,p-1,则设可将矩阵A、B、C分成p个大小为mxm的子

块，其中A=(Aij)m

×m ，B=(Bij)m

×m

，和C=(Cij)m

×m

，其中A¬ij，Bij和Cij是n×n矩

阵。同时假设和。

定义对角块矩阵，则

其中，。利用此关系式，将节点编号从一维映射到二维，数据，，存放在中，可得到下面的在处理机

结点上的算法。该算法数据交量

算法流程如下：

流程图如下所示：

2、程序代码：

#include

#include

#include

#include

#include

#include

float **A, **B, **C;

float *a, *b, *c, *tmp_a, *tmp_b;

int dg=1000, dl, dl2,p, sp;

int my_rank, my_row, my_col;

MPI_Status status;

int get_index(int row, int col, int sp)

{

return ((row+sp)%sp)*sp + (col+sp)%sp; }

void random_A_B()

{

int i,j;

srand((unsigned int)time(NULL));

for(i=0; i

for(j=0; j

{

A[i][j] = rand();

B[i][j] = rand();

C[i][j] = 0.0;

}

}

void scatter_A_B()

{

int i,j,k,l;

int p_imin,p_imax,p_jmin,p_jmax;

for(k=0; 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;

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++;

}

}

if(k==0)

{

memcpy(a, tmp_a, dl2 * sizeof(float));

memcpy(b, tmp_b, dl2 * sizeof(float));

} else

{

MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD);

MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD);

}

}

}

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) );

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) );

}

void main_shift()

{

int i,j,k,l;

for(l=0; l

{

for(i=0; i

for(j=0; j

for(k=0; k

c[i*dl+j] += a[i*dl+k]*b[k*dl+j];

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);

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);

}

}

void collect_C()

{

int i,j,i2,j2,k;

int p_imin,p_imax,p_jmin,p_jmax;

for (i=0;i

for(j=0;j

C[i][j]=c[i*dl+j];