当前位置:首页 >> 建筑/土木 >>

QR分解 实验报告


并行计算课程考核 实验报告
考核题目:QR 分解算法的并行实现 并行实现方式:OpenMP 1. 问题概述 QR 分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交 相似变化成为 householder 矩阵,然后再应用 QR 方法求特征值和特征向量。它是将矩阵分 解成一个正规正交矩阵 Q 与上三角形矩阵 R,所以称为 QR 分解法。 2. 串

行代码描述 串行主要代码如下: #include <stdio.h> #include <math.h> #include <stdlib.h> #include <time.h> #define n 100 //maxN void Matrix_print(double A[n][n]) { for (int i=0;i<n;i++) { for (int j=0;j<n;j++) printf("%8.4lf\t",A[i][j]); printf("\n"); } } double Matrix_norm(double a[n]) { double d=0; for (int i=0;i<n;i++) d+=a[i]*a[i]; return sqrt(d); } void Matrix_multiply(double A[n][n],double B[n][n],double C[n][n]) { for (int i=0;i<n;i++) for (int j=0;j<n;j++) { C[i][j]=0; for (int t=0;t<n;t++) C[i][j]+=A[i][t]*B[t][j]; }

} void Matrix_copy(double A[n][n],double B[n][n]) { for(int i=0;i<n;i++) for(int j=0;j<n;j++) A[i][j]=B[i][j]; } void householder_trans(double A[n][n],int k,double Q[n][n]) { double a[n]; for(int i=0;i<n-k;i++) a[i]=0; for(int i=n-k;i<n;i++) a[i]=A[i][n-k]; a[n-k]-=Matrix_norm(a); double d=Matrix_norm(a); for(int i=0;i<n;i++) a[i]=a[i]/d; double H[n][n]; for(int i=0;i<n;i++) { for(int j=0;j<n;j++) H[i][j]=-2*a[i]*a[j]; H[i][i]++; } //μ?H=I-2vvT double temp[n][n]; Matrix_multiply(H,A,temp); Matrix_copy(A,temp); Matrix_multiply(Q,H,temp); Matrix_copy(Q,temp); } void Matrix_input(double A[n][n]) { srand( (unsigned)time( NULL ) ); for(int i=0 ;i<n;i++) for(int j=0;j<n;j++) { /*printf("a[%d][%d]=",i,j); scanf("%lf",&A[i][j]);*/ A[i][j] = rand(); } } void main()

{ double Q[n][n]; double A[n][n]; Matrix_input(A); printf("A: \n"); Matrix_print(A); for (int i=0;i<n;i++) { for (int j=0;j<n;j++) Q[i][j]=0; Q[i][i]=1; } clock_t start = clock(); for (int i=n;i>=2;i--) householder_trans(A,i,Q); clock_t end = clock(); printf("\n\nR: \n"); Matrix_print(A); printf("Q: \n"); Matrix_print(Q); printf("串行时间 time=%.0f ms",double(end - start)); system("pause"); } 3. 并行化设计思路 a. 将有 for 循环且没有数据关系、计算量大的计算进行并行优化; b. 将毫不相关且可独立运行的代码块进行并行优化; c. 程序主要优化是对 householder 变换里的乘法和矩阵赋值进行优化。 4. 并行代码描述 #include <stdio.h> #include <math.h> #include <stdlib.h> #include <omp.h> #include <time.h> #define n 4 //maxN void Matrix_print(double A[n][n]) { for (int i=0;i<n;i++) { for (int j=0;j<n;j++) printf("%8.4lf\t",A[i][j]); printf("\n"); } }

void Matrix_input(double A[n][n]) { srand(unsigned(time(NULL))); for(int i=0 ;i<n;i++) for(int j=0;j<n;j++) { A[i][j] = rand(); } /*for(int i=0 ;i<n;i++) for(int j=0;j<n;j++) { printf("a[%d][%d]=",i,j); scanf("%f",&a[i][j]); }*/ } double Matrix_norm(double a[n]) { double d=0; #pragma omp parallel for num_threads(2) reduction(+:d) for (int i=0;i<n;i++) d+=a[i]*a[i]; return sqrt(d); } void Matrix_multiply(double A[n][n],double B[n][n],double C[n][n]) { 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; C[i][j]=0; for (int t=0;t<n;t++) C[i][j]+=A[i][t]*B[t][j]; } } void Matrix_copy(double A[n][n],double B[n][n]) { int i,j; int border = n*n; #pragma omp parallel for num_threads(2)

for (int index = 0; index < border; index++) { i = index/n; j = index%n; A[i][j]=B[i][j]; } } void householder_trans(double A[n][n],int k,double Q[n][n]) { double a[n]; #pragma omp parallel num_threads(2) { #pragma omp for for(int i=0;i<n-k;i++) a[i]=0; #pragma omp for for(int i=n-k;i<n;i++) a[i]=A[i][n-k]; } a[n-k]-=Matrix_norm(a); double d=Matrix_norm(a); #pragma omp parallel for num_threads(2) for(int i=0;i<n;i++) a[i]=a[i]/d; double H[n][n]; 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; H[i][j] = -2*a[i]*a[j]; if (j == n-1) { H[i][i]++; } } double temp[n][n]; #pragma omp parallel sections private(temp) num_threads(2) { #pragma omp section

{ Matrix_multiply(H,A,temp); Matrix_copy(A,temp); } #pragma omp section { Matrix_multiply(Q,H,temp); Matrix_copy(Q,temp); } } } int main() { double Q[n][n]; double A[n][n]; 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 < border; index ++) { i = index/n; j = index%n; Q[i][j] = 0; if (j == n-1) { Q[i][i] = 1; } } clock_t start = clock(); for (int i=n;i>=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 内存:4G CPU:AMD A6-4400M APU with Radeon(tm) HD 双核 2) 测试结果描述 (串行运行时间、并行运行时间) 数据量 80 100 120 140 160 串行时间(ms) 562 1310 2543 4992 9329 并行时间(ms) 515 1186 2153 3978 7285 加速比 1.091 1.105 1.181 1.255 1.281

3) 加速比分析 加速比=1/5(1.091+1.105+1.181+1.255+1.281) =1.1826


相关文章:
实验三 矩阵的QR分解
实验三 1、原理 矩阵的 QR 分解 矩阵 QR 分解是一种特殊的三角分解,在解决矩阵特征值的计算、最小二乘法等问题中 起到重要作用。 ×设 A∈Cm n, m≥n ...
QR分解的数值效果报告
§ 7.参考文献与附注……….20 2 QR 分解数值效果分析报告基于数值试验的 gram-schmit 方法,modified gram-schmit 方法,householder 变换,givens 变换计 算矩阵 Q...
QR分解及其应用
QR分解及其应用_数学_自然科学_专业资料。矩阵QR分解及其应用 《矩阵分析与应用》专题报告 ——QR 分解及应用 学生姓名:卢楠、胡河群、朱浩 2015 年 11 月 25 ...
数值代数qr方法实验报告
数值代数qr方法实验报告_数学_自然科学_专业资料。数值代数有关qr方法的实验报告...b ) 2 ,为方便求解,需要 Q T A 是上三角矩阵,这样引入 QR 分解就比较求解...
LU和QR分解法解线性方程组
LU和QR分解法解线性方程组_IT/计算机_专业资料。介绍LU和QR分解法解线性方程组的原理,步骤,上机代码。实验报告LU 和 QR 法解线性方程组 一、 问题描述 ?4 ?...
QR分解解方程的MATLAB程序
QR分解解方程的MATLAB程序_理学_高等教育_教育专区。function [Q,R,x]=QR(A,b) n=length(A);Q=eye(n);R=A; for i=1:n-1 x=R(i:n,i); if ...
数值分析大作业QR分解
[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];dp_qrfenjie(A); 得到用带位移的 QR 分解法所用的迭代次数 i,舒尔矩阵 Ak,以及...
QR分解、SVD分解、最小二乘问题数值上机报告
QR 分解、SVD 分解、最小二乘问题数值上机报告 练习 6.13 随机构造一个可逆...实验步骤: 1、首先分别编写出 CGS、MGS、Householder、Givens 等不同方法进行...
浅谈矩阵的LU分解和QR分解及其应用
浅谈矩阵的LU分解和QR分解及其应用_数学_自然科学_专业资料。浅谈矩阵的 LU 分解和 QR 分解及其应用基于理论研究和计算的需要, 往往有必要把矩阵分解为具有某种特性...
数值分析实验报告——基本QR算法求全部特征值
数值分析实验报告专业 信息与计算科学 班级 信计 101 姓名 学号 协作队员 实验...用基本 QR 算法求全部特征值(可用 Matlab 函数“qr”实现 矩阵的 QR 分解) ...
更多相关标签:
过氧化氢分解实验报告 | 分解电压测定实验报告 | 铵盐的热分解实验报告 | qr分解 | 矩阵qr分解 | matlab qr分解 | 矩阵的qr分解 | householder qr分解 |