矩阵QR分解并行实现

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

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

实验原理：设A是n阶非奇异矩阵，则存在正交（酉）矩阵Q与实（复）非奇异 上三角矩阵R使得 A=Q*R，且除去相差一个对角元素的绝对值（模）全为1的对角因子外，上述分解唯一。

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

源代码：ch_qr.c #include \#include \#include \#include \

#define a(x,y) a[x*M+y] #define q(x,y) q[x*M+y] #define A(x,y) A[x*M+y] #define Q(x,y) Q[x*M+y] #define R(x,y) R[x*M+y]

float temp; float *A; float *R; float *Q;

double starttime; double time1; double time2; int p;

MPI_Status status;

void Environment_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); }

int main(int argc, char **argv) {

int M,N,m; int z;

int i,j,k,my_rank,group_size; float *ai,*qi,*aj,*qj; float c,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(\ fscanf(fdA,\

if(M != N) {

puts(\ 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 < M; j ++) fscanf(fdA, \

}

fclose(fdA);

for(i=0; i

for(j=0; j

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(\}

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); }