关于矩阵的特征值和特征向量求解，大部分的数学运算库都进行了提供，下面是使用MTL库的接口进行封装。

#include 
    
     #include 
     
      #include 
      
       #include 
       
        #include 
        
         /*! 对角阵 */typedef mtl::matrix 
         
          , mtl::dense<>, mtl::row_major>::type DiagMatrix;/*! 对称阵 */typedef mtl::matrix 
          
           , mtl::packed
           
            , mtl::column_major>::type SymmMatrix;/*! 标准矩阵，数值存在矩阵中 */typedef mtl::matrix 
            
             , mtl::dense<>, mtl::row_major>::type Matrix;/*! 标准矩阵，数值存在外部 */typedef mtl::matrix 
             
              , mtl::dense
              
               , mtl::row_major>::type ExtMatrix;/*! 标准向量，数值存在矩阵中 */typedef mtl::dense1D
               
                 Vector;/*! 标准向量，数值存在外部 */typedef mtl::external_vec
                
                  ExtVector;/*** @brief 求实对称矩阵的特征值及特征向量* @param MMatrix 所求的矩阵* @param EigenValues 矩阵的特征值* @param EigenVectors 矩阵的特征向量（按照列来存），如果特征值的序号为1，那么对应的特征向量为第1列* @param Percent 矩阵的特征值百分比* @param AccPercent 矩阵的特征值累积百分比* @param deps 累计误差* @return 返回值小于0表示超过迭代jt次仍未达到精度要求，返回值大于0表示正常返回 */ bool GetMatrixEigen(Matrix MMatrix, Vector &EigenValues, Matrix &EigenVectors, Vector *Percent = NULL, Vector *AccPercent = NULL, double deps = 0.00000001);
                
               
              
             
            
           
          
         
        
       
      
     
    

下面是求解矩阵特征值和特征向量的实现：

/*** @brief 求实对称矩阵的特征值及特征向量的雅格比法 * 利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 * @param a		长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 * @param n		矩阵的阶数 * @param v		长度为n*n的数组，返回特征向量(按列存储) * @param eps	控制精度要求 * @param jt		整型变量，控制最大迭代次数 * @return 返回false表示超过迭代jt次仍未达到精度要求，返回true表示正常返回 */ bool Eejcb(double a[], int n, double v[], double eps, int jt){ 	int i,j,p,q,u,w,t,s,l; 	double fm,cn,sn,omega,x,y,d; 	l=1; 	//初始化特征向量矩阵使其全为0 	for(i=0; i<=n-1; i++) 	{   		v[i*n+i] = 1.0; 		for(j=0; j<=n-1; j++) 		{ 			if(i != j)   				v[i*n+j]=0.0; 		} 	} 	while(true)   //循环 	{   		fm = 0.0; 		for(i=0; i<=n-1; i++)   //找出,矩阵a(特征值),中除对角线外其他元素的最大绝对值 		{ 			//这个最大值是位于a[p][q] ,等于fm 			for(j=0; j<=n-1; j++) 			{ 				d = fabs(a[i*n+j]); 				if((i!=j) && (d> fm)) 				{ 					fm = d;   					p = i;   					q = j; 				} 			} 		} 		if(fm < eps)     //精度复合要求 			return true; //正常返回 		if(l > jt)       //迭代次数太多 			return false;//失败返回 		l ++;       //   迭代计数器 		u = p*n + q; 		w = p*n + p;   		t = q*n + p;   		s = q*n + q; 		x = -a[u]; 		y = (a[s]-a[w])/2.0;		//x y的求法不同 		omega = x/sqrt(x*x+y*y);	//sin2θ 		//tan2θ=x/y = -2.0*a[u]/(a[s]-a[w]) 		if(y < 0.0) 			omega=-omega; 		sn = 1.0 + sqrt(1.0-omega*omega);   		sn = omega /sqrt(2.0*sn);		//sinθ 		cn = sqrt(1.0-sn*sn);			//cosθ 		fm = a[w];   //   变换前的a[w]   a[p][p] 		a[w] = fm*cn*cn + a[s]*sn*sn + a[u]*omega; 		a[s] = fm*sn*sn + a[s]*cn*cn - a[u]*omega; 		a[u] = 0.0; 		a[t] = 0.0; 		//   以下是旋转矩阵,旋转了了p行,q行,p列,q列 		//   但是四个特殊点没有旋转(这四个点在上述语句中发生了变化) 		//   其他不在这些行和列的点也没变 		//   旋转矩阵,旋转p行和q行 		for(j=0; j<=n-1; j++) 		{ 			if((j!=p) && (j!=q)) 			{ 				u = p*n + j; 				w = q*n + j; 				fm = a[u]; 				a[u] = a[w]*sn + fm*cn; 				a[w] = a[w]*cn - fm*sn; 			} 		} 		//旋转矩阵,旋转p列和q列 		for(i=0; i<=n-1; i++) 		{ 			if((i!=p) && (i!=q)) 			{ 				u = i*n + p;   				w = i*n + q; 				fm = a[u]; 				a[u]= a[w]*sn + fm*cn; 				a[w]= a[w]*cn - fm*sn; 			} 		} 		//记录旋转矩阵特征向量 		for(i=0; i<=n-1; i++) 		{ 			u = i*n + p;   			w = i*n + q; 			fm = v[u]; 			v[u] =v[w]*sn + fm*cn; 			v[w] =v[w]*cn - fm*sn; 		} 	} 	return true; }bool GetMatrixEigen(Matrix MMatrix, Vector &EigenValues, Matrix &EigenVectors, 					Vector *Percent, Vector *AccPercent, double deps){	int nDem = MMatrix.nrows();	double *mat = new double[nDem*nDem];	//输入矩阵，计算完成之后保存特征值	double *eiv = new double[nDem*nDem];	//计算完成之后保存特征向量	for (int i=0; i
    
      mapEigen;	for (size_t index=0; index
     
      ::reverse_iterator iter = mapEigen.rbegin();	for(size_t ii=0; ii
      
       first;		int index = iter->second;	//获取该特征值对应的特征向量对应的列号		if(Percent != NULL)	//计算百分比以及累积百分比		{			double dTemp = iter->first / dSumEigenValue;			(*Percent)[ii] = dTemp;			if(AccPercent != NULL)			{				if(ii != 0)					(*AccPercent)[ii] = (*AccPercent)[ii-1] + dTemp;				else					(*AccPercent)[ii] = dTemp;			}		}		double dMaxValue = ABS(TmpEigenVectors(0, index));		int iNum = 0;		for(int iRow=0; iRow
      
     
    

求解最核心的方法是雅可比方法。关于雅可比方法网上有很多资料，这里不再进行说明。

参考资料：

[1]

[2]

[3]

[4]