数值分析大作业QR分解
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
数值分析大作业
第二次计算实验：SVD及其应用梁杰存2014310739航博1431.方法求矩阵A奇异值分解一个途径是求解A T A的特征值，但因为舍入误差容易丢掉小奇异值。
因此通常先将矩阵上双对角化，即构造正交阵Q和W,使得Q T AW=B(upper−bidiagnal)。
这一过程可以通过逐次Household变换或逐次Given’s变换完成，还有一种基于待定系数法思想的Lanczos算法。
由于Linpack中SVD算法需要输入上双对角矩阵，本文采用Lanczos 算法实现上双对角化。
1.1.隐式零移位QR法（implicit zero-shift QR）与传统移位QR迭代算法不同，隐式零移位QR算法不进行移位，并且第一步构造右乘Given’s变换矩阵GR(1,2)将上双对角矩阵B的(1,2)位置上的元素12b消零，而不是传统方法中引入一个非零元素21b。
但这一步可能会使原来为零的b12变为非零。
第二步左乘Given’s阵GL(1,2)使得12b为0，但可能会使为零b13变为非零。
与上述步骤类似，将b13变为0后可能会使b23非零。
如下图所示，重复上述步骤最终将恢复为上双对角矩阵，即完成一步隐式零移位QR迭代。
反复迭代，矩阵B将趋近与对角阵阵，对角元即特征值。
图1隐式QR迭代1.2.分而治之(Divide-and-conquer)分而治之算法将上双对角阵B分成有两个互相独立对角块矩阵与另一矩阵之和，即：B=B100B2+b m vvT=Q1Σ1Q1T00Q2Σ1Q2T+b m vv T =Q100Q2(Σ100Σ1+b m uuT)Q1T00Q2T所以矩阵B的特征值与矩阵D+ρu u T的特征值相同，其中D=Σ100Σ1为对角阵，又：det D+ρu u T−λI=det((D−λI)(I+ρD−λ−1u u T))由于D−λI非奇异，则det I+ρD−λ−1u u T=1+ρu T D−λ−1u=1+ρu i2d i−λ=0ni=1在每个d i与d i+1之间分布着一个特征值，可用牛顿法快速找到该特征值。
北航数值分析第二次大作业--QR分解
《数值分析A》计算实习题目二姓名学号联系方式班级指导教师2012年10月一、算法设计方案整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。
因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。
1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。
2.对经过拟上三角化处理的矩阵进行QR分解。
3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。
计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。
二、源程序代码#include<stdio.h>#include<math.h>#include<string.h>int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s;double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10];double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12;void main(){void Quasiuppertriangular(double A[][10]);void QRdecomposition(double A[][10]);void DoublestepsQR(double A[][10]);int i,j;for(i=0;i<10;i++){for(j=0;j<10;j++){A[i][j]=sin(0.5*(i+1)+0.2*(j+1));Q[i][j]=0;AA[i][j]=A[i][j];}A[i][i]=1.5*cos(2.2*(i+1));AA[i][i]=A[i][i];}Quasiuppertriangular(A); //调用拟上三角化函数printf( "\n A经过拟上三角化矩阵为：\n\n");for(i=0;i<10;i++) //输出拟上三角化矩阵{for(j=0;j<10;j++){printf("%.12e ",A[i][j]); //输出拟上三角化矩阵}printf( "\n\n");}QRdecomposition(A); //调用QR分解函数printf( " 进行QR分解后，R矩阵为：\n\n"); //输出R矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",R[i][j]);}printf( "\n\n");}printf( " Q矩阵为:\n\n"); //输出Q矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",Q[i][j]);}printf( "\n\n");}printf( " RQ矩阵为：\n\n"); //输出RQ矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",RQ[i][j]);}printf( "\n\n");}DoublestepsQR(A); //调用双步位移函数printf( "\n\n 特征值实部依次为：\n\n"); //输出特征值实部for(j=0;j<10;j++){printf("%.12e ",Re[j]);}printf("\n\n 特征值虚部依次为：\n\n "); //输出特征值虚部for(j=0;j<10;j++){printf("%.12e ",Im[j]);}//按行输出特征向量printf( "\n\n 按行输出实特征根相应特征向量为：\n\n");for(i=0;i<10;i++){if(i==1||i==2||i==5||i==6){continue;}for(j=0;j<10;j++){printf("%.12e ",X[i][j]);}printf( "\n\n");}getchar();}//拟上三角化函数void Quasiuppertriangular(double A[][10]) {for(j=0;j<8;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;T[i]=0;W[i]=0;}m=0;for(i=j+2;i<10;i++){if(A[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j+1;i<10;i++){d=d+pow(A[i][j],2);}d=sqrt(d);c=-d;if(A[j+1][j]<=0){c=d;}h=c*(c-A[j+1][j]);U[j+1]=A[j+1][j]-c;for(i=j+2;i<10;i++){U[i]=A[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*A[k][i];}P[i]=P[i]/h;}t=0;for(i=0;i<10;i++){for(k=0;k<10;k++){T[i]=T[i]+U[k]*A[i][k];}T[i]=T[i]/h;t=t+P[i]*U[i];}t=t/h;for(i=0;i<10;i++){W[i]=T[i]-t*U[i];for(k=0;k<10;k++){A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];if(abs(A[i][k])<1e-12){A[i][k]=0;}}}}}//QR分解函数void QRdecomposition(double A[][10]) {for(i=0;i<10;i++){for(j=0;j<10;j++){RQ[i][j]=0;Q[i][j]=0;R[i][j]=A[i][j];}Q[i][i]=1;}for(j=0;j<9;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;W[i]=0;}m=0;for(i=j+1;i<10;i++){if(R[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j;i<10;i++){d=d+pow(R[i][j],2);}d=sqrt(d);c=-d;if(R[j][j]<=0){c=d;}h=c*(c-R[j][j]);U[j]=R[j][j]-c;for(i=j+1;i<10;i++){U[i]=R[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){W[i]=W[i]+U[k]*Q[i][k];}}for(i=0;i<10;i++){for(k=0;k<10;k++){Q[i][k]=Q[i][k]-((W[i]*U[k])/h);}}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*R[k][i];}P[i]=P[i]/h;}for(i=0;i<10;i++){for(k=0;k<10;k++){R[i][k]=R[i][k]-U[i]*P[k];if(abs(R[i][k])<epsilon){R[i][k]=0;}}}}for(i=0;i<10;i++) //计算A(n+1)=RQ {for(j=0;j<10;j++){for(k=0;k<10;k++){RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];}}}}//双步位移法计算特征值特征向量函数void DoublestepsQR(double A[][10]){int L=1000,m=9; //定义最大循环次数for(i=0;i<L;i++){for(;m>-1;){if(abs(A[m][m-1])<=epsilon){Re[m]=A[m][m];m=m-1; //降阶if(m==0) //4{Re[0]=A[0][0];break;}if(m==-1){break;}if(m>1){continue;}}b=-A[m][m]-A[m-1][m-1]; //5c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];if(m==1) //6{if((b*b-4*c)>=0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-1; //循环出口条件break;}if((m>1)&&(abs(A[m-1][m-2])>epsilon)) //8{if(i==L-1){printf("No results! \n");m=0; //循环出口条件break;}break;}if((m>1)&&(abs(A[m-1][m-2])<=epsilon)) //7 {if((b*b-4*c)>0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-2; //降阶if(m>0){continue;}if(m==0){Re[0]=A[0][0];break;}}}if(m<=0){break;}s=A[m-1][m-1]+A[m][m]; //9t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=0;Q[j][k]=0;M[j][k]=0;X[j][k]=0;Y[j][k]=0;}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){M[j][k]=M[j][k]+A[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){M[j][k]=M[j][k]-s*A[j][k];}M[j][j]=M[j][j]+t;}//调用QR分解函数对M矩阵进行分解并传递参数矩阵QQRdecomposition(M);for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=Q[k][j];}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){X[j][k]=X[j][k]+Qt[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];}}}for(j=0;j<10;j++){{A[j][k]=Y[j][k];}}}//应用列主元高斯消元法计算实部特征向量for(l=0;l<10;l++){if(l==1||l==2||l==5||l==6){continue;}for(k=0;k<10;k++){for(m=0;m<10;m++){A[k][m]=AA[k][m];}A[k][k]=A[k][k]-Re[l];}for(j=0;j<9;j++){m=j;for(i=j+1;i<10;i++){if(abs(A[i][j])>abs(A[m][j])){m=i;}}{Y[j][k]=A[j][k];A[j][k]=A[m][k];A[m][k]=Y[j][k];}for(k=j+1;k<10;k++){b=A[k][j]/A[j][j];for(i=j;i<10;i++){A[k][i]=A[k][i]-A[j][i]*b;}}}X[l][9]=1;for(i=8;i>=0;i--){c=0;for(j=i+1;j<10;j++){c=c+A[i][j]*X[l][j];}X[l][i]=-c/A[i][i];}}}三、程序输出结果1819。
数值分析7.2矩阵的正交分解与求矩阵全部特征值的QR方法
但 x y ,则存在householder阵
2
2
UU T
H I2 U 2
2
使Hx y，其中U x y。 W
x
x y
y
证：若设W U ，则有 W 1，因此
U
2
I
H 2
2
I 2WW T
(x y) x y 2
( xT
yT
I )
UU T 2 U2
2
Hx
x
2
( x2 y) x y 2
2
k
)
,
H
k
(k 2
)
,
,
H
k
(k n
)
A(k 1)
1(
k
1)
,
(k 2
1)
,
,
(k n
1)
a(2) 11 0
a(2) 1k
Hk
A(k )
Hk
0
a(k) kk
0
0
a(k) nk
a(2) 1n
a(k kn
)
a(k) nn
H
(
k1
k
)
,
H
k
(k 2
)
,
,
H
k
(k n
)
a1(12) 0 0 0
迭代格式
Ak Qk Rk Ak 1 RkQk
(k 1, 2, ).
将A A1化成相似的上三角阵（或分块上三角阵），
从而求出矩阵A的全部特征值与特征向量。
由A A1 Q1R1 ,即Q11 A R1。 于是A2 R1Q1 Q1 AQ1 ,即A2与A相似。
同理可得，Ak A (k 2, 3, )。 故它们有相同的特征值。
北航数值分析大作业二
数值分析第二题1 算法设计方案要想得出该题的答案首先要将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解。
要得出矩阵A 的全部特征值先对A 进行QR 的双步位移得出特征值。
采用列主元的高斯消元法求解特征向量。
1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。
具体算法如下所示，记A A =)1(，并记)(r A 的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。
对于2,,2,1-=n r 执行1. 若()n r r i a r ir,,3,2)( ++=全为零，则令)()1(r r A A =+,转5；否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()n Tr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。
4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωTrr T r r r r p u u A A --=+ω)()1(5. 继续。
1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r 执行1.若()n r r i a r ir,,3,1)( ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。
2.计算()∑==nri r irr a d 2)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,2r r r r r r a c c h -=3.令()n Tr nrr r r r r r r r R a a c a u ∈-=+)()(,2)(,,,,,0,,0 。
数值分析大作业之2
数值分析大作业二一、算法设计方案(一)算法流程1.将要求解的矩阵进行Householder变换，进行拟上三角化，并输出拟上三角化的结果；2.将拟上三角化后的矩阵进行带双步位移的QR分解(最大迭代次数1000次)，逐个求出特征值，输出QR法结束后得到的Q、R、RQ阵，输出求出的特征值（用实部和虚部表示）3.对于求出来的实特征值，使用带原点平移的反幂法求出对应的特征向量，输出这些特征向量。
（二）程序设计流程1. 定义精度和最大迭代次数；2. 使用容器vector进行编程（方便增减元素），使用传引用调用（提高执行效率）；3. 将各个步骤用到的数学算法封装成函数，方便调用。
具体需要的函数如下：double VectMod(vector< double > &p)：求向量的模vector< double > NumbMultiVect(vector< double > &vect, double a)：数乘向量double VectMultiVect(vector< double > &y, vector< double > &u)：求两个向量的积vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector<double > &u)：矩阵的转置乘向量vector< double > ArraMultiVect(vector< vector< double > > &B, vector< double > &u)：矩阵乘向量vector< vector< double > > VectMultiCovVect(vector< double > &a, vector< double >&b)：向量乘向量转置得到矩阵void ArraSubs(vector< vector< double > > &C, vector< vector < double > > D) ：两个矩阵相减Vector< double > VectSubs(vector< double > &a, vector< double > &b)：两个向量相减vector<vector<double>> GetM(vector<vector<double>> &A)：求解矩阵M （用于双部位移QR迭代用）double GetMode(vector < vector < double > > &B, const int r)：求解矩阵B的r列向量的模double GetModeH(vector < vector < double > > &B, const int r)：求解矩阵B的第r列向量的模，用于拟对角化vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)：一个实数乘矩阵bool IsBirZeroH(vector< vector< double > > &B, const int r)：判断B[i][r]对角线下是否为零void GausElim(vector< vector< double > > a)：列主元高斯消元法求齐次方程解向量void Stop(vector< vector< double > > &Ar)：停止，结束程序void SolutS1S2(complex< double > &s1, complex< double > &s2, vector< vector<double > > &A)：求解二阶子阵的特征值s1，s2；void Save2(complex< double > &s1, complex< double > &s2)：保存两个特征值void Save1(complex< double > &s)：保存一个特征值void JudgemBelow2(vector< vector< double > > &A, vector< vector< double > > Abk)：对于m == 1 及m == 0 的处理void Hessenberg(vector< vector< double > > &A)：矩阵拟上三角化void QRMethod(vector< vector< double > > A)：矩阵QR分解void CalculatAk(vector< vector < double > > &Ak)：带双步位移QR迭代法二、源程序#include "stdafx.h"#include <vector>#include <iostream>#include <math.h>#include <complex>#include <fstream>using namespace std;const double epsion = 1e-12;const int L = 1000;int m,n;int k = 1;vector< vector< double > > I;vector< complex< double > > Lambda;///////////////////////////以下为自定义的算法流程中用到的函数double VectMod(vector< double > &p) //求向量的模{double value = 0.0;vector< double >::size_type i,j;j = p.size();for (i=1; i<j; i++){value += p[i] * p[i];}value = sqrt(value);return value;}vector< double > NumbMultiVect(vector< double > &vect, double a) //数乘向量{int j = vect.size();vector< double > b(j, 0);for (int i=1; i<j; i++){b[i] = a * vect[i];}return b;}double VectMultiVect(vector< double > &y, vector< double > &u)//两个向量相乘{vector< double >::size_type a = y.size();double value = 0;for (vector< double >::size_type i=1; i<a; i++){value += y[i] * u[i];}return value;}vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector< double > &u)//矩阵的转置乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[j][i] * u[j];}}return vec;}}vector< double > ArraMultiVect(vector< vector< double > > &B,vector< double > &u) //矩阵乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[i][j] * u[j];}}return vec;}}vector< vector<double> > GetM(vector< vector <double> > &A){int a = A.size();double s = A[a-2][a-2] + A[a-1][a-1];double t = A[a-2][a-2] * A[a-1][a-1] - A[a-1][a-2] * A[a-2][a-1];vector<vector<double>> D(a, vector< double >(a, 0));for (int i=1; i<a; i++){{double sum = 0;for (int k=1; k<a; k++){sum += A[i][k] * A[k][j];}D[i][j] = sum - s * A[i][j] + t *(i==j ? 1.0 : 0);}}return D;}double GetMode(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}double GetModeH(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r+1; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)//数乘//向量{int b = D.size();vector< vector< double > > U(b, vector< double >(b, 0));for (int i=1; i<b; i++){{U[i][j] = a *D[i][j];}}return U;}bool IsBirZero(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+1; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}bool IsBirZeroH(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+2; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}vector< vector< double > > VectMultiCovVect(vector< double > &a,vector< double > &b)//向量乘向量的转置得到矩阵{int s1 = a.size();int s2 = b.size();if (s1 != s2){cerr << "Vectors not match in size ! ";}else{vector< vector< double > > U(s1, vector< double >(s1, 0));for (int i=1; i<s1; i++){for (int j=1; j<s1; j++){U[i][j] = a[i] * b[j];}}return U;}}void ArraSubs(vector< vector< double > > &C,vector< vector < double > > D){int a = C.size();int b = D.size();int c = C[0].size();int d = D[0].size();if (a!=b || c!=d){cerr << "Vectors not match in size !";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){C[i][j] -= D[i][j];}}}}vector< double > VectSubs(vector< double > &a,vector< double > &b){int s1 = a.size();int s2 = b.size();if(s1 != s2){cerr << "Vectors not match in size !";}else{vector< double > value(s1,0);for (int i=1;i<s1;i++){value[i] = a[i] - b[i];}return value;}}void GausElim(vector< vector< double > > a) //高斯消元{int sz = a.size();vector< int > fx, ufx,ufxp, record;vector< double > x(sz, 0);vector< vector< double > > ret;//*****消元过程********for (int k = 1; k < sz-1; k++){double max = a[k][k];int p=0;for (int i=k+1; i<sz; i++){if (abs(a[i][k]) > abs(max)){p =i;max = a[i][k];}}//选出主元行if (abs(max) >= epsion){if (p != 0){for (int j=k; j<sz; j++){double temp = a[k][j];a[k][j] = a[p][j];a[p][j] = temp;}//交换主元行}for (int i=k+1; i<sz; i++){double tt = a[i][k] / a[k][k];for (int j=k+1; j<sz; j++){a[i][j] = a[i][j] - tt * a[k][j];}}}// end of if (abs(max) >= epsion)}// end of for (int k = 1; k < sz; k++)for (int i=1; i<sz; i++){int p=0;for (int j=1; j<=i;j++){if(abs(a[j][i]) >= 10*epsion)//不为零{p=j;}}record.push_back(p);}for (int i=1; i<sz; i++){int p=0;for (int j=sz-2; j>=0; j--){if (record[j] == i){p=j + 1;}}if (p != 0){ufxp.push_back(p);ufx.push_back(i);}}for (int i=1; i<sz; i++){int p = 0;for (int j=0; j < ufxp.size(); j++){if (ufxp[j] == i){p = 1;}}if (p == 0){fx.push_back(i);}}//end of for (int i=1; i<sz; i++)int c = fx.size();//***************************************************** for (int i=0; i<c; i++){for (int j=0; j<c; j++){if (i == j){x[fx.at(j)] = 1.0 + epsion;}else{x[fx.at(j)] = 0;}}int b = ufxp.size();for (int s=b-1; s>=0; s--){double temp=0;for(int t=ufxp.at(s)+1; t<sz; t++){temp += x[t] * a[ufx.at(s)][t];}x[ufxp.at(s)] = (0 - temp) / a[ufx.at(s)][ufxp.at(s)];}ret.push_back(x);}//end of for (int i=0; i<c; i++)//******************************************************* int sz1 = ret.size();int sz2 = ret[0].size();ofstream result2("CharactVector .txt", ios::app);result2.precision(12);for (int i=0; i<sz1; i++){for (int j=1;j<sz2; j++){result2 << scientific << ret[i][j] << endl;}result2 << endl << endl;}result2 << "下一个特征值的特征向量！" << endl;}void Stop(vector< vector< double > > &Ar){int a = Lambda.size();for (int i=0; i<a; i++){if (abs(Lambda[i].imag())<=epsion){vector< vector< double > > An=Ar;ArraSubs(An,NumbMultiArra(I,Lambda[i].real()));GausElim(An);}}ofstream result("result.txt");result.precision(12);vector<complex< double > >::iterator i,j;j = Lambda.end();for (i = Lambda.begin(); i<j; i++){result << scientific <<(*i) << endl;}exit(0);}void SolutS1S2(complex< double > &s1, complex< double > &s2,vector< vector< double > > &A){double s = A[m-1][m-1] + A[m][m];double t = A[m-1][m-1] * A[m][m] - A[m][m-1] * A[m-1][m];double det = s *s - 4.0 * t;if (det > 0){s1.imag(0);s1.real ((s + sqrt(det)) / 2);s2.imag (0);s2.real((s - sqrt(det)) / 2);}else{s1.imag(sqrt(0 - det) / 2);s1.real(s / 2);s2.imag(0 - sqrt(0 - det) / 2);s2.real(s / 2);}}void Save2(complex< double > &s1, complex< double > &s2)//存储特征值s1，s2 {Lambda.push_back(s1);Lambda.push_back(s2);}void Save1(complex< double > &s){Lambda.push_back(s);}void JudgemBelow2(vector< vector< double > > &A,vector< vector< double > > Abk){if (m == 1){complex< double > lbd;lbd.imag(0);lbd.real(A[1][1]);Save1(lbd);Stop(Abk);}else if (m == 0){Stop(Abk);}}//拟上三角化void Hessenberg(vector< vector< double > > &A){int a = A.size();vector< double > u(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;for (int r=1; r<a-2; r++){if (!IsBirZeroH(A, r)){d = GetModeH(A,r);c = (A[r+1][r] > 0) ?(0-d) : d;h = c * c - c * A[r+1][r];for (int j=1; j<a; j++){if (j < r+1){u[j] = 0;}else if (j == r+1){u[j] = A[j][r] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)p = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);q = NumbMultiVect(ArraMultiVect(A, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(A, VectMultiCovVect(w, u));ArraSubs(A, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("NISHANGDANJIAO.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}}void QRMethod(vector< vector< double > > A){int a = A.size();vector< vector< double > > Q(a,vector < double > (a,0));for (int i=1; i<a; i++){Q[i][i] = 1.0;}vector< vector< double > > RQ(A);//vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);vector< double > l(a,0);double d, c, h, t;for (int r=1; r<a-1; r++){if (!IsBirZero(A, r)){d = GetMode(A,r);c = (A[r][r] > 0) ?(0-d) : d;h = c * c - c * A[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = A[j][j] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)l = ArraMultiVect(Q, u);ArraSubs(Q, VectMultiCovVect(l,NumbMultiVect(u, 1/h)));v = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);ArraSubs(A, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(RQ, u), 1 / h);q = NumbMultiVect(ArraMultiVect(RQ, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(RQ, VectMultiCovVect(w, u));ArraSubs(RQ, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("QR.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << Q[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << RQ[i][j] << " ";}result << endl;}}void CalculatAk(vector< vector < double > > &Ak){int a = Ak.size();vector< vector < double > > B(a,vector < double > (a,0));vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;B = GetM(Ak);for (int r=1; r<a-1; r++){if (!IsBirZero(B, r)){d = GetMode(B,r);c = (B[r][r] > 0) ?(0-d) : d;h = c * c - c * B[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = B[j][j] - c;}else{u[j] = B[j][r];}}// end of for (int j=1; j<=m; j++)v = NumbMultiVect(ConveArraMultiVect(B, u), 1 / h);ArraSubs(B, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(Ak, u), 1 / h);q = NumbMultiVect(ArraMultiVect(Ak, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(Ak, VectMultiCovVect(w, u));ArraSubs(Ak, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)}int _tmain(int argc, _TCHAR* argv[]){vector< vector< double > > A(11, vector< double >(11,0));vector< vector< double > > Abk;I=A;for (int i=1; i<11; i++){vector< double > temp;for (int j=1; j<11; j++){if(i != j){A[i][j] = sin(0.5 * i + 0.2 * j);I[i][j] = 0;}else{A[i][j] = 1.5 * cos(i + 1.2 *j);I[i][j] = 1.0;}}}Abk = A;n = A.size() - 1;m = n;//初始化问题Hessenberg(A);QRMethod(A);while (1){if (abs(A[m][m-1]) <= epsion){complex< double > lbdm;lbdm.imag(0);lbdm.real(A[m][m]);Save1(lbdm);m -= 1;//对A进行降维处理！！！！A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();}JudgemBelow2(A, Abk);}else{complex< double > va1, va2;SolutS1S2(va1, va2, A);if (m == 2){Save2(va1,va2);Stop(Abk);}//end of if (m == 2)if ( abs(A[m-1][m-2]) <= epsion){Save2(va1,va2);m = m - 2;//矩阵降维A.pop_back();A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();A[i].pop_back();}JudgemBelow2(A, Abk);}else{if (k == L){cerr << "Stop without solution";exit(-1);}else{CalculatAk(A);k += 1;}}}// end of if (abs(A[m][m-1]) >= epsion)}}三、实验结果（1）A经过拟上三角化后得到的矩阵-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000-7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000-8.782436382928e-002 -9.255889387184e-001 6.032599440534e-0011.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+0003.237804101546e-001 2.205798440320e-001 2.102692662546e+0001.816138086098e-001 1.278839089990e+000 -6.380578124405e-001-4.154075603804e-0011.0556********e-016 1.728274599968e+000 -1.171467642785e+000-1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+0002.924947206124e-001 -6.412830068395e-001 9.783997621285e-0022.557763574160e-001-5.393383812774e-017 -3.165873865181e-017 -1.291669534130e+000-1.111603513396e+000 1.171346824096e+000 -1.307356030021e+0001.803699177750e-001 -4.246385358369e-001 7.988955239304e-0021.608819928069e-0011.533464996622e-017 5.963321406181e-017 0.000000000000e+0001.560126298527e+000 8.125049397515e-001 4.421756832923e-001-3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0021.300562737229e-016 -3.097060010889e-017 0.000000000000e+0000.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000-3.042843176799e-001 2.528712446035e-001 -6.709925401449e-0012.544619929082e-0011.610216724767e-016 -2.211571837369e-016 0.000000000000e+0000.000000000000e+000 6.483933100712e-017 -7.463453456938e-001-2.708365157019e-002 -9.486521893682e-001 1.195871081495e-0011.929265617952e-0021.368550186199e-016 7.151513190805e-017 0.000000000000e+0000.000000000000e+000 -2.522454384246e-017 1.072074718358e-016-7.701801374364e-001 -4.697623990618e-001 4.988259468008e-0011.137691603776e-001-2.780851300718e-017 -6.708630788363e-017 0.000000000000e+0000.000000000000e+000 -3.282796821065e-017 4.879774287475e-0171.854829286220e-016 7.013167092107e-001 1.582180688475e-0013.862594614233e-001-2.124604440055e-017 -1.707979758930e-016 0.000000000000e+0000.000000000000e+000 1.013496750890e-016 -4.153758325241e-0171.222621691887e-016 -5.551115123126e-017 4.843807602783e-0013.992777995177e-001（2）Q-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-0016.486200753651e-002 3.709718861896e-001 1.855847143605e-001-1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-0026.124472142963e-003 4.440505443139e-002 1.975907909728e-0022.062836970533e-002-4.208697095111e-017 -8.810520554692e-001 -3.989762796959e-0013.720308728479e-002 2.127794064090e-001 1.064463557221e-001-9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002-2.150178169911e-017 4.009681353156e-017 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001 -3.533456368496e-0013.101438948264e-002 2.248676491598e-001 1.000601783527e-0011.044622748702e-0016.113458775639e-018 -3.721194648197e-0177.068175586628e-0189.892235468621e-001 -1.239411731211e-001 -6.200358589825e-0025.442272839461e-003 3.945881637235e-002 1.755813350011e-0021.833059462907e-0025.184948268593e-017 -4.198303559531e-017 3.255071298412e-017-1.527665297025e-017 5.323610690264e-001 -6.733900344896e-0015.910581205868e-002 4.285425323867e-001 1.906901343193e-0011.990794495295e-0016.419444583601e-017 4.121668945107e-017 3.096964582901e-017-5.050360323651e-017 -7.078737686239e-017 -6.0597********e-001 -9.165783032818e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-0015.455993559780e-017 -9.724896332186e-017 3.956603694870e-0171.913857001710e-018 -4.316583456405e-0172.797376665557e-0179.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002-4.501694411183e-002-1.108640877071e-017 4.655237056115e-017 -1.0722********e-017 -9.470236121745e-018 4.277993227986e-017 8.866601870176e-017 -2.516725033065e-016 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001-8.470152033092e-018 9.650816729410e-017 -1.429362211392e-017 -2.789222052040e-017 -3.690793770141e-017 -8.090765462521e-017 -1.964050295697e-016 -6.772274683437e-017 -7.221591336735e-0016.917269588876e-001（3）R2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000-3.558787493836e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-0013.353802899665e-0012.100627753398e-016 -1.961603277854e+000 2.407523727633e-0017.054714572823e-001 5.957204318279e-001 5.526978774676e-001-3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002-3.300197935770e-016 -1.872873410421e-016 2.404534601993e+0001.706758096328e+000 -4.239566704091e-001 3.405332305815e+000-1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-0013.0773********e-017 1.746386351950e-017 0.000000000000e+0001.577122080722e+000 6.399535133956e-001 3.468127872427e-001-5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-0021.760039865880e-016 9.988285339980e-017 0.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000-2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-0021.996629079956e-0018.804885435596e-017 4.996802050991e-017 0.000000000000e+0000.000000000000e+000 8.831633340975e-017 1.231641451542e+0001.619701003419e-001 1.962638275504e-001 5.350035621760e-001-1.509273424767e-001-7.728357669400e-018 -4.385868928757e-018 0.000000000000e+0000.000000000000e+000 -7.751835246841e-018 -1.279231078565e-017-7.753441914209e-001 -3.464514508821e-001 4.312226803504e-0011.234643696237e-001-5.603391361152e-017 -3.179943413314e-017 0.000000000000e+0000.000000000000e+000 -5.620413613517e-017 -9.274974944455e-0170.000000000000e+000 1.296312940612e+000 -4.288053318338e-0012.737334158165e-001-2.493361499851e-017 -1.414991023727e-017 0.000000000000e+0000.000000000000e+000 -2.500935953597e-017 -4.127119443934e-0170.000000000000e+000 0.000000000000e+000 -6.707396440648e-001-4.842320121884e-001-2.603055667460e-017 -1.477242832192e-017 0.000000000000e+0000.000000000000e+000 -2.610963355436e-017 -4.308689959101e-0170.000000000000e+000 0.000000000000e+000 -1.110223024625e-0167.168323926323e-002（4）RQ1.163074414164e+0002.632670934508e+000 -1.772796003272e+000-8.668899138521e-002 3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-0011.836115060851e+000 1.144286420080e-001 -9.880381403133e-0015.589725694767e-001 4.694190067101e-002 -2.978478237007e-0011.617130577649e-002 6.936977702522e-001 1.367670571405e-0011.419099231519e-0025.598200524418e-016 -2.118520153533e+000 -1.876189745783e+000-5.407071940597e-001 1.171538359721e+000 -2.550323020223e+0001.691577936540e+000 1.229951613262e+000 1.387947777212e+0008.667502917242e-001-3.179059303049e-017 -5.261714527977e-017 -8.471995127808e-0014.382910468318e-001 -1.008632199185e+000 -7.959374261495e-0014.769258865577e-001 4.072683083890e-001 4.096390493527e-0013.363378940862e-001-3.514195999269e-016 3.391949386044e-017 -2.160938214545e-016 -1.432244342447e+000 -5.742284908055e-001 1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002-3.704735750656e-017 -1.0998********e-016 -1.953241363386e-0178.269089741494e-017 6.556779598004e-001 -9.275250974463e-0012.529079844053e-001 6.905949216976e-001 -2.374430675823e-002-2.429781119781e-001-6.420051822287e-017 3.865817713597e-017 -3.806718328740e-0172.129734126613e-017 7.853*********e-017 4.698400884876e-001-2.730776009527e-001 7.821296259798e-001 -9.580964936399e-0027.846239841323e-0021.478496020372e-016 -8.383952267131e-017 9.577540382744e-017-7.120338911078e-017 -1.220876056602e-016 -1.854471832043e-0161.287679058937e+000 -3.576058900348e-001 -4.116725408807e-0033.914268216423e-0014.477158378610e-017 -6.204017427568e-017 3.360322998507e-017-1.133882337819e-017 -2.759056827456e-017 -9.217582575819e-0172.214639994444e-016 -3.628760503545e-001 7.398980975354e-0017.241608309576e-0023.408891482791e-017 2.353495494255e-017 1.932283926688e-017-3.456910153081e-017 -2.015177337156e-017 -8.084422691634e-017 -5.839476980893e-017 5.551115123126e-017 -5.176670596524e-0024.958522909877e-002（5）特征值(6.360627875745e-001,0.000000000000e+000)(5.650488993501e-002,0.000000000000e+000)(9.355889078188e-001,0.000000000000e+000)(-9.805309562902e-001,1.139489127430e-001)(-9.805309562902e-001,-1.139489127430e-001)(1.577548557113e+000,0.000000000000e+000)(-1.484039822259e+000,0.000000000000e+000)(-2.323496210212e+000,8.930405177200e-001)(-2.323496210212e+000,-8.930405177200e-001)(3.383039617436e+000,0.000000000000e+000)（6）实特征值的特征向量4.745031936539e+0003.157868541750e+0001.729946912420e+001-1.980049331455e+000-3.187521973524e+0017.794009603201e+000-1.004255685845e+0011.670757770480e+0011.310524273054e+0011.000000000001e+000下一个实特征值对应的特征向量：-5.105003830625e+000-4.886319842360e+0009.505161576151e+000-6.788331788214e-001-9.604334110499e+000-3.0457********e+0001.574873885602e+001-7.395037126277e+000-7.109654943661e+0001.000000000001e+000下一个实特征值对应的特征向量：2.792418944529e+0001.598236841512e+000-5.207507040911e-001-1.667886451702e+000-1.225708535859e+0017.241214790777e+000-5.398214501433e+0002.841008912974e+001-1.216518754416e+0011.000000000001e+000下一个实特征值对应的特征向量：6.217350824581e-001-1.115111815226e-001-2.483773580804e+000-1.306860840421e+000-3.815605442533e+0008.117305509395e+000-1.239170883675e+000-6.800309586197e-0012.691900144840e+0001.000000000001e+000下一个实特征值对应的特征向量：-1.730784582112e+0012.408192534967e+0014.133852365119e-001-8.572044074538e+0009.287334657928e-002-7.832726042776e-002-6.374274025716e-001-3.403204760832e-001。
北航数值分析大作业 第二题 QR分解
《数值分析B》课计算实习第一题设计文档与源程序姓名：杨彦杰学号：SY10171341 算法的设计方案（1）运行平台操作系统：Windows XP；开发平台：VC6.0++；工程类型：文档视图类；工程名：Numanalysis；（2）开发描述首先新建类CMetrix，该类完成矩阵之间的相关运算，包括相乘、加减等，以主程序方便调用；题目的解算过程在视图类CNumanalysisView中实现，解算结果在视图界面中显示；（3）运行流程（4）运行界面2、全部源代码（1）类CMetrixMetrix.h文件：class CMetrix{public:double** MetrixMultiplyConst(double**A,int nRow,int nCol,double nConst);//矩阵乘常数double** MetrixMultiplyMetrix(double**A,double**mA,int nRow,int nCol);//矩阵相乘double** MetrixSubtractMetrix(double **A, double **subA, int nRow,int nCol);//矩阵减矩阵double VectorMultiplyVector(double*V,double*mulV,int nV);//向量点积double** VectorMultiplyVectortoMetrix(double*V,double*VT,int nV);//向量相乘为矩阵double* VectorSubtractVector(double*V,double*subV,int nV);//向量相减double* VectorMultiplyConst(double *V, int nV, double nConst);//向量乘常数double LengthofVector(double *V,int nV);//求向量的长度double* MetrixMultiplyVector(double**A,int nRow,int nCol,double*V,int nV);//矩阵与向量相乘double** AtoAT(double **A,int Row,int Col);//矩阵转置运算void FreeMem();CMetrix(int nRow,int nCol);uCMetrix();virtual ~CMetrix();double* vector; //过渡向量double** B; //过渡矩阵};Metrix.cpp文件：CMetrix::CMetrix(int nRow, int nCol){B = new double*[nRow];for (int i = 0;i < nCol;i++){B[i] = new double[nCol];}vector = new double[nRow];}CMetrix::~CMetrix(){delete vector;B = NULL;delete B;}double** CMetrix::AtoAT(double **A, int nRow, int nCol){for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[col][row] = A[row][col];}}return B;}double* CMetrix::MetrixMultiplyVector(double **A, int nRow, int nCol, double *V, int nV) {if (nCol != nV){AfxMessageBox("矩阵列数和向量维数不等,不能相乘!");return 0;}double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){sum += A[row][col]*V[col];}vector[row] = sum;sum = 0.0;}return vector;}double CMetrix::LengthofVector(double *V, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*V[col];}return length;}double* CMetrix::VectorMultiplyConst(double *V, int nV, double nConst){for (int col = 0;col < nV;col++){vector[col] = V[col]*nConst;}return vector;}double* CMetrix::VectorSubtractVector(double *V, double *subV, int nV){for (int col = 0;col < nV;col++){vector[col] = V[col]-subV[col];}return vector;}double** CMetrix::VectorMultiplyVectortoMetrix(double*V, double *VT, int nV){for (int row = 0;row < nV;row++){for (int col = 0;col < nV;col++){B[row][col] = V[row]*VT[col];}}return B;}double CMetrix::VectorMultiplyVector(double *V, double *mulV, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*mulV[col];}return length;}double** CMetrix::MetrixSubtractMetrix(double **A, double **subA, int nRow, int nCol) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]-subA[row][col];}}return B;}double** CMetrix::MetrixMultiplyMetrix(double **A, double **mA, int nRow, int nCol) {double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){for(int n = 0;n < nCol;n++){sum += A[row][n]*mA[n][col];}B[row][col] = sum;sum = 0.0;}}return B;}double** CMetrix::MetrixMultiplyConst(double **A, int nRow, int nCol, double nConst) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]*nConst;}}return B;}（2）类CNumanalysisViewNumanalysisview.hclass CNumanalysisView : public CEditView{…………public:double Sign(double x);void DisplayVector(double*V,int nV); // 显示向量数据void DisplayMetrix(double **A,int Row,int Col); //显示矩阵void DisplayText(CString str); //显示文本protected://{{AFX_MSG(CNumanalysisView)afx_msg void OnQRanalyze(); //运行主函数…………};Numanalysisview.cppvoid CNumanalysisView::OnQRanalyze(){//开辟空间int nRow = 10;int nCol = 10;CString str;CMetrix Metrix(nRow,nCol);double tempa = 0.0;double *V = new double[nCol]; //分配10*10矩阵空间double *ur = new double[nCol];double *pr = new double[nCol];double *qr = new double[nCol];double *wr = new double[nCol];double *tempV = new double[nCol];double **Ar = new double*[nRow];double **C = new double*[nRow];double **Cr = new double*[nRow];double **tempA = new double*[nRow];double **A = new double*[nRow];double **R = new double*[nRow];for (int col = 0;col < nRow;col++){A[col] = new double[nCol];Ar[col] = new double[nCol];C[col] = new double[nCol];Cr[col] = new double[nCol];tempA[col] = new double[nCol];R[col] = new double[nCol];}//矩阵A求解for (int i = 0;i < nRow;i++){for (int j = 0;j < nCol;j++){if(i == j)A[i][j] = 1.5*cos((i+1.0)+1.2*(j+1.0));elseA[i][j] = sin(0.5*(i+1.0)+0.2*(j+1.0));}}//--------------------拟上三角化-------------------------// double dr = 0.0,cr = 0.0,hr = 0.0,tr = 0.0;for (int r = 0;r < nCol - 2;r++){dr = 0.0;for (i = r+1;i < nCol;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+2;i < nCol;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r+1][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r+1][r])*dr;hr = cr*cr - cr*A[r+1][r]; //hrstr.Format("dr = %.6e, cr = %.6e, hr = %.6e",dr,cr,hr);for (int row = 0;row < nRow;row++) //ur{if (row < r+1)ur[row] = 0.0;else if (row == r+1)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,nRow,nCol,ur,nCol); //pr memcpy(pr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(pr,nCol,1.0/hr);memcpy(pr,tempV,nCol*8);tempV = Metrix.MetrixMultiplyVector(A,nRow,nCol,ur,nCol); //qr memcpy(qr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(qr,nCol,1.0/hr);memcpy(qr,tempV,nCol*8);tr = Metrix.VectorMultiplyVector(pr,ur,nCol)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,nCol,tr); //wr memcpy(wr,tempV,nCol*8);tempV = Metrix.VectorSubtractVector(qr,wr,nCol);memcpy(wr,tempV,nCol*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,nCol); //Arfor (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)A[row][col] = tempA[row][col];}tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}}DisplayText("矩阵A拟上三角化后所得的矩阵为：");DisplayMetrix(A,nRow,nCol);for (int row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)R[row][col] = A[row][col];}// -------------------------------------------------////--------------------带双步位移的QR分解-------------------------// int m = nCol;struct EigenVal //定义特征值结构，实数和虚数{double Realnum;double Imagnum;};EigenVal *eigenvalue = new EigenVal[m];EigenVal tmpEigen1,tmpEigen2;double b = 0.0,c = 0.0,delta = 0.0,s = 0.0,t = 0.0;double *vr = new double[m];for (int k = 1;k < 100; k++){//m代表矩阵阶数，判断式中直接用，运算中需要-1while (m > 1 && fabs(A[m-1][m-2]) <= IPSLEN)//第三步和第四步{eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;m = m - 1;}if (m == 1){eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;DisplayText("已求出A的全部特征值：");break;}b = -(A[m-2][m-2]+A[m-1][m-1]); //第五步求一元二次方程式的根s1,s2c = A[m-2][m-2]*A[m-1][m-1]-A[m-2][m-1]*A[m-1][m-2];delta =b*b - 4*c;if (delta >= 0.0){tmpEigen1.Realnum = (-b-sqrt(delta))/2;tmpEigen1.Imagnum = 0.0;tmpEigen2.Realnum = (-b+sqrt(delta))/2;tmpEigen2.Imagnum = 0.0;}else{tmpEigen1.Realnum = -b/2;tmpEigen1.Imagnum = -sqrt(fabs(delta))/2 ;tmpEigen2.Realnum = -b/2;tmpEigen2.Imagnum = sqrt(fabs(delta))/2;}if (m == 2) //第六步 m=2时结束运算{eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;DisplayText("已求出A的全部特征值：");break;}else //第七步 m > 1{if (fabs(A[m-2][m-3]) <= IPSLEN){eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;m = m - 2;continue;}}for (int row = 0;row < m;row++) //Mk求之前需要把A付给C{for (int col = 0;col < m;col++)C[row][col] = A[row][col];}double **I = new double*[m]; //第九步求Mk和Mk的QR分解for (int i = 0;i < m;i++) //求单位矩阵I,分配m*m矩阵空间{I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}s = A[m-2][m-2]+A[m-1][m-1];t = A[m-2][m-2]*A[m-1][m-1] - A[m-2][m-1]*A[m-2][m-1];tempA = Metrix.MetrixMultiplyMetrix(A,A,m,m);//A*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixMultiplyConst(A,m,m,s);//s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(Ar,A,m,m);//A*A-s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col]; }tempA = Metrix.MetrixMultiplyConst(I,m,m,-1*t);//-t*Ifor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m);//A*A - s*A + r*I for (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}delete I;//Mk的QR分解for (int r = 0;r < m - 1;r++){dr = 0.0;for (i = r;i < m;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+1;i < m;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r][r])*dr;hr = cr*cr - cr*A[r][r]; //hrfor (int row = 0;row < m;row++) //ur{if (row < r)ur[row] = 0.0;else if (row == r)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,m,m); //Btfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,m,m,ur,m); //Bt*ur memcpy(vr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(vr,m,1.0/hr); //vr = Bt*ur/hr memcpy(vr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(ur,vr,m);//Ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m); //Br-ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}tempA = Metrix.AtoAT(C,m,m); //Ctfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempV = Metrix.MetrixMultiplyVector(Cr,m,m,ur,m); //pr memcpy(pr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(pr,m,1.0/hr);memcpy(pr,tempV,m*8);tempV = Metrix.MetrixMultiplyVector(C,m,m,ur,m); //qr memcpy(qr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(qr,m,1.0/hr);memcpy(qr,tempV,m*8);tr = Metrix.VectorMultiplyVector(pr,ur,m)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,m,tr); //wr memcpy(wr,tempV,m*8);tempV = Metrix.VectorSubtractVector(qr,wr,m);memcpy(wr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,m);//Cr+1for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)C[row][col] = tempA[row][col]; }tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++){C[row][col] = tempA[row][col];if (fabs(C[row][col]) < IPSLEN){C[row][col] = 0.0;}}}}str.Format("矩阵A%d QR分解结束后所得到的矩阵为：",m);//计算结果输出DisplayText(str);DisplayMetrix(A,m,m);for (row = 0;row < m;row++) //Mk的QR分解后需要把C付给A{for (col = 0;col < m;col++)A[row][col] = C[row][col];}str.Format("迭代完成后的矩阵A%d = ",k);DisplayText(str);DisplayMetrix(A,m,m);}DisplayText("矩阵A的全体特征值如下: ");for (i = 0;i<nCol;i++){str.Format("%.6e + j%.6e",eigenvalue[i].Realnum,eigenvalue[i].Imagnum);DisplayText(str);}// -------------------------------------------------//求实特征值的特征向量，在拟上三角矩阵基础上直接求解即可////(A-egiI)X = 0.0;m = nRow;for (row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)A[row][col] = R[row][col];}double **I = new double*[m]; //求单位矩阵I,分配m*m矩阵空间double sum = 0.0;for (i = 0;i < m;i++){I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}for (i = 0;i < nRow;i++){if (eigenvalue[i].Imagnum != 0.0){str.Format("特征值%.6e+j%.6e为虚数，不需要求特征向量。
数值分析大作业QR分解
数值分析大作业QR分解题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵20010-1-24A=0-2131431?? ? ? ? ???为例进行具体的求解计算。
一、算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。
其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。
二、算法设计：1、原始实矩阵A Rn n∈拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)rij a i n j r r n ==+。
对于r=1,2,…,n -2执行(1) 若()(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。
(2) 计算令()2()()1,1,1,sgn(0,sgn()=1)r r r r r rr r r r r c aa a ρ+++=-=，（若则取(3) 令-0=r n r u u ?? ?，()()()-1,2,1(,,...,)r r r Tn r r r r r r nr u a c a a ρ++=- (4) 计算(r)(r)(r)Tn-r r+1,r r r+2,r nr r Tn-r T n-r n-r n-r n-r r+1r 1u =(a -c ,a ,...,a )ρI H =I -2uu =H H =I -2u u A =HA H(5) 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。
2、拟上三角矩阵QR 分解的求原矩阵的全部特征值记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)A必然变为上三角矩阵k-1R ，即k-1n-121k-1-1=-1-1-1k-1n-12112n-1R =P P P A =P P P P P PQ （）k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q(2) 计算k k-1k-1A =R Q(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。
QR分解的数值效果报告
基于Householder、givens、gram-schmit、modified gram-schmit方法的QR分解数值效果分析报告班级：0811101学号：081110215姓名：桑树浩目录§1. 摘要 (3)§2. 关键词 (3)§3. 引言 (3)§4. 理论基础 (4)§5. 数值实验 (8)§6. 总结与分析 (19)§7.参考文献与附注 (20)QR分解数值效果分析报告基于数值试验的gram-schmit方法，modified gram-schmit方法，householder变换，givens变换计算矩阵QR分解。
一．摘要：由于工程需要引入了QR分解。
可以利用多次正交变换（householder变换、givens变换、gram-schmit方法，modified gram-schmit方法）来实现QR分解，各个方法各有优劣和其意义，在此就分析一下其各自效果。
经过本次数值实验知道givens变换所用时间较长，但是精确度还是比较高的，householder变换所用时间最短，gram-schmit方法，modified gram-schmit方法所用时间较短，精度最高。
可以根据不同需要选用不同算法。
二．关键词：householder变换givens变换gram-schmit方法modified gram-schmit 方法QR分解数值效果比较三．引言：很多情况下在工程中抽象出来的数学方程组是超定的，没有精确解，这样就需要找一个在某种意义下最接近精确解的解。
设A是m*n的实矩阵，||Ax-b||2=||Q’(Ax-b)||2,这样min||Ay-b||2就等价与min||Q’Ax-Q’b||2,为方便求解，需要Q’A是上三角矩阵，这样引入QR分解就比较求解方便。
可以利用多次正交变换（householder变换、givens变换、gram-schmit方法，modified gram-schmit方法）来实现QR分解1. 豪斯霍尔德变换（Householder transformation）又称初等反射（Elementaryreflection），最初由A.C Aitken在1932年提出[1]。
数值分析1030求矩阵全部特征值的 QR 方法
(2)计算Hk 1 A A j k , k 1, , n n 1 （1）t j ( ul( k 1) alj )
k 1 k 1 ( k 1 ak 1，k ) U ( k 1) (0, ..., 0, k 1 ak 1，k , ak 2,k , ..., ank )
T
xi2 xi21
T
Ri ,i 2 Ri ,i 1 x = x1 , ..., xi-1 , r , 0, 0, xi+3 , ..., xn , r
r
xi2 xi21 xi2 2
T 2 xi2 xn
Ri ,n Ri ,i 2 Ri ,i 1 x = x1 , ..., xi-1 , r , 0, ..., 0 ,
(3)用Ri , j 对矩阵A作变换得到的结论 Ri , j 左乘A只改变A的第i，j行。 RiT, j 右乘A只改变A的第i，j列。 T Ri , j ARi , j 只改变A的第i，j行和第i，j列。
数值分析
例 已知向量x = 2,1, 4 ，试用Givens变换将x约
T
数值分析
分块上三角阵）。即主对角线（或主对角线子块） 及其以下元素均收敛，主对角线（或主对角线子 块）以上元素可以不收敛。特别的，如果A是实对 称阵，则｛Ak ｝ “基本”收敛于对角矩阵。
数值分析
数值分析
QR方法的实际计算步骤
第一步 A 上Hessenberg阵 B
用Householder阵 作正交相似变换
数值分析
数值分析
uuT H2 I 2 T u u 2 2 1 0 4 2 0 1 2 4 0 1 0 0 1 0 H 2 0 0 2 2 2 0 0 2
数值分析矩阵的正交分解(QR分解)
§10 矩阵的正交分解(QR 分解)设nm RA ⨯∈，则存在初等反射阵s H H 1使得)1(2+=s s A A H H (上梯形)[]nmn m m n n a a a a a a a a a a a a A ,,,21212222111211=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=(按列分块) (1)(1)第1步：当01=a 时，取I H =1这一步不需约化，不妨设01≠a ，于是有初等反射阵1H 使1111e a H σ-=，其中Tu u I H 11111--=β。
于是],,,[21)1(1n Ha Ha Ha A H =⎥⎦⎤⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-=222)2(121)2()2(2)2(2)2(22)2(2)2(121000D c B a a a a a a a mn m n n σσ)2(A =其中)2()1(21)2(2)2(222,),,(-⨯--∈∈=n m m T m R D Ra a c (2)第k 步：设已完成对A 上述第1步~第k-1步约化，即存在初等反射阵11,,-k H H 使)(121k k A A H H H =-其中 ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-)()()()(12)2(1)2(1)2(121)(k mn k mkk knk kkk n k k a a a a a a a Aσσσ ⎥⎦⎤⎢⎣⎡=k k k k k D c B r R 0其中)()1(1)()(,],,[k n k m k k n T k m k k kk k R D Ra c c -⨯+-=+-∈∈= ，为 EMBED Equation.3阶上三角阵。
如果0=k c ，这一步不需约化，取I H k =。
不妨设0≠k c ，于是存在初等反射阵kH '使 1e c H k k kσ-='计算kH '的公式： T k k k ku u I H ''-='-1β ⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧+='=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡='=+=∑)(21)()()(22)(,)(,1)(2)()(k kk k k k k k k m k k k k kk k m k i k ik k kk k a u a a a u a a sign σσβσ ………………(2) 令mm k k m k k k R H I H ⨯-+--∈⎥⎦⎤⎢⎣⎡'=111第k 步约化：)1(1)(+==k k k k A A H H A H⎥⎦⎤⎢⎣⎡''⎥⎦⎤⎢⎣⎡'=-k k k k k k kk k D H c H B r R H I 01 )1(121+-=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡'----=k k k k k k k AD H B r σσσσ方框内为第k 步约化需要计算的部分，其中)1(+k A 左上角子阵，1+k R 为k阶上三角阵，这样就使A 三角化过程前进了一步。
北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。
总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。
为了防止矩阵在后面的计算中被破坏保存A[][]。
2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。
由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。
这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。
3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。
用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。
4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。
为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。
若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。
这里用函数intTwoStepDisplacement_QR()来实现。
这是用来存储计算得到的特征值的二维数组。
考虑到特征值可能为复数，因此将所有特征值均当成复数处理。
此函数中，QR分解部分用子函数void QR_decompositionMk()实现。
这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。
5）计算特征向量得到特征值后，计算实特征值相应的特征向量。
这里判断所得特征值的虚数部分是否为零。
求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。
因此先初始化矩阵Array，计算(A-λI)，再求解方程组。
北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。
总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。
为了防止矩阵在后面的计算中被破坏保存A[][]。
2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。
由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。
这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。
3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。
用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。
4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。
为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。
若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。
这里用函数intTwoStepDisplacement_QR()来实现。
这是用来存储计算得到的特征值的二维数组。
考虑到特征值可能为复数，因此将所有特征值均当成复数处理。
此函数中，QR分解部分用子函数void QR_decompositionMk()实现。
这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。
5）计算特征向量得到特征值后，计算实特征值相应的特征向量。
这里判断所得特征值的虚数部分是否为零。
求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。
因此先初始化矩阵Array，计算(A-λI)，再求解方程组。
QR分解 实验报告
并行计算课程考核实验报告考核题目：QR分解算法的并行实现并行实现方式：OpenMP1. 问题概述QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法，一般矩阵先经过正交相似变化成为householder矩阵，然后再应用QR方法求特征值和特征向量。
它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，所以称为QR分解法。
2. 串行代码描述串行主要代码如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#include <time.h>#define n 100 //maxNvoid 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-2vvTdouble 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变换里的乘法和矩阵赋值进行优化。
qr分解例题及答案
qr分解例题及答案QR分解是一种矩阵分解的方法，可以将一个矩阵分解成两个矩阵的乘积，其中一个矩阵为正交矩阵，另一个矩阵为上三角矩阵。
QR分解在数值计算、信号处理、统计学等领域中有广泛的应用。
以下是一个例题及答案，详细说明了如何进行QR分解。
例题：将矩阵A进行QR分解。
A = [1, 2, 3; 4, 5, 6; 7, 8, 7]答案：1. 首先对A的列向量进行正交化。
a1 = [1; 4; 7]v1 = a1 / ||a1|| = [1/√66; 4/√66; 7/√66]a2 = [2; 5; 8]v2 = a2 - (v1' * a2) * v1 = [-0.2673; 0.8729; -0.4099]a3 = [3; 6; 7]v3 = a3 - (v1' * a3) * v1 - (v2' * a3) * v2 = [-0.8018; -0.4364; 0.4082]2. 构造正交矩阵Q，将v1、v2、v3作为Q的列向量。
Q = [v1, v2, v3] =[0.1491, -0.2673, -0.8018;0.5964, 0.8729, -0.4364;0.7887, -0.4099, 0.4082]3. 计算上三角矩阵R，使A = QR。
R = Q' * A =[6.7082, 9.6695, 6.7082;0, 0.8739, 0.8739;0, 0, 1.7321]因此，矩阵A的QR分解为:A = [1, 2, 3; 4, 5, 6; 7, 8, 7] = QR =[0.1491, -0.2673, -0.8018; 0.5964, 0.8729, -0.4364; 0.7887, -0.4099, 0.4082] *[6.7082, 9.6695, 6.7082; 0, 0.8739, 0.8739; 0, 0, 1.7321]QR分解可以应用于线性最小二乘问题、特征值问题、奇异值分解等领域，是求解复杂计算问题的重要方法之一。
数值分析(11)QR方法)
k k ( k a ) (k ) (k ) (k ) (k ) U (0, ..., 0, k akk , ak 1,k , ..., ank )
计算A( k 1) H k A( k )，
( k 1) ( k 1) ( k 1) (k ) (k ) (k) 即 , , , H , H , , H 2 n k 2 k n 1 k 1
1 T Hk Hk H k , Q 1 Q T H n 1 H n 2 1 1 1 2 1 n 1
H 2 H1
数值分析
数值分析
化矩阵A R nn为上三角阵,A 只须依次将各列对角线下元素化为零 (1) (1) (1) 记A A(1) , , , 2 n 1
1
U ( k ) (U ( k ) )T ) (j k )
H
k 1 (k ) (k ) T (k ) (k ) j U (U ) j k 1 (k ) j ( (U ( k ) )T (j k ) )U ( k )， (j k , k 1, , n) k j k , k 1, , n 1 n (k ) (k ) (1)t j ( ul alj )
1、用Householder变换对A作QR分解
有两种情况
数值分析
（1）A R nn非奇异 构造Householder 阵 H k R nn (k 1, 2, , n 1) 则H n1 H n2 H 2 H1 A R （上三角阵）
A H H H R H 1 H 2 H n 1 R QR 其中 Q H 1 H 2 H n 1 R nn为正交阵 R Q 1 A Q T A H n 1 H n 2 H 2 H 1 A
数值分析——带双步位移的QR分解求特征值算法
数 值 分 析（B ） 大 作 业（二）1、算法设计：①矩阵的拟上三角化：对实矩阵A 进行相似变换化为拟上三角矩阵(1)An -，其变换矩阵采用Householder 矩阵，变换过程如下：若()a 0(2,,)r ir i r n ==+，则r H I =；否则，(r)(r)T r r+1,r n,r s =(0,,0,a ,a )， (r)r r+1,r r 2c = -sgn(a )||s ||，()()()r r r r+11,2,u =s -c e (0,,0,,,,)r r r r r r r r nr a c a a ++=-， 2Tr r r 2H =I-2u u /r u ，(1)()r r r r A H A H +=。
当2r n =-时，得(1)(2)(1)222112n n n n n n AH A H H H A H H ------===，令12n-2P=H H H 又r H 是对称正交矩阵，于是n-1T A =P AP 成立，因而n-1A 与 A 相似。
②矩阵的QR 分解：矩阵的QR 分解过程与拟上三角化过程相似，在这里不再重复其原理。
③求全部特征值矩阵拟上三角化后利用带双步位移的QR 方法，采用书本Page 63页具体算法实现。
为了使编程方便，并体会goto 语句使用的灵活性，程序中的跳转均使用goto Loop*实现。
④求A 的相应于实特征值的特征向量求实特征值对应的特征向量，即是求解线性方程组(λI-A)=0i i x ，(1,,)i n =。
因此，为得到全部实特征值对应的特征向量，解线性方程组的过程要循环n 次（n 为矩阵阶数）。
线性方程组的求解采用列主元素Gauss 消去法。
#include <stdio.h>#include <math.h>#define ERR 1.0e-12 //误差限#define N 10 //矩阵行列数#define L 1.0e5 //最大迭代次数double A[N][N]={0};void Init_A() //初始化矩阵{int i,j;for(i=0;i<N;i++)for(j=0;j<N;j++){if(i==j)A[i][j]=1.5*cos((i+1)+1.2*(j+1));elseA[i][j]=sin(0.5*(i+1)+0.2*(j+1));}}/*void Display_A(){int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++)printf("%8.4f",A[i][j]);printf("\n");}}*/int Sgn(double a){if(a>=0)return 1;elsereturn -1;}void On_To_The_Triangle() //矩阵拟上三角化{int i,j,r,flag=0;double cr,dr,hr,tr,temp;double ur[N],pr[N],qr[N],wr[N];for(r=1;r<=N-2;r++){flag=0;for(i=r+2;i<=N;i++)if(A[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r+1;i<=N;i++)dr+=A[i-1][r-1]*A[i-1][r-1];dr=sqrt(dr);if(0==A[r][r-1])cr=dr;else cr=-Sgn(A[r][r-1])*dr;hr=cr*cr-cr*A[r][r-1];for(i=1;i<=r;i++)ur[i-1]=0;ur[r]=A[r][r-1]-cr;for(i=r+2;i<=N;i++)ur[i-1]=A[i-1][r-1];for(i=1;i<=N;i++){temp=0;for(j=1;j<=N;j++)temp+=A[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i<=N;i++){temp=0;for(j=1;j<=N;j++)temp+=A[i-1][j-1]*ur[j-1];qr[i-1]=temp/hr;}temp=0;for(i=1;i<=N;i++){temp+=pr[i-1]*ur[i-1];tr=temp/hr;}for(i=1;i<=N;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i<=N;i++)for(j=1;j<=N;j++)A[i-1][j-1]=A[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}}void Get_Roots(double eigenvalue[][2],int m,double ss,double tt) //求一元二次方程的根{double discriminant=ss*ss-4*tt; //if(discriminant<0){*(*(eigenvalue+m-2))=0.5*ss;*(*(eigenvalue+m-2)+1)=0.5*sqrt(-discriminant);*(*(eigenvalue+m-1))=0.5*ss;*(*(eigenvalue+m-1)+1)=-0.5*sqrt(-discriminant);}else{*(*(eigenvalue+m-2))=0.5*(ss+sqrt(discriminant));*(*(eigenvalue+m-2)+1)=0;*(*(eigenvalue+m-1))=0.5*(ss-sqrt(discriminant));*(*(eigenvalue+m-1)+1)=0;}}void Get_Mk(double mk[][N],int m,double ss,double tt) //获取Mk,用于带双步位移的QR分解{int i,j,k;for(i=0;i<m;i++)for(j=0;j<m;j++)*(*(mk+i)+j)=0;for(i=0;i<m;i++)for(j=0;j<m;j++){for(k=0;k<m;k++)*(*(mk+i)+j)+=A[i][k]*A[k][j];*(*(mk+i)+j)-=ss*A[i][j];if(j==i)*(*(mk+i)+j)+=tt;}}void QR_Reslove(double mk[][N],int m) //QR分解{int i,j,r,flag=0;double cr,dr,hr,tr,temp;double ur[N],vr[N],pr[N],qr[N],wr[N];double B[N][N],C[N][N];for(i=0;i<m;i++)for(j=0;j<m;j++){B[i][j]=*(*(mk+i)+j);C[i][j]=A[i][j];}for(r=1;r<=m-1;r++){flag=0;for(i=r+1;i<=m;i++)if(B[i-1][r-1]!=0){flag=1;break;}if(0==flag)continue;dr=0;for(i=r;i<=m;i++)dr+=B[i-1][r-1]*B[i-1][r-1];dr=sqrt(dr);if(0==B[r-1][r-1])cr=dr;else cr=-Sgn(B[r-1][r-1])*dr;hr=cr*cr-cr*B[r-1][r-1];for(i=1;i<r;i++)ur[i-1]=0;ur[r-1]=B[r-1][r-1]-cr;for(i=r+1;i<=m;i++)ur[i-1]=B[i-1][r-1];for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=B[j-1][i-1]*ur[j-1];vr[i-1]=temp/hr;}for(i=0;i<m;i++)for(j=0;j<m;j++)B[i][j]-=ur[i]*vr[j];for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=C[j-1][i-1]*ur[j-1];pr[i-1]=temp/hr;}for(i=1;i<=m;i++){temp=0;for(j=1;j<=m;j++)temp+=C[i-1][j-1]*ur[j-1];qr[i-1]=temp/hr;}temp=0;for(i=1;i<=m;i++){temp+=pr[i-1]*ur[i-1];tr=temp/hr;}for(i=1;i<=m;i++){wr[i-1]=qr[i-1]-tr*ur[i-1];}for(i=1;i<=m;i++)for(j=1;j<=m;j++)C[i-1][j-1]=C[i-1][j-1]-wr[i-1]*ur[j-1]-ur[i-1]*pr[j-1];}for(i=0;i<m;i++)for(j=0;j<m;j++)A[i][j]=C[i][j];}void Display_Eigenvalue(double value[][2]) //显示特征值{int i;for(i=0;i<N;i++){printf("λ%d=%8.4f",i+1,*(*(value+i)));if(*(*(value+i)+1)>0)printf("+%8.4f",*(*(value+i)+1));else if(*(*(value+i)+1)<0)printf("%8.4f",*(*(value+i)+1));printf("\n");}printf("\n");}int QR_With_Double_Step_Displacement(double eigenvalue[][2]) //带双步位移QR分解求特征值{int i,j,k=1,m=N;double s,t;double Mk[N][N];for(i=0;i<N;i++)for(j=0;j<2;j++)eigenvalue[i][j]=0;do{k++;if(m==1){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}else if(m==2){s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-1][m-2];Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根m=0;continue;}else if(m==0)return 0;else if(fabs(A[m-1][m-2])<=ERR){eigenvalue[m-1][0]=A[m-1][m-1];m--;continue;}else{s=A[m-2][m-2]+A[m-1][m-1];t=A[m-2][m-2]*A[m-1][m-1]-A[m-1][m-2]*A[m-2][m-1];if(fabs(A[m-2][m-3])<=ERR){Get_Roots(eigenvalue,m,s,t); //求一元二次方程的根m-=2;continue;}else{Get_Mk(Mk,m,s,t); //获取Mk,用于带双步位移的QR分解QR_Reslove(Mk,m); //QR分解}}}while(k<L);printf("Errer\n"); //迭代过程不成功，报错}void Select_Principal_Element(int k) //Gauss消元法求解方程组之选主元{int i,max=k;double trans[N];for(i=k+1;i<N;i++)if(fabs(A[i][k])>fabs(A[max][k]))max=i;if(k==max)return;else{for(i=k;i<N;i++){trans[i]=A[k][i];A[k][i]=A[max][i];A[max][i]=trans[i];}}}void Eliminant(int k) //Gauss消元法求解方程组之消元{int i,j;double rate;for(i=k+1;i<N;i++){rate=A[i][k]/A[k][k];for(j=k;j<N;j++)A[i][j]=A[i][j]-rate*A[k][j];}}void Back_Substitution(double Eigenvector[]) //Gauss消元法求解方程组之回代{int i,j,k;double Temp_M[N][N+1];for(i=0;i<N-1;i++){for(j=i;j<N;j++)Temp_M[i][j]=A[i][j];Temp_M[i][N]=0;}Temp_M[N-1][N-1]=1;Temp_M[N-1][N]=1;for(k=N-1;k>=0;k--){*(Eigenvector+k)=Temp_M[k][N];for(i=k;i<N-1;i++)*(Eigenvector+k)=*(Eigenvector+k)-*(Eigenvector+i+1)*Temp_M[k][i+1];*(Eigenvector+k)=*(Eigenvector+k)/Temp_M[k][k];}}void Display_Eigenvector(double Eigenvector[],double eigenvalue) //显示特征值对应特征向量{int i;printf("When λ=%8.4f\nEigenvector=\n",eigenvalue);for(i=0;i<N;i++)printf("%8.4f",*(Eigenvector+i));printf("\n");}void Get_Eigenvector(double value[][2]) //利用Gauss消元法求解特征向量{int i,j;double eigenvalue;double Eigenvector[N];for(i=0;i<N;i++)if(*(*(value+i)+1)==0){Init_A(); //初始化矩阵eigenvalue=*(*(value+i));for(j=0;j<N;j++)A[j][j]-=eigenvalue;for(j=0;j<N-1;j++){Select_Principal_Element(j); //Gauss消元法求解方程组之选主元Eliminant(j); //Gauss消元法求解方程组之消元}Back_Substitution(Eigenvector); //Gauss消元法求解方程组之回代Display_Eigenvector(Eigenvector,eigenvalue); //显示特征值对应特征向量}}main(){double eigenvalue[N][2];Init_A(); //初始化矩阵On_To_The_Triangle(); //矩阵上三角化QR_With_Double_Step_Displacement(eigenvalue); //带双步位移QR分解求特征值printf("Contactme:****************\n\n");Display_Eigenvalue(eigenvalue); //显示特征值Get_Eigenvector(eigenvalue); //利用Gauss消元法求解特征向量return 0;}PS: Forfurtherdetails,pleasecontactme:****************。
QR分解、SVD分解、最小二乘问题数值上机报告
QR分解、SVD分解、最小二乘问题数值上机报告练习6.13 随机构造一个可逆方阵，利用不同的方法给出它的QR分解，观察所得列向量的正交性、CPU时间和所谓的向后稳定性。
解：本题中所得列向量的正交性由||Q T Q−I||2来进行刻画。
向后稳定性由||A−QR||2来进行刻画。
实验步骤：1、首先分别编写出CGS、MGS、Householder、Givens 等不同方法进行矩阵的QR分解的程序。
2、在主程序中调用这些子程序。
其中在编写利用Householder矩阵进行矩阵的QR分解时，可以利用H矩阵作为秩一修正矩阵在计算复杂度方面的优势，即 H∗g=g−b−1(u T g)u,只需2n+1次乘除法运算。
实验数据：（1）取实验矩阵A为满秩的200*200的随机矩阵时得到的数据如下：（2）取实验矩阵A为满秩的400*400的随机矩阵时得到的数据如下：实验数据分析：1.从向后稳定性上来看，CGS和MGS给出的QR乘积均能很好地近似原始矩阵A，因此都是向后稳定的算法。
而Householder方法和Givens方法的向后稳定性略差。
2.从所求得列向量的正交性来看，CGS方法最差，MGS比CGS稍好，Householder方法和Givens方法的表现最好。
Householder方法给出的Q具有更好的列正交性。
3.从CPU时间上来看，Givens方法所耗时间最长，这是合理的，因为对应于一个Householder矩阵，Givens方法中要构造一系列Givens矩阵。
理论分析，对于稠密阵，Householder方法所需的乘除法次数是Givens方法的一半。
且由讲义所指出，有如下结果：N opt House≈∑2(n+1−k)(m+1−k)≈mn2−13n3nk=1N opt GS ≈∑2(k −1)m ≈mn 2n k=1.而由Householder 镜像变换法的计算复杂度略低于MGS 方法，Householder 方法所耗的CPU 时间略低于MGS 方法，故是合理的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目:设计求取n n ⨯实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵20010-1-24A=0-2131431⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭为例进行具体的求解计算。
一、 算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。
其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。
二、 算法设计：1、 原始实矩阵A Rn n⨯∈拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)rij a i n j r r n ==+。
对于r=1,2,…,n -2执行(1) 若()(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。
(2) 计算 令()2()()1,1,1,sgn(0,sgn()=1)r r r r r rr r r r r c aa a ρ+++=-=，（若则取(3) 令-0=r n r u u ⎛⎫ ⎪⎝⎭，()()()-1,2,1(,,...,)r r r Tn r r r r r r nr u a c a a ρ++=- (4) 计算(r)(r)(r)Tn-r r+1,r r r+2,r nr r Tn-r T n-r n-r n-r n-r r+1r 1u =(a -c ,a ,...,a )ρI H =I -2uu =H H =I -2u u A =HA H⎛⎫⎪⎝⎭(5) 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。
2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)A必然变为上三角矩阵k-1R ，即k-1n-121k-1-1=-1-1-1k-1n-12112n-1R =P P P A =P P P P P PQ （）k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q(2) 计算k k-1k-1A =R Q(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。
3、 带位移的QR 方法求实矩阵A Rn n⨯∈全部特征值一般情况下，QR 分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR 算法所得矩阵序列k A 的收敛。
0A =A （A 是原始矩阵的拟上三角矩阵）(1) 分解k-1-1k-1k-1A -u I=R k Q ，其中-1u k 即位移量，一般取A 的某一特征值的近似值；实际计算通常取为k-1A 的右下角元素(k-1)n,n a ，或取为右下角22⨯矩阵特征之中接近(k-1)n,n a 者。
(2) 计算k k-1k-1-1A =+R u I k Q 。
(3) 重复过程（1）（2）直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。
4、 反幂法求实矩阵A Rn n⨯∈的特征向量记通过QR 分解得到A 的某一特征值的近似值为p ，反幂法步骤如下： (1) 三角分解A -pI =LR 。
(2) 对k=1，2，…，做①当k=1时，令Tu =(1 11) 当k 1≠时，解-1Lu=z k②解Ry =u k③求-1y k 绝对值最大的元素k m ④令k k k z =y /m⑤当k k-1m -m 或k k-1z -z 小于规定的误差界时，停止计算，k z 即为所求的特征向量，k p +1/m 即为对应特征值的更为精确的取值。
三、 程序设计：1、 对生成的矩阵A 进行拟上三角化利用hohd 函数对A 进行householder 变换，得到A 的拟上三角矩阵。
2、 对拟上三角化后的矩阵进行QR 分解和带位移的QR 分解，求矩阵的全部特征值在绝对误差界为-41102⨯的条件下，利用qrfenjie 函数对拟上三角矩阵进行QR 分解，利用dp_qrfenjie 函数进行带位移的QR 分解，比较两者的收敛速度。
3、 反幂法求更精确的特征值和特征向量利用vect 函数和（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为-61102⨯。
4、 输出相关结果。
四、 源程序：1、 hohd 函数function[A]=hohd(a) n=length(a); for i=1:n-2b=a(i+1:n,i); if b(1)>=0c=-sqrt(sum(b.^2)); elsec=sqrt(sum(b.^2)); endrho=sqrt(2*c*(c-b(1))); u1=(b-c*eye(n-i,1))/rho; H1=eye(n-i)-2*u1*u1'; H=eye(n-i);H(i+1:n,i+1:n)=H1; a=H*a*(H'); endA=a2、qrfenjie函数function A=qrfenjie(a)n=length(a);A=a;for i=1:100theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;for j=1:n-1theta(j)=atan(A(j+1,j)/A(j,j));c(j)=cos(theta(j));s(j)=sin(theta(j));P=eye(n);P(j,j)=c(j);P(j+1,j)=-s(j);P(j,j+1)=s(j);P(j+1,j+1)=c(j);invP=eye(n);invP(j,j)=c(j);invP(j+1,j)=s(j);invP(j,j+1)=-s(j);invP(j+1,j+1)=c(j);Q=Q*invP;R=P*A;A=R;endA=R*Q;tor=max(abs(diag(A)-diag(a)));if tor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'3、dp_qrfenjie函数function A=dp_qrfenjie(a)n=length(a);A=a;A=a-a(n,n)*eye(n);theta(n-1)=0;c(n-1)=0;s(n-1)=0;Q=1;for j=1:n-1theta(j)=atan(A(j+1,j)/A(j,j));c(j)=cos(theta(j));s(j)=sin(theta(j));P=eye(n);P(j,j)=c(j);P(j+1,j)=-s(j);P(j,j+1)=s(j);P(j+1,j+1)=c(j);invP=eye(n);invP(j,j)=c(j);invP(j+1,j)=s(j);invP(j,j+1)=-s(j);invP(j+1,j+1)=c(j);Q=Q*invP;R=P*A;A=R;endA=R*Q+a(n,n)*eye(n);tor=max(abs(diag(A)-diag(a)));if tor>5*10^-5a=A;elsebreak;endendi,Ak=A,lanmda=diag(A)'4、vect函数function z=vect(a0,p)n=length(a0);Z=zeros(n,n+2);for x=1:na=a0-p(x)*eye(n);r=zeros(n,n);r(1,1:n)=a(1,1:n);l(2:n,1)=a(2:n,1)/a(1,1);for i=2:nfor j=i:nM=0;for k=1:i-1M=l(i,k)*r(k,j)+M;endr(i,j)=a(i,j)-M;endif i<nfor j=i+1:nN=0;for k=1:i-1N=l(j,k)*r(k,i)+N;endl(j,i)=(a(j,i)-N)/r(i,i);endelsebreak;endendR=r;L=l;u=ones(n,1);mo=0;invr=inv(R);invl=inv(L);for i=1:100y=invr*u;p1=max(y);p2=max(-y);if p1>p2m=p1;elsem=-p2;endz=y/m;tor=abs(m-mo);if tor<5*10^-7break;elsemo=m;endu=invl*z;endZ(x,3:n+2)=z';Z(x,1)=i;Z(x,2)=p(x)+1/m;endZ五、运行结果：1、执行命令：>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1];hohd(a);得到原始实数矩阵的拟上三角矩阵A。
A =2.0000 -1.0000 0.0000 0.0000-1.0000 1.0000 5.0000 -0.00000.0000 5.0000 -2.2000 -0.40000.0000 -0.0000 -0.4000 2.20002、执行命令：>> clear,clc,A=[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];qrfenjie(A);得到迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。
i =38Ak =-5.9068 0.0315 -0.0000 0.00000.0315 4.8972 -0.0000 -0.00000.0000 -0.0000 2.2138 0.00070.0000 0.0000 0.0007 1.7958lanmda =-5.9068 4.8972 2.2138 1.79583、执行命令：>> clear,clc,A=[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，以及A的特征值矩阵lamda。
i =8Ak =-5.9068 -0.0056 -0.0000 -0.0000-0.0056 4.8973 0.0000 0.0000-0.0000 0.0000 1.7958 0.00000.0000 0.0000 0.0000 2.2138lanmda =-5.9068 4.8973 1.7958 2.21384、执行命令：>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1],p=[-5.9068 4.8972 2.2138 1.7958];vect(a,p);得到Z矩阵，Z矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。