矩阵QR分解并行实现Word格式.docx
《矩阵QR分解并行实现Word格式.docx》由会员分享，可在线阅读，更多相关《矩阵QR分解并行实现Word格式.docx（16页珍藏版）》请在冰点文库上搜索。
doubletime1;
doubletime2;
intp;
MPI_Statusstatus;
voidEnvironment_Finalize（float*a,float*q,float*v,float*f,float*R,
float*Q,float*ai,float*aj,float*qi,float*qj）
{
free（a）;
free（q）;
free（v）;
free（f）;
free（R）;
free（Q）;
free（ai）;
free（aj）;
free（qi）;
free（qj）;
}
intmain（intargc,char**argv）
intM,N,m;
intz;
inti,j,k,my_rank,group_size;
float*ai,*qi,*aj,*qj;
floatc,s,sp;
float*f,*v;
float*a,*q;
FILE*fdA;
MPI_Init（&
argc,&
argv）;
MPI_Comm_rank（MPI_COMM_WORLD,&
my_rank）;
MPI_Comm_size（MPI_COMM_WORLD,&
group_size）;
p=group_size;
starttime=MPI_Wtime（）;
if（my_rank==p-1）
{
fdA=fopen（"
dataIn.txt"
"
r"
）;
fscanf（fdA,"
%d%d"
&
M,&
N）;
if（M!
=N）
puts（"
Theinputiserror!
"
exit（0）;
}
A=（float*）malloc（sizeof（float）*M*M）;
Q=（float*）malloc（sizeof（float）*M*M）;
R=（float*）malloc（sizeof（float）*M*M）;
for（i=0;
i<
M;
i++）
for（j=0;
j<
j++）fscanf（fdA,"
%f"
A+i*M+j）;
fclose（fdA）;
for（i=0;
i<
M;
i++）
for（j=0;
j<
j++）if（i==j）Q（i,j）=1.0;
elseQ（i,j）=0.0;
MPI_Bcast（&
M,1,MPI_INT,p-1,MPI_COMM_WORLD）;
m=M/p;
if（M%p!
=0）m++;
qi=（float*）malloc（sizeof（float）*M）;
qj=（float*）malloc（sizeof（float）*M）;
aj=（float*）malloc（sizeof（float）*M）;
ai=（float*）malloc（sizeof（float）*M）;
v=（float*）malloc（sizeof（float）*M）;
f=（float*）malloc（sizeof（float）*M）;
a=（float*）malloc（sizeof（float）*m*M）;
q=（float*）malloc（sizeof（float）*m*M）;
if（a==NULL||q==NULL||f==NULL||v==NULL||qi==NULL||qj==NULL||
ai==NULL||aj==NULL）
printf（"
memoryallocationiswrong\n"
if（my_rank==p-1）
m;
j++）
a（i,j）=A（（my_rank*m+i）,j）;
q（i,j）=Q（（my_rank*m+i）,j）;
i<
p-1;
i++）
MPI_Send（&
A（m*i,0）,m*M,MPI_FLOAT,i,i,MPI_COMM_WORLD）;
Q（m*i,0）,m*M,MPI_FLOAT,i,i,MPI_COMM_WORLD）;
free（A）;
else{
MPI_Recv（a,m*M,MPI_FLOAT,p-1,my_rank,MPI_COMM_WORLD,&
status）;
MPI_Recv（q,m*M,MPI_FLOAT,p-1,my_rank,MPI_COMM_WORLD,&
time1=MPI_Wtime（）;
if（p>
1）
if（my_rank==0）
j<
m-1;
j++）
for（i=j+1;
sp=sqrt（a（j,j）*a（j,j）+a（i,j）*a（i,j））;
c=a（j,j）/sp;
s=a（i,j）/sp;
for（k=0;
k<
k++）
aj[k]=c*a（j,k）+s*a（i,k）;
qj[k]=c*q（j,k）+s*q（i,k）;
ai[k]=-s*a（j,k）+c*a（i,k）;
qi[k]=-s*q（j,k）+c*q（i,k）;
a（j,k）=aj[k];
q（j,k）=qj[k];
a（i,k）=ai[k];
q（i,k）=qi[k];
}/*i*/
f[k]=a（j,k）;
v[k]=q（j,k）;
f[0],M,MPI_FLOAT,1,j,MPI_COMM_WORLD）;
v[0],M,MPI_FLOAT,1,j,MPI_COMM_WORLD）;
}/*forj*/
f[k]=a（（m-1）,k）;
v[k]=q（（m-1）,k）;
f[0],M,MPI_FLOAT,1,m-1,MPI_COMM_WORLD）;
v[0],M,MPI_FLOAT,1,m-1,MPI_COMM_WORLD）;
}/*my_rank==0*/
else/*my_rank!
=0*/
if（my_rank!
=（group_size-1））
my_rank*m;
MPI_Recv（&
f[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&
v[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&
sp=sqrt（f[j]*f[j]+a（i,j）*a（i,j））;
c=f[j]/sp;
aj[k]=c*f[k]+s*a（i,k）;
qj[k]=c*v[k]+s*q（i,k）;
ai[k]=-s*f[k]+c*a（i,k）;
qi[k]=-s*v[k]+c*q（i,k）;
f[k]=aj[k];
v[k]=qj[k];
f[0],M,MPI_FLOAT,（my_rank+1）,j,MPI_COMM_WORLD）;
v[0],M,MPI_FLOAT,（my_rank+1）,j,MPI_COMM_WORLD）;
sp=sqrt（a（j,（my_rank*m+j））*a（j,（my_rank*m+j））+a（i,（my_rank*m+j））*a（i,（my_rank*m+j）））;
c=a（j,（my_rank*m+j））/sp;
s=a（i,（my_rank*m+j））/sp;
f[0],M,MPI_FLOAT,my_rank+1,my_rank*m+j,MPI_COMM_WORLD）;
v[0],M,MPI_FLOAT,my_rank+1,my_rank*m+j,MPI_COMM_WORLD）;
f[0],M,MPI_FLOAT,my_rank+1,my_rank*m+m-1,MPI_COMM_WORLD）;
MPI_Send（&
v[0],M,MPI_FLOAT,my_rank+1,my_rank*m+m-1,MPI_COMM_WORLD）;
}/*my_rank!
=groupsize-1*/
if（my_rank==（group_size-1））
MPI_Recv（&
}/*fori*/
Q（j,k）=v[k];
R（j,k）=f[k];
Q（（my_rank*m+j）,k）=q（j,k）;
R（（my_rank*m+j）,k）=a（j,k）;
Q（（my_rank*m+m-1）,k）=q（（m-1）,k）;
R（（my_rank*m+m-1）,k）=a（（m-1）,k）;
}/*formy_rank==groupsize-1*/
}/*elsemy_rank!
}/*ifp>
1*/
if（p==1）
for（j=0;
for（i=j+1;
sp=sqrt（a（j,j）*a（j,j）+a（i,j）*a（i,j））;
for（k=0;
aj[k]=c*a（j,k）+s*a（i,k）;
qj[k]=c*q（j,k）+s*q（i,k）;
ai[k]=（-s）*a（j,k）+c*a（i,k）;
qi[k]=（-s）*q（j,k）+c*q（i,k）;
}/*for*/
R（i,j）=a（i,j）;
Q（i,j）=q（i,j）;
}/*ifp==1*/
printf（"
Inputoffile\"
dataIn.txt\"
\n"
%d\t%d\n"
M,N）;
N;
j++）printf（"
%f\t"
A（i,j））;
\nOutputofQRoperation\n"
MatrixR:
R（i,j））;
for（j=i+1;
temp=Q（i,j）;
Q（i,j）=Q（j,i）;
Q（j,i）=temp;
MatrixQ:
Q（i,j））;
time2=MPI_Wtime（）;
Wholerunningtime=%fseconds\n"
time2-starttime）;
Distributedatatime=%fseconds\n"
time1-starttime）;
Parallelcomputetime=%fseconds\n"
time2-time1）;
MPI_Barrier（MPI_COMM_WORLD）;
MPI_Finalize（）;
Environment_Finalize（a,q,v,f,R,Q,ai,aj,qi,qj）;
return（0）;
编译：
mpiccch_qr.c-oqr-lm
运行：
mpirun-np4qr
dataIn.txt里数据：
1.0000002.0000001.000000
2.0000001.0000001.000000
实验结果：
1、编译成功：
2、运行成果：
个人总结：
矩阵QR分解是研究生“随机过程”中较为重要的一个内容，试验中将涉及到的学生知识转化成程序有一定的难度。
试验中遇到了许多问题，由于这方面的经验比较缺少，特别是并行实验平台上进行编程，更是一大挑战。
所以摸索的过程有点艰苦。
通过此次实验，对并行编程方法有了一定的了解。
希望能够通过今后的学习来取得大的进步，对KD60实验平台有更深入的了解。
升级会员 