收藏
下载资源下载资源 加入VIP,免费下载

矩阵QR分解并行实现.docx

上传人：b****2 文档编号：1662295 上传时间：2023-05-01 格式：DOCX 页数：16 大小：67KB
下载 相关 举报
矩阵QR分解并行实现.docx_第1页
第1页 / 共16页
矩阵QR分解并行实现.docx_第2页
第2页 / 共16页
矩阵QR分解并行实现.docx_第3页
第3页 / 共16页
矩阵QR分解并行实现.docx_第4页
第4页 / 共16页
矩阵QR分解并行实现.docx_第5页
第5页 / 共16页
矩阵QR分解并行实现.docx_第6页
第6页 / 共16页
矩阵QR分解并行实现.docx_第7页
第7页 / 共16页
矩阵QR分解并行实现.docx_第8页
第8页 / 共16页
矩阵QR分解并行实现.docx_第9页
第9页 / 共16页
矩阵QR分解并行实现.docx_第10页
第10页 / 共16页
矩阵QR分解并行实现.docx_第11页
第11页 / 共16页
矩阵QR分解并行实现.docx_第12页
第12页 / 共16页
矩阵QR分解并行实现.docx_第13页
第13页 / 共16页
矩阵QR分解并行实现.docx_第14页
第14页 / 共16页
矩阵QR分解并行实现.docx_第15页
第15页 / 共16页
矩阵QR分解并行实现.docx_第16页
第16页 / 共16页
亲，该文档总共16页，全部预览完了，如果喜欢就下载吧！
下载资源
资源描述

矩阵QR分解并行实现.docx

《矩阵QR分解并行实现.docx》由会员分享，可在线阅读，更多相关《矩阵QR分解并行实现.docx（16页珍藏版）》请在冰点文库上搜索。

矩阵QR分解并行实现.docx

矩阵QR分解并行实现

实验18.1矩阵QR分解并行实现

实验要求：

对于一个给定的N*N矩阵，运用Gram-schmidt正交化方法将其进行QR分解。

实验原理：

设A是n阶非奇异矩阵，则存在正交（酉）矩阵Q与实（复）非奇异

上三角矩阵R使得A=Q*R，且除去相差一个对角元素的绝对值（模）全为1的对角因子外，上述分解唯一。

实验目的：

在理解基于消息传递模型的并行程序设计思想，把纯数学的理论—矩阵的QR分解，通过并行编程在计算机上实现，更深层次的了解并行编程方法，并能熟练使用mpi编写并行程序。

源代码：

ch_qr.c

#include"stdio.h"

#include"stdlib.h"

#include"math.h"

#include"mpi.h"

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

#defineq（x,y）q[x*M+y]

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

#defineQ（x,y）Q[x*M+y]

#defineR（x,y）R[x*M+y]

floattemp;

float*A;

float*R;

float*Q;

doublestarttime;

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

{

for（j=0;j

}

fclose（fdA）;

for（i=0;i

for（j=0;j

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）

{

for（i=0;i

for（j=0;j

{

a（i,j）=A（（my_rank*m+i）,j）;

q（i,j）=Q（（my_rank*m+i）,j）;

}

}

if（my_rank==p-1）

{

for（i=0;i

{

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

MPI_Send（&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,&status）;

}

time1=MPI_Wtime（）;

if（p>1）

{

if（my_rank==0）

{

for（j=0;j

{

for（i=j+1;i

{

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

{

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

{

a（j,k）=aj[k];

q（j,k）=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}/*i*/

for（k=0;k

{

f[k]=a（j,k）;

v[k]=q（j,k）;

}

MPI_Send（&f[0],M,MPI_FLOAT,1,j,MPI_COMM_WORLD）;

MPI_Send（&v[0],M,MPI_FLOAT,1,j,MPI_COMM_WORLD）;

}/*forj*/

for（k=0;k

{

f[k]=a（（m-1）,k）;

v[k]=q（（m-1）,k）;

}

MPI_Send（&f[0],M,MPI_FLOAT,1,m-1,MPI_COMM_WORLD）;

MPI_Send（&v[0],M,MPI_FLOAT,1,m-1,MPI_COMM_WORLD）;

}/*my_rank==0*/

else/*my_rank!

=0*/

{

if（my_rank!

=（group_size-1））

{

for（j=0;j

{

MPI_Recv（&f[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&status）;

MPI_Recv（&v[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&status）;

for（i=0;i

{

sp=sqrt（f[j]*f[j]+a（i,j）*a（i,j））;

c=f[j]/sp;s=a（i,j）/sp;

for（k=0;k

{

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

}

for（k=0;k

{

f[k]=aj[k];

v[k]=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}

MPI_Send（&f[0],M,MPI_FLOAT,（my_rank+1）,j,MPI_COMM_WORLD）;

MPI_Send（&v[0],M,MPI_FLOAT,（my_rank+1）,j,MPI_COMM_WORLD）;

}/*forj*/

for（j=0;j

{

for（i=j+1;i

{

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;

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

}

for（k=0;k

{

a（j,k）=aj[k];

q（j,k）=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}

for（k=0;k

{

f[k]=a（j,k）;

v[k]=q（j,k）;

}

MPI_Send（&f[0],M,MPI_FLOAT,my_rank+1,my_rank*m+j,MPI_COMM_WORLD）;

MPI_Send（&v[0],M,MPI_FLOAT,my_rank+1,my_rank*m+j,MPI_COMM_WORLD）;

}/*forj*/

for（k=0;k

{

f[k]=a（（m-1）,k）;

v[k]=q（（m-1）,k）;

}

MPI_Send（&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））

{

for（j=0;j

{

MPI_Recv（&f[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&status）;MPI_Recv（&v[0],M,MPI_FLOAT,（my_rank-1）,j,MPI_COMM_WORLD,&status）;

for（i=0;i

{

sp=sqrt（f[j]*f[j]+a（i,j）*a（i,j））;

c=f[j]/sp;s=a（i,j）/sp;

for（k=0;k

{

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

}

for（k=0;k

{

f[k]=aj[k];

v[k]=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}/*fori*/

for（k=0;k

{

Q（j,k）=v[k];

R（j,k）=f[k];

}

}/*forj*/

for（j=0;j

{

for（i=j+1;i

{

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;

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

}

for（k=0;k

{

a（j,k）=aj[k];

q（j,k）=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}/*fori*/

for（k=0;k

{

Q（（my_rank*m+j）,k）=q（j,k）;

R（（my_rank*m+j）,k）=a（j,k）;

}

}/*forj*/

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

=0*/

}/*ifp>1*/

if（p==1）

{

for（j=0;j

for（i=j+1;i

{

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

{

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

{

a（j,k）=aj[k];

q（j,k）=qj[k];

a（i,k）=ai[k];

q（i,k）=qi[k];

}

}/*for*/

for（i=0;i

for（j=0;j

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

for（i=0;i

for（j=0;j

Q（i,j）=q（i,j）;

}/*ifp==1*/

if（my_rank==p-1）

{

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

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

for（i=0;i

{

for（j=0;j

printf（"\n"）;

}

printf（"\nOutputofQRoperation\n"）;

printf（"MatrixR:

\n"）;

for（i=0;i

{

for（j=0;j

printf（"%f\t",R（i,j））;

printf（"\n"）;

}

for（i=0;i

for（j=i+1;j

{

temp=Q（i,j）;

Q（i,j）=Q（j,i）;

Q（j,i）=temp;

}

printf（"MatrixQ:

\n"）;

for（i=0;i

{

for（j=0;j

printf（"%f\t",Q（i,j））;

printf（"\n"）;

}

}

time2=MPI_Wtime（）;

if（my_rank==0）

{

printf（"\n"）;

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

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

printf（"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.0000002.0000001.000000

实验结果：

1、编译成功：

2、运行成果：

个人总结：

矩阵QR分解是研究生“随机过程”中较为重要的一个内容，试验中将涉及到的学生知识转化成程序有一定的难度。

试验中遇到了许多问题，由于这方面的经验比较缺少，特别是并行实验平台上进行编程，更是一大挑战。

所以摸索的过程有点艰苦。

通过此次实验，对并行编程方法有了一定的了解。

希望能够通过今后的学习来取得大的进步，对KD60实验平台有更深入的了解。

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置：首页 > 人文社科 > 法律资料

copyright@ 2008-2023 冰点文库 网站版权所有

经营许可证编号:鄂ICP备19020893号-2