雅克比过关法求特征值和特征向量

1.//////////////////////////////////////////////////////////////////////
2.// 求实对称矩阵特征值与特征向量的雅可比法
3.//
4.// 参数：
5.// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值
6.// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与
7.// 数组dblEigenValue中第j个特征值对应的特征向量
8.// 3. int nMaxIt - 迭代次数，默认值为60
9.// 4. double eps - 计算精度，默认值为0.000001
10.//
11.// 返回值：BOOL型，求解是否成功
12.//////////////////////////////////////////////////////////////////////
13.BOOL CMatrix::JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
14.{
15.int i,j,p,q,u,w,t,s,l;
16.double fm,cn,sn,omega,x,y,d;
17.
18.if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
19.return FALSE;
20.
21.l=1;
22.for (i=0; i<=m_nNumColumns-1; i++)
23.{
24.mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
25.for (j=0; j<=m_nNumColumns-1; j++)
26.if (i!=j)
27.mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;//单位矩阵
28.}
29.
30.while (TRUE)
31.{
32.fm=0.0;
33.for (i=1; i<=m_nNumColumns-1; i++)
34.{
35.for (j=0; j<=i-1; j++)
36.{
37.d=fabs(m_pData[i*m_nNumColumns+j]);
38.if ((i!=j)&&(d>fm))
39.{
40.fm=d;
41.p=i;
42.q=j;
}//取绝对值最大的非对角线元素，并记住位置
43.}
44.}
45.
46.if (fm<eps)
47.{
48.for (i=0; i<m_nNumColumns; ++i)
49.dblEigenValue = GetElement(i,i);
50.return TRUE;
51.}
52.
53.if (l>nMaxIt)
54.return FALSE;
55.
56.l=l+1;
57.u=p*m_nNumColumns+q;
58.w=p*m_nNumColumns+p;
59.t=q*m_nNumColumns+p;
60.s=q*m_nNumColumns+q;
61.x=-m_pData;
62.y=(m_pData[s]-m_pData[w])/2.0;
63.omega=x/sqrt(x*x+y*y);
64.
65.if (y<0.0)
66.omega=-omega;
67.
68.sn=1.0+sqrt(1.0-omega*omega);
69.sn=omega/sqrt(2.0*sn);//sin
=sqrt(1.0-sn*sn);//cos
71.fm=m_pData[w];
72.m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData*omega;
73.m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData*omega;
74.m_pData=0.0;
75.m_pData[t]=0.0;
76.for (j=0; j<=m_nNumColumns-1; j++)
77.{
78.if ((j!=p)&&(j!=q))
79.{
80.u=p*m_nNumColumns+j; w=q*m_nNumColumns+j;
81.fm=m_pData;
82.m_pData=fm*cn+m_pData[w]*sn;
83.m_pData[w]=-fm*sn+m_pData[w]*cn;
84.}
85.}
86.for (i=0; i<=m_nNumColumns-1; i++)
87.{
88.if ((i!=p)&&(i!=q))
89.{
90.u=i*m_nNumColumns+p;
91.w=i*m_nNumColumns+q;
92.fm=m_pData;
93.m_pData=fm*cn+m_pData[w]*sn;
94.m_pData[w]=-fm*sn+m_pData[w]*cn;
95.}
96.}
97.
98.for (i=0; i<=m_nNumColumns-1; i++)
99.{
100.u=i*m_nNumColumns+p;
101.w=i*m_nNumColumns+q;
102.fm=mtxEigenVector.m_pData;
103.mtxEigenVector.m_pData=fm*cn+mtxEigenVector.m_pData[w]*sn; 104.mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn; 105.}
106.
107.}
108.
109.for (i=0; i<m_nNumColumns; ++i)
110.dblEigenValue = GetElement(i,i);
111.
112.return TRUE;
113.}。