北航数值分析大作业 第二题 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