QR分解-实验报告(共7页)

1、精选优质文档-倾情为你奉上并行计算课程考核 实验报告考核题目：QR分解算法的并行实现并行实现方式：OpenMP1. 问题概述QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为householder矩阵，然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法。2. 串行代码描述串行主要代码如下：#include #include #include #include #define n 100 /maxNvoid Matrix_print(double Ann)for (int i=0;in;i+)fo

2、r (int j=0;jn;j+)printf(%8.4lft,Aij);printf(n);double Matrix_norm(double an) double d=0;for (int i=0;in;i+)d+=ai*ai;return sqrt(d);void Matrix_multiply(double Ann,double Bnn,double Cnn)for (int i=0;in;i+)for (int j=0;jn;j+)Cij=0;for (int t=0;tn;t+)Cij+=Ait*Btj;void Matrix_copy(double Ann,double Bnn)

3、for(int i=0;in;i+)for(int j=0;jn;j+)Aij=Bij;void householder_trans(double Ann,int k,double Qnn)double an;for(int i=0;in-k;i+)ai=0;for(int i=n-k;in;i+)ai=Ain-k; an-k-=Matrix_norm(a);double d=Matrix_norm(a);for(int i=0;in;i+)ai=ai/d;double Hnn;for(int i=0;in;i+)for(int j=0;jn;j+)Hij=-2*ai*aj;Hii+; /?H

4、=I-2vvTdouble tempnn;Matrix_multiply(H,A,temp);Matrix_copy(A,temp);Matrix_multiply(Q,H,temp);Matrix_copy(Q,temp);void Matrix_input(double Ann) srand( (unsigned)time( NULL ) ); for(int i=0 ;in;i+)for(int j=0;jn;j+)/*printf(a%d%d=,i,j);scanf(%lf,&Aij);*/Aij = rand();void main()double Qnn;double Ann;Ma

5、trix_input(A);printf(A: n);Matrix_print(A);for (int i=0;in;i+)for (int j=0;j=2;i-)householder_trans(A,i,Q);clock_t end = clock();printf(nnR: n);Matrix_print(A);printf(Q: n);Matrix_print(Q);printf(串行时间time=%.0f ms,double(end - start);system(pause);3. 并行化设计思路a. 将有for循环且没有数据关系、计算量大的计算进行并行优化；b. 将毫不相关且可独

6、立运行的代码块进行并行优化;c. 程序主要优化是对householder变换里的乘法和矩阵赋值进行优化。4. 并行代码描述#include #include #include #include #include #define n 4 /maxNvoid Matrix_print(double Ann)for (int i=0;in;i+)for (int j=0;jn;j+)printf(%8.4lft,Aij);printf(n);void Matrix_input(double Ann)srand(unsigned(time(NULL);for(int i=0 ;in;i+)for(in

7、t j=0;jn;j+)Aij = rand();/*for(int i=0 ;in;i+)for(int j=0;jn;j+)printf(a%d%d=,i,j);scanf(%f,&aij);*/double Matrix_norm(double an) double d=0;#pragma omp parallel for num_threads(2) reduction(+:d)for (int i=0;in;i+)d+=ai*ai;return sqrt(d);void Matrix_multiply(double Ann,double Bnn,double Cnn)int bord

8、er = n*n;int i,j;#pragma omp parallel for num_threads(2)for (int index = 0; index border; index+)i = index/n;j = index%n;Cij=0;for (int t=0;tn;t+)Cij+=Ait*Btj;void Matrix_copy(double Ann,double Bnn)int i,j;int border = n*n;#pragma omp parallel for num_threads(2)for (int index = 0; index border; inde

9、x+)i = index/n;j = index%n;Aij=Bij;void householder_trans(double Ann,int k,double Qnn)double an;#pragma omp parallel num_threads(2)#pragma omp forfor(int i=0;in-k;i+)ai=0;#pragma omp forfor(int i=n-k;in;i+)ai=Ain-k;an-k-=Matrix_norm(a);double d=Matrix_norm(a);#pragma omp parallel for num_threads(2)f

10、or(int i=0;in;i+)ai=ai/d;double Hnn;int border = n*n;int i,j;#pragma omp parallel for num_threads(2)for (int index = 0; index border; index +)i = index/n;j = index%n;Hij = -2*ai*aj;if (j = n-1)Hii+;double tempnn;#pragma omp parallel sections private(temp) num_threads(2)#pragma omp sectionMatrix_mult

11、iply(H,A,temp);Matrix_copy(A,temp);#pragma omp sectionMatrix_multiply(Q,H,temp);Matrix_copy(Q,temp);int main()double Qnn;double Ann;Matrix_input(A);/printf(A: n);/Matrix_print(A);int border = n*n;int i,j;#pragma omp parallel for num_threads(2)for (int index = 0; index =2;i-)householder_trans(A,i,Q);clock_t end = clock();double time = (double)end - start;printf(串行结果nR: n);Matrix_print(A);printf(Q: n);Matrix_print(Q);printf(n并行时间：time = %.0f ms,time);system(pause);return 0;5. 实验结果分析1) 测试平台描述系统：Windows 7旗舰版编译器：Visual Studio 2012内存：4GCPU：AMD A6-4400M APU with Radeon(tm) HD