并行计算矩阵分块乘法.docx

并行计算矩阵分块乘法.docx

《并行计算矩阵分块乘法.docx》由会员分享，可在线阅读，更多相关《并行计算矩阵分块乘法.docx（32页珍藏版）》请在冰点文库上搜索。

并行计算矩阵分块乘法

一、题目及要求

1、题目

简单并行分块乘法：

（1）情形1:

超立方连接；

（2）情形2：

二维环绕网孔连接

已知

求

。

2、要求

（1）题目分析、查阅与题目相关的资料；

（2）设计算法；

（3）算法实现、代码编写；

（4）结果分析和性能分析与改进；

（5）论文撰写、答辩；

二、设计算法、算法原理

要考虑的计算问题是C=AB,其中A与B分别是

矩阵。

A、B和C分成

的方块阵

和

大小均为

，p个处理器编号为

存放

和

。

通讯:

每行处理器进行A矩阵块的多到多播送（得到

k=0~

）

每列处理器进行B矩阵块的多到多播送（得到

k=0~

）

乘-加运算:

做

三、算法描述、设计流程

3.1算法描述

超立方情形下矩阵的简单并行分块算法

输入：

待选路的信包在源处理器中

输出：

将原处理器中的信包送至其目的地

Begin

（1）fori=1tondo

endfor

（2）

（3）while

do

（3.1）if

then从当前节点V选路到节点为V

1

（3.2）

endwhile

End

二维网孔情形下矩阵的简单并行分块算法

输入：

待选路的信包处于源处理器中

输出：

将各信包送至各自的目的地中

Begin

（1）沿x维将信包向左或向右选路至目的地的处理器所在的列

（2）沿y维将信包向上或向下选路至目的地的处理器所在的行

分块乘法算法

//输入:

;子快大小均为

输出:

n

Begin

（1）fori=0to

do

forallpar-do

ifi>kthen

endif

ifj>kthen

B（i+1）mod,j

endif

endfor

endfor

fori=0to

do

forall

par-do

=

+

endfor

Endfor

End

3.2设计流程

以下是二维网孔与超立方连接设计流程。

如图3-1

二维网孔

步骤：

（1）先进行行播送;

（2）再同时进行列播送;

3

7

6

2

5

1

4

9

8

0

4

4

4

4

4

4

4

4

4

2

2

3

3

3

3

1

11

12

13

14

15

10

图3-1二维网孔示意图

超立方

步骤：

依次从低维到高维播送,d-立方,d=0,1,2,3,4…;

算法流程如图所示：

包括mpi的头文件

相关变量声明

MPI_INIT（）

MPI_COMM_RANK（）

MPI_COMM_SIZE（）

进入MPI系统

矩阵内部的通信

应用控制实体：

矩阵内部的计算程序

MPI_FINALIZE（）

退出MPI系统

结束

开始

循环直至结束

图3-2算法流程

四、源程序代码及运行结果

1、超立方

1.1超立方的源程序代码

#include"stdio.h"

#include"stdlib.h"

#include"mpi.h"

#defineintsizesizeof（int）

#definefloatsizesizeof（float）

#definecharsizesizeof（char）

#defineA（x,y）A[x*K+y]

#defineB（x,y）B[x*N+y]

#defineC（x,y）C[x*N+y]

#definea（x,y）a[x*K+y]

#defineb（x,y）b[x*n+y]

#definebuffer（x,y）buffer[x*n+y]

#definec（l,x,y）c[x*N+y+l*n]

float*a,*b,*c,*buffer;

ints;

float*A,*B,*C;

intM,N,K,P;

intm,n;

intmyid;

intp;

FILE*dataFile;

MPI_Statusstatus;

doubletime1;

doublestarttime,endtime;

voidreadData（）

{

inti,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

{

for（j=0;j

{

fscanf（dataFile,"%f",A+i*K+j）;

}

}

fscanf（dataFile,"%d%d",&P,&N）;

if（K!

=P）

{

printf（"theinputiswrong\n"）;

exit

（1）;

}

B=（float*）malloc（floatsize*K*N）;

for（i=0;i

{

for（j=0;j

{

fscanf（dataFile,"%f",B+i*N+j）;

}

}

fclose（dataFile）;

printf（"Inputoffile\"yin.txt\"\n"）;

printf（"%d\t%d\n",M,K）;

for（i=0;i

{

for（j=0;j

printf（"\n"）;

}

printf（"%d\t%d\n",K,N）;

for（i=0;i

{

for（j=0;j

printf（"\n"）;

}

C=（float*）malloc（floatsize*M*N）;

}

intgcd（intM,intN,intgroup_size）

{

inti;

for（i=M;i>0;i--）

{

if（（M%i==0）&&（N%i==0）&&（i<=group_size））

returni;

}

return1;

}

voidprintResult（）

{

inti,j;

printf（"\nOutputofMatrixC=AB\n"）;

for（i=0;i

{

for（j=0;j

printf（"\n"）;

}

endtime=MPI_Wtime（）;

printf（"\n"）;

printf（"Wholerunningtime=%fseconds\n",endtime-starttime）;

printf（"Distributedatatime=%fseconds\n",time1-starttime）;

printf（"Parallelcomputetime=%fseconds\n",endtime-time1）;

}

intmain（intargc,char**argv）

{

inti,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

{

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

{

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（"Allocatespacefora,borcfail!

"）;

if（myid==0）

{

for（i=0;i

for（j=0;j

a（i,j）=A（i,j）;

for（i=0;i

for（j=0;j

b（i,j）=B（i,j）;

}

if（myid==0）

{

for（i=1;i

{

MPI_Send（&A（m*i,0）,K*m,MPI_FLOAT,i,i,MPI_COMM_WORLD）;

for（j=0;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

MPI_Recv（&b（j,0）,n,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status）;

}

if（myid==0）

time1=MPI_Wtime（）;

for（i=0;i

{

l=（i+myid）%p;

for（k=0;k

for（j=0;j

for（c（l,k,j）=0,s=0;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

for（j=0;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

for（j=0;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

{

MPI_Recv（c,m*N,MPI_FLOAT,k,k,MPI_COMM_WORLD,&status）;

for（i=0;i

for（j=0;j

C（（k*m+i）,j）=*（c+i*N+j）;

}

}

if（myid==0）

printResult（）;

}

MPI_Finalize（）;

if（myid

{

free（a）;

free（b）;

free（c）;

if（myid==0）

free（C）;

if（myid%2!

=0）

free（buffer）;

}

return（0）;

}

1.2运行结果

图4.14个处理器的运行结果

2、网孔连接

2.1源程序代码

#include

#include

#include

#include

#include

#include

/*全局变量声明*/

float**A,**B,**C;/*总矩阵,C=A*B*/

float*a,*b,*c,*tmp_a,*tmp_b;/*a、b、c表分块，tmp_a、tmp_b表缓冲区*/

intdg,dl,dl2,p,sp;/*dg:

总矩阵维数;dl:

矩阵块维数;dl2=dl*dl;p:

处理器个数;sp＝sqrt（p）*/

intmy_rank,my_row,my_col;/*my_rank:

处理器ID;（my_row,my_col）:

处理器逻辑阵列坐标*/

MPI_Statusstatus;

/*

*函数名:

get_index

*功能：

处理器逻辑阵列坐标至rank号的转换

*输入：

坐标、逻辑阵列维数

*输出：

rank号

*/

intget_index（introw,intcol,intsp）

{

return（（row+sp）%sp）*sp+（col+sp）%sp;

}

/*

*函数名：

random_A_B

*功能：

随机生成矩阵A和B

*/

voidrandom_A_B（）

{

inti,j;

floatm;

//srand（（unsignedint）time（NULL））;/*设随机数种子*

/*随机生成A,B,并初始化C*/

for（i=0;i

for（j=0;j

{

scanf（"%f",&m）;

A[i][j]=m;

C[i][j]=0.0;

m=0;

}

for（i=0;i

for（j=0;j

{

scanf（"%f",&m）;

B[i][j]=m;

m=0;

}

}

/*函数名：

scatter_A_B

*功能：

rank为0的处理器向其他处理器发送A、B矩阵的相关块

*/

voidscatter_A_B（）

{

inti,j,k,l;

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

/*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初始对准

*/

voidinit_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

*/

voidmain_shift（）

{

inti,j,k,l;

for（l=0;l

{

/*矩阵块相乘，c+=a*b*/

for（i=0;i

for（j=0;j

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

*/

voidcollect_C（）

{

inti,j,i2,j2,k;

intp_imin,p_imax,p_jmin,p_jmax;/*分块矩阵在总矩阵中顶点边界值*/

/*将rank为0的处理器中分块矩阵c结果赋给总矩阵C对应位置*/

for（i=0;i

for（j=0;j

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

for（k=1;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

*功能：

打印矩阵

*输入：

指向矩阵指针的指针，字符串

*/

voidprint（float**m,char*str）

{

inti,j;

printf（"%s",str）;

/*打印矩阵m*/

for（i=0;i

{

for（j=0;j

printf（"%15.0f",m[i][j]）;

printf（"\n"）;

}

printf（"\n"）;

}

/*

*函数名：

main

*功能：

主过程，Cannon算法，矩阵相乘

*输入：

argc为命令行参数个数，argv为每个命令行参数组成的字符串数组

*/

intmain（intargc,char*argv[]）

{

inti;

MPI_Init（&argc,&argv）;/*启动MPI计算*/

MPI_Comm_size（MPI_COMM_WORLD,&p）;/*确定处理器个数*/

MPI_Comm_rank（MPI_COMM_WORLD,&my_rank）;/*确定各自的处理器标识符*/

sp=sqrt（p）;

/*确保处理器个数是完全平方数，否则打印错误信息，程序退出*/

if（