矩阵QR分解并行实现.docx
- 文档编号:1662295
- 上传时间:2023-05-01
- 格式:DOCX
- 页数:16
- 大小:67KB
矩阵QR分解并行实现.docx
《矩阵QR分解并行实现.docx》由会员分享,可在线阅读,更多相关《矩阵QR分解并行实现.docx(16页珍藏版)》请在冰点文库上搜索。
矩阵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实验平台有更深入的了解。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 QR 分解 并行 实现