北航数值分析大作业 第二题 QR分解文档格式.docx
《北航数值分析大作业 第二题 QR分解文档格式.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业 第二题 QR分解文档格式.docx（35页珍藏版）》请在冰豆网上搜索。
判断B[i][r]对角线下是否为零
doubleVectMod（vector<
p）
求向量的模
ConveArraMultiVect（vector<
B,vector<
矩阵的转置乘向量
ArraMultiVect（vector<
矩阵乘向量
NumbMultiVect（vector<
vect,doublea）
数乘向量
VectMultiCovVect（vector<
a,vector<
b）
向量乘向量转置得到矩阵
voidArraSubs（vector<
C,
D）
两个矩阵相减
VectSubs（vector<
两个向量相减
voidGausElim（vector<
a）
列主元高斯消元法求齐次方程解向量
voidStop（vector<
Ar）
停止，结束程序
voidSolutS1S2（complex<
s1,complex<
s2,
求解s1，s2；
voidSave2（complex<
s2）
保存两个特征值
voidSave1（complex<
s）
保存一个特征值
voidJudgemBelow2（vector<
A,
Abk）
对于m==1及m==0的处理
voidHessenberg（vector<
矩阵拟上三角化
voidQRMethod（vector<
A）
矩阵QR分解
voidCalculatAk（vector<
Ak）
带双步位移QR迭代法
**********************************************************/
3.//QRmethod.cpp:
4.//
5.
6.#include"
stdafx.h"
7.#include<
vector>
8.#include<
iostream>
9.#include<
math.h>
10.#include<
complex>
11.#include<
fstream>
12.usingnamespacestd;
13.
14.
15.constdoubleepsion=1e-12;
16.constintL=20000;
17.intm,n;
18.intk=1;
19.
20.vector<
I;
21.vector<
complex<
Lambda;
22.
23.
24.
25.
26.vector<
27.{
28.inta=A.size（）;
29.doubles=A[a-2][a-2]+A[a-1][a-1];
30.doublet=A[a-2][a-2]*A[a-1][a-1]-A[a-1][a-2]*A[a-2][a-1];
31.vector<
D（a,vector<
（a,0））;
32.for（inti=1;
i<
a;
i++）
33.{
34.for（intj=1;
j<
j++）
35.{
36.doublesum=0;
37.for（intk=1;
k<
k++）
38.{
39.sum+=A[i][k]*A[k][j];
40.}
41.D[i][j]=sum-s*A[i][j]+t*（i==j?
1.0:
0）;
42.}
43.}
44.returnD;
45.}
46.
47.doubleVectMultiVect（vector<
48.{
49.vector<
:
size_typea=y.size（）;
50.doublevalue=0;
51.for（vector<
size_typei=1;
52.{
53.value+=y[i]*u[i];
54.}
55.returnvalue;
56.}
57.
58.doubleGetMode（vector<
59.{
60.doublevalue=0;
61.inta=B.size（）;
62.for（intk=r;
63.{
64.value+=B[k][r]*B[k][r];
65.}
66.value=sqrt（value）;
67.returnvalue;
68.}
69.
70.doubleGetModeH（vector<
71.{
72.doublevalue=0;
73.inta=B.size（）;
74.for（intk=r+1;
75.{
76.value+=B[k][r]*B[k][r];
77.}
78.value=sqrt（value）;
79.returnvalue;
80.}
81.
82.vector<
83.{
84.intb=D.size（）;
85.vector<
U（b,vector<
（b,0））;
86.for（inti=1;
b;
87.{
88.for（intj=1;
89.{
90.U[i][j]=a*D[i][j];
91.}
92.}
93.returnU;
94.}
95.
96.boolIsBirZero（vector<
97.{
98.boolb=true;
99.inta=B.size（）;
100.for（inti=r+1;
101.{
102.if（abs（B[i][r]）>
epsion）
103.{
104.b=false;
105.}
106.}
107.returnb;
108.}
109.
110.boolIsBirZeroH（vector<
111.{
112.boolb=true;
113.inta=B.size（）;
114.for（inti=r+2;
115.{
116.if（abs（B[i][r]）>
117.{
118.b=false;
119.}
120.}
121.returnb;
122.}
123.
124.doubleVectMod（vector<
125.{
126.doublevalue=0.0;
127.vector<
size_typei,j;
128.j=p.size（）;
129.for（i=1;
j;
130.{
131.value+=p[i]*p[i];
132.}
133.value=sqrt（value）;
134.returnvalue;
135.}
136.
137.vector<
138.{
139.inta=B.size（）;
140.intb=u.size（）;
141.vector<
vec（a,0）;
142.if（a!
=b）
143.{
144.cerr<
<
"
ArrayandVectornotmatchinsize!
"
;
145.}
146.else
147.{
148.for（inti=1;
149.{
150.for（intj=1;
151.{
152.vec[i]+=B[j][i]*u[j];
153.}
154.}
155.returnvec;
156.}
157.}
158.
159.vector<
B,
160.vector<
161.{
162.inta=B.size（）;
163.intb=u.size（）;
164.vector<
165.if（a!
166.{
167.cerr<
168.}
169.else
170.{
171.for（inti=1;
172.{
173.for（intj=1;
174.{
175.vec[i]+=B[i][j]*u[j];
176.}
177.}
178.returnvec;
179.}
180.}
181.
182.vector<
183.{
184.intj=vect.size（）;
185.vector<
b（j,0）;
186.for（inti=1;
187.{
188.b[i]=a*vect[i];
189.}
190.returnb;
191.}
192.
193.vector<
a,
194.vector<
195.{
196.ints1=a.size（）;
197.ints2=b.size（）;
198.if（s1!
=s2）
199.{
200.cerr<
Vectorsnotmatchinsize!
201.}
202.else
203.{
204.vector<
U（s1,vector<
（s1,0））;
205.for（inti=1;
s1;
206.{
207.for（intj=1;
208.{
209.U[i][j]=a[i]*b[j];
210.}
211.}
212.returnU;
213.}
214.}
215.
216.voidArraSubs（vector<
217.vector<
218.{
219.inta=C.size（）;
220.intb=D.size（）;
221.intc=C[0].size（）;
222.intd=D[0].size（）;
223.if（a!
=b||c!
=d）
224.{
225.cerr<
226.}
227.else
228.{
229.for（inti=1;
230.{
231.for（intj=1;
232.{
233.C[i][j]-=D[i][j];
234.}
235.}
236.}
237.}
238.
239.
240.vector<
241.vector<
242.{
243.ints1=a.size（）;
244.ints2=b.size（）;
245.if（s1!
246.{
247.cerr<
248.}
249.else
250.{
251.vector<
value（s1,0）;
252.for（inti=1;
i<
i++）
253.{
254.value[i]=a[i]-b[i];
255.}
256.returnvalue;
257.}
258.}
259.
260.voidGausElim（vector<
261.{
262.intsz=a.size（）;
263.vector<
int>
fx,ufx,ufxp,record;
264.vector<
x（sz,0）;
265.vector<
ret;
266.//*****消元过程********
267.for（intk=1;
k<
sz-1;
268.{
269.doublemax=a[k][k];
270.intp=0;
271.for（inti=k+1;
sz;
272.{
273.if（abs（a[i][k]）>
abs（max））
274.{
275.p=i;
276.max=a[i][k];
277.}
278.}//选出主元行
279.if（abs（max）>
=epsion）
280.{
281.if（p!
=0）
282.{
283.for（intj=k;
284.{
285.doubletemp=a[k][j];
286.a[k][j]=a[p][j];
287.a[p][j]=temp;
288.}//交换主元行
289.}
290.
291.for（inti=k+1;
292.{
293.doublett=a[i][k]/a[k][k];
294.for（intj=k+1;
295.{
296.a[i][j]=a[i][j]-tt*a[k][j];
297.}
298.}
299.
300.}//endofif（abs（max）>
301.}//endoffor（intk=1;
sz;
302.
303.for（inti=1;
304.{
305.intp=0;
306.for（intj=1;
=i;
j++）
307.{
308.
309.if（abs（a[j][i]）>
=10*epsion）//不为零
310.{
311.p=j;
312.}
313.}
314.record.push_back（p）;
315.}
316.
317.for（inti=1;
318.{
319.intp=0;
320.for（intj=sz-2;
j>
=0;
j--）
321.{
322.if（record[j]==i）
323.{
324.p=j+1;
325.}
326.}
327.if（p!
328.{
329.ufxp.push_back（p）;
330.ufx.push_back（i）;
331.}
332.}
333.
334.for（inti=1;
335.{
336.intp=0;
337.for（intj=0;
j<
ufxp.size（）;
338.{
339.if（ufxp[j]==i）
340.{
341.p=1;
342.}
343.}
344.if（p==0）
345.{
346.fx.push_back（i）;
347.}
348.}//endoffor（inti=1;
349.intc=fx.size（）;
350.//*****************************************************
351.for（inti=0;
c;
352.{
353.for（intj=0;
354.{
355.if（i==j）
356.{
357.x[fx.at（j）]=1.0+epsion;
358.}
359.else
360.{
361.x[fx.at（j）]=0;
362.}
363.}
364.
365.intb=ufxp.size（）;
366.for
升级会员 