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

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

47.{

48.for (i=0; 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

110.dblEigenValue = GetElement(i,i);

111.

112.return TRUE;

113.}