并行计算矩阵分块乘法.docx
《并行计算矩阵分块乘法.docx》由会员分享,可在线阅读,更多相关《并行计算矩阵分块乘法.docx(32页珍藏版)》请在冰豆网上搜索。
并行计算矩阵分块乘法
目录
一、题目及要求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:
二维环绕网孔连接
已知
求
。
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-1二维网孔示意图
超立方
步骤:
依次从低维到高维播送,d-立方,d=0,1,2,3,4…;
算法流程如图所示:
图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;jprintf("\n");
}
printf("%d\t%d\n",K,N);
for(i=0;i{
for(j=0;jprintf("\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;jprintf("\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;ifor(j=0;ja(i,j)=A(i,j);
for(i=0;ifor(j=0;jb(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;jMPI_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;jMPI_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;kfor(j=0;jfor(c(l,k,j)=0,s=0;sc(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;kfor(j=0;jbuffer(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;ifor(j=0;jC(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;ifor(j=0;jC((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;ifor(j=0;j{
scanf("%f",&m);
A[i][j]=m;
C[i][j]=0.0;
m=0;
}
for(i=0;ifor(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;jprintf("%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);
/