雅可比算法求方阵的全部特征值和特征向量

ValVect.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)/2int 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,qfor(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';}//计算旋转变换之后的特征向量方阵Qdouble **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';}}

main

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