求解對稱矩陣全部特征值和特征向量算法(雅可比)

ValVece.h

#include <iostream>

class ValVect {
public :
	ValVect(void);
	void clear(void);
	//~ValVect(void);

public :
	void rdOrMatrix(int _dim, double _e); //從文件中讀取對稱方陣,將其存在矩valAry中(每一列為一樣本向量,即行數為維數)._dim為方陣的階,_e為計算精度
	void itVtAry(void); //初始化樣本向量矩陣vectAry(每一列為一特征向量)
	void calValVect(void); //計算特征值和特征向量,特征值為valAry個對角線元素,對應的特征向量為vectAry各列
	void stValVect(void); //按特征值由大到小將特征值和相應的特征向量排序
	void showValVect(void); //輸出特征值和相應的特征向量
	void sortVect(void); //輸出排序后的特征向量矩陣

private :
	int dim_; //方陣的階數
	double **valAry, **vectAry; //初始方陣和特征向量方陣
	double e_;
};

 ValVectFile.h

#include "ValVect.h"
#include <math.h>
#include <malloc.h>
#include <iostream>

ValVect::ValVect() {
	dim_ = 0;  e_ = 0;
	valAry = NULL;  vectAry = NULL;
}

void ValVect::rdOrMatrix(int _dim, double _e) { //讀取待計算的方陣
	dim_ = _dim;  e_ = _e;
	valAry = (double **)malloc(dim_*sizeof(double));
	for(int i=0; i<dim_; i++)
		*(valAry+i) = (double *)malloc(dim_*sizeof(double)); //valAry為初始方陣,階為dim_,每一列為一個樣本向量,列數為樣本個數,以后在矩陣本身計算特征值

	std::cout << "read the original vectors" << '\n';
	FILE *fp;
	char ch;
	if((fp=fopen("ValVectOrMatrix.txt","rt"))==NULL) { //ValVectOrMatrix文件中每一行為一個樣本向量,列數代表維數                             
		printf("Cannot open file strike any key exit!");
		exit(0);
	}

	ch=fgetc(fp);
	std::string str=""; //這里str是類string的一個對象
	int k=0;
	while(ch!=EOF) {
		for(int i=0; i<dim_; i++) {
			while(ch==' ' || ch=='v' || ch=='\n')
				ch=fgetc(fp);
			while(ch!=' ' && ch!='\n' && ch!='v') {
				str+=ch;
				ch=fgetc(fp);
			}
			double aa=atof(str.c_str()); //將字符串轉換成double型
			*(*(valAry+i)+k) = aa; //文本中的一行為矩陣valAry中的一列
			str="";
		}
		ch=fgetc(fp);
		k++;
	}
}

void ValVect::itVtAry() { //初始化特征向量矩陣
	vectAry = (double **)malloc(dim_*sizeof(double));
	for(int i=0; i<dim_; i++)
		*(vectAry+i) = (double *)malloc(dim_*sizeof(double)); //特征向量矩陣為vectAry,每一列為響應特征值對應的特征向量;矩陣的初始狀態為單位矩陣
	for(i=0; i<dim_; i++)
		for(int j=0; j<dim_; j++) {
			if(i==j)
				*(*(vectAry+i)+j) = 1;
			else *(*(vectAry+i)+j) = 0;
		}
}

void ValVect::calValVect() { //計算方陣的特征值和特征向量矩陣
	//int p, q; //滿足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)
	double b, c, cn, sn; //b=(*(*(valAry+p)+p)- *(*(valAry+q)+q))/2, c=sqrt(b^2+(*(*(valAry+p)+q))^2)
								//cn=(1/2+|b|/(2c))*(1/2), sn=+(*(*(valAry+p)+q)/(2c*cn))(b<=0) sn=-(*(*(valAry+p)+q)/(2c*cn))(b>0)
								//但是當*(*(valAry+p)+p)=*(*(valAry+q)+q)是cn=sn=sqrt(2)/2
	int count =1; //記錄轉換次數
	double sum = 0.0; //計算非對角線元素平方和,與e_比較判斷結束計算條件
	for(int i=0; i<dim_; i++)
		for(int j=0; j<i; j++)
			if(i!=j)
				sum += *(*(valAry+i)+j)* (*(*(valAry+i)+j));
	while(2*sum > e_) {
		double max = 0.0;
		int p=0, q=0; //滿足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)
		//找到滿足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)的p,q
		for(i=0; i<dim_; i++)
			for(int j=0; j<dim_; j++)
				if(i!=j)
					if(fabs(*(*(valAry+i)+j)) > fabs(max)) {
						max = *(*(valAry+i)+j);
						p = i;	q = j;
					}
		std::cout << "p=" << p << ";q=" << q <<'\n';
		//計算用於建立旋轉矩陣Q的 cn,sn值
		std::cout << *(*(valAry+p)+p) << ',' << *(*(valAry+q)+q)  <<'\n'; 
		if( (double)(*(*(valAry+p)+p))== (double)(*(*(valAry+q)+q)) || (int)(*(*(valAry+p)+p))== (int)(*(*(valAry+q)+q))) {
			cn = sqrt(2)/2;
			sn = sqrt(2)/2;
		}
		else {
			b = ( *(*(valAry+p)+p)- (*(*(valAry+q)+q)) )/2;
			c = sqrt(b*b- *(*(valAry+p)+q)* (*(*(valAry+p)+q)));
			cn = 1.0/2+ fabs(b)/(2*c);
			if(b<=0)
				sn = *(*(valAry+p)+q)/(2*c*cn);
			else sn = -1* (*(*(valAry+p)+q)/(2*c*cn));
		}
		std::cout << "cn=" << cn << ";sn=" << sn << '\n';

		//計算旋轉變化后的valAry方陣
		double nwPP = *(*(valAry+p)+p)* (cn*cn) + *(*(valAry+q)+q)* (sn*sn) - *(*(valAry+p)+q)*sn*cn*2;
		double nwQQ = *(*(valAry+p)+p)* (sn*sn) + *(*(valAry+q)+q)* (cn*cn) + *(*(valAry+p)+q)*sn*cn*2;
		double nwPQ = ( *(*(valAry+p)+p) - *(*(valAry+q)+q) )*sn*cn + *(*(valAry+p)+q)*(cn*cn-sn*sn);
		if(fabs(nwPP) < e_) nwPP = 0.0;
		if(fabs(nwQQ) < e_) nwQQ = 0.0;
		if(fabs(nwPQ) < e_) nwPQ = 0.0;

		double **temp; //存儲兩行臨時值p,q行
		temp = (double **)malloc(2*sizeof(double));
		for(i=0; i<2; i++)
			*(temp+i) = (double *)malloc(dim_*sizeof(double));
		for(i=0; i<dim_; i++)
			if(i!=p && i!=q) { 
				*(*(temp+0)+i) = *(*(valAry+p)+i)*cn - *(*(valAry+q)+i)*sn; //計算新p行元素
				//if(fabs(*(*(temp+0)+i)) < e_) *(*(temp+0)+i) = 0.0;
				*(*(temp+1)+i) = *(*(valAry+p)+i)*sn + *(*(valAry+q)+i)*cn; //計算新q行元素
				//if(fabs(*(*(temp+1)+i)) < e_) *(*(temp+1)+i) = 0.0;
			}
		*(*(valAry+p)+p) = nwPP; //更新PP行	
		*(*(valAry+q)+q) = nwQQ; //更新QQ行
		*(*(valAry+p)+q) = nwPQ; //更新PQ行	
		*(*(valAry+q)+p) = nwPQ; //更新QP行
		for(i=0; i<dim_; i++)
			if(i!=p && i!=q) {
				*(*(valAry+p)+i) = *(*(temp+0)+i); //更新p行
				*(*(valAry+i)+p) = *(*(temp+0)+i); //更行p列
				*(*(valAry+q)+i) = *(*(temp+1)+i); //更新q行
				*(*(valAry+i)+q) = *(*(temp+1)+i); //更行q列
			}

		//檢驗更新是否正確
		std::cout << "轉換" << count << "次之后的特征值矩陣" << '\n';
		for(i=0; i<dim_; i++) {
			for(int j=0; j<dim_; j++)
				std::cout << *(*(valAry+i)+j) << "  ";
			std::cout << '\n';
		}

		//計算旋轉變換之后的特征向量方陣Q
		double **tpVect; //存儲兩列臨時值p,q列
		tpVect = (double **)malloc(dim_*sizeof(double));
		for(i=0; i<dim_; i++)
			*(tpVect+i) = (double *)malloc(2*sizeof(double));
		for(i=0; i<dim_; i++) {
			*(*(tpVect+i)+0) = *(*(vectAry+i)+p)*cn - *(*(vectAry+i)+q)*sn; //存放新p列
			//if(fabs(*(*(tpVect+i)+0)) < e_) *(*(tpVect+i)+0) = 0.0;
			*(*(tpVect+i)+1) = *(*(vectAry+i)+p)*sn + *(*(vectAry+i)+q)*cn; //存放新q列
			//if(fabs(*(*(tpVect+i)+1)) < e_) *(*(tpVect+i)+1) = 0.0;
		}
		for(i=0; i<dim_; i++) {
			*(*(vectAry+i)+p) = *(*(tpVect+i)+0); //更新p列
			*(*(vectAry+i)+q) = *(*(tpVect+i)+1); //更新p列
		}

		//檢驗一下更行后的特征向量矩陣
		std::cout << "轉換" << count << "次之后的特征向量矩陣" << '\n';
		for(i=0; i<dim_; i++) {
			for(int j=0; j<dim_; j++)
				std::cout << *(*(vectAry+i)+j) << "  ";
			std::cout << '\n';
		}
		count++;

		//進行下一輪運算
		sum = 0.0;
		for(i=0; i<dim_; i++) //計算非對角線元素平方和,與e_比較判斷結束計算條件
			for(int j=0; j<i; j++)
				if(i!=j)
					sum += *(*(valAry+i)+j)* (*(*(valAry+i)+j));

	}//end while
}

void ValVect::showValVect() { //輸出特征值和其對應的特征向量
	std::cout << "output valAndVect :" << '\n';
	for(int i=0; i<dim_; i++) {
		std::cout << *(*(valAry+i)+i) << "  "; //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		for(int j=0; j<dim_; j++)
			std::cout << *(*(vectAry+j)+i) << "  ";
		std::cout << '\n';
		}
}

void ValVect::stValVect() { //按特征值由大到小排序並將相應的熱證向量矩陣排序
	double *sort; //特征向量排序時存儲數據的臨時數組
	sort = (double *)malloc(dim_*sizeof(double));
	for(int i=0; i<dim_-1; i++) {
		int flag=i;
		double max = *(*(valAry+i)+i);
		for(int j=i+1; j<dim_; j++)
			if(*(*(valAry+j)+j) > max) {
				max = *(*(valAry+j)+j);
				flag=j;
			}
		if(flag!=i) {
			for(int k=0; k<dim_; k++) { //將相應的特征向量交換
				*(sort+k)= *(*(vectAry+k)+i);
				*(*(vectAry+k)+i)= *(*(vectAry+k)+flag);
				*(*(vectAry+k)+flag)= *(sort+k);
			}
			//將相應的特征值進行交換,否則下一次比較會出錯
			double tmp = *(*(valAry+i)+i);
			*(*(valAry+i)+i)= *(*(valAry+flag)+flag);
			*(*(valAry+flag)+flag)= tmp;
		}
	}
}

void ValVect::sortVect() { //輸出排序后的特征向量矩陣
	std::cout << "排序后的特征向量矩陣(每一列為一個特征向量) :" << '\n';
	for(int i=0; i<dim_; i++) {
		for(int j=0; j<dim_; j++)
			std::cout << *(*(vectAry+i)+j) << "  ";
		std::cout << '\n';
		}
}

 ValVect.cpp

#include <iostream>
#include "ValVectFile.h"
#define er 1.0E-10
int main() {
	ValVect valVect;
	valVect.rdOrMatrix(3,er);
	valVect.itVtAry();
	valVect.calValVect();
	valVect.showValVect();
	valVect.stValVect();
	valVect.sortVect();

	return 0;
}