QR分解-實(shí)驗(yàn)報(bào)告(共7頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上并行計(jì)算課程考核 實(shí)驗(yàn)報(bào)告考核題目：QR分解算法的并行實(shí)現(xiàn)并行實(shí)現(xiàn)方式：OpenMP1. 問(wèn)題概述QR分解法是目前求一般矩陣全部特征值的最有效并廣泛應(yīng)用的方法，一般矩陣先經(jīng)過(guò)正交相似變化成為householder矩陣，然后再應(yīng)用QR方法求特征值和特征向量。它是將矩陣分解成一個(gè)正規(guī)正交矩陣Q與上三角形矩陣R，所以稱(chēng)為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(串行時(shí)間time=%.0f ms,double(end - start);system(pause);3. 并行化設(shè)計(jì)思路a. 將有for循環(huán)且沒(méi)有數(shù)據(jù)關(guān)系、計(jì)算量大的計(jì)算進(jìn)行并行優(yōu)化；b. 將毫不相關(guān)且可獨(dú)

6、立運(yùn)行的代碼塊進(jìn)行并行優(yōu)化;c. 程序主要優(yōu)化是對(duì)householder變換里的乘法和矩陣賦值進(jìn)行優(yōu)化。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(串行結(jié)果nR: n);Matrix_print(A);printf(Q: n);Matrix_print(Q);printf(n并行時(shí)間：time = %.0f ms,time);system(pause);return 0;5. 實(shí)驗(yàn)結(jié)果分析1) 測(cè)試平臺(tái)描述系統(tǒng)：Windows 7旗艦版編譯器：Visual Studio 2012內(nèi)存：4GCPU：AMD A6-4400M APU with Radeon(tm) HD