求解对称矩阵全部特征值和特征向量算法(雅可比)

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;
}