实对称矩阵特征值和特征向量求解
// face.cpp : Defines the entry point for the console application.
//
#include
#include
#include
#include
#include
#include
using namespace std;
int main(int argc, char* argv[])
{
vector
vector
//读取矩阵里的数据给矩阵A
FILE *fp;
char ch;
if((fp=fopen("C:\\matrix.txt","rt"))==NULL)
{
printf("Cannot open file strike any key exit!");
getch();
exit(1);
}
//写下读取时矩阵的行和列到TEXT2文本中，验证读取时是否正确
FILE *fp3;
if((fp3=fopen("C:\\TEXT2.txt","w"))==NULL)
{
printf("无法建立！");
exit(0);
}
ch=fgetc(fp);
string str="";
int num=0;
int allnum=0;
int ia=0,ia2=0,temp=0;
int tempnum=0;
while(ch!=EOF)
{
if(ch==' '||ch=='\n'){
//printf("是空格或换行");
}else
{
str+=ch;
num++;
if(num==6)
{
allnum++;
//printf(" ");
//printf(str.c_str());
double aa=atof(str.c_str());
fprintf(fp3, "%d%s%d%s",ia," ",ia2," # ");
if(ia2<117)
{
A[ia][ia2]=aa;
if(allnum<153){
//cout<
ia=(allnum/13)%120;
//if(ia<120)
//{
ia2++;
ia2=ia2%13+temp*13;
if(ia==119)
{
tempnum++;
cout<
if(tempnum==13)
{
temp++;
tempnum=0;
}
//}
//else
//{
//ia2++;
//temp=ia2/13;
//cout<
//}
}
else if(ia2>116)
{
cout<
ia=((allnum-14041)/3)%120;
//cout<
ia2++;
ia2=ia2%3+temp*13;
}
//cout<
num=0;
}
// putchar(ch);
}
ch=fgetc(fp);
}
cout<
//for(int iac=0;iac<120;iac++){
//for(int iac2=0;iac2<120;iac2++)
//{
//cout<
//cout<
//getch();
//------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------
//下面实现求特征值和特征向量的作用
for(int iv=0;iv<120;iv++){
for(int iv2=0;iv2<120;iv2++)
{
if(iv==iv2){
V[iv][iv2]=1;
}
else
{
V[iv][iv2]=0;
}
}
}
double *eigsv=new double[120];
for(int ieigsv=0;ieigsv<120;ieigsv++)
{
eigsv[ieigsv]=0;
}
double epsl=0.0001;
int maxt=10;
int n=120;
int success=0; // 函数返回值
double tao, t, cn, sn; // 临时变量
double maxa;
//写特征向量到C:\\tezhengxiangliang.txt
FILE *fp2;
if((fp2=fopen("C:\\tezhengxiangliang.txt","w"))==NU
LL)
{
printf("无法建立！");
exit(0);
}
//写特征值到C:\\tezhengzhi.txt
FILE *fp4;
if((fp4=fopen("C:\\tezhengzhi.txt","w"))==NULL)
{
printf("无法建立！");
exit(0);
}
for (int it=0; it
maxa=0;
for (int p=0; p
for (int q=p+1; q
if (fabs(A[p][q])>maxa) // 记录非对角线元素最大值
{
maxa=fabs(A[p][q]);
}
//cout<<"fabs(A[p][q]"<
{
// 计算Givens旋转矩阵的重要元素:cos(theta), 即cn, sin(theta), 即sn
tao=0.5*(A[q][q]-A[p][p])/A[p][q];
if (tao>=0) // t为方程 t^2 + 2*t*tao - 1 = 0的根, 取绝对值较小的根为t
{
t=-tao+sqrt(1+tao*tao);
}
else
{
t=-tao-sqrt(1+tao*tao);
}
cn=1/sqrt(1+t*t);
sn=t*cn;
// Givens旋转矩阵之转置左乘A, 即更新A的p行和q行
for (int j=0; j
double apj=A[p][j];
double aqj=A[q][j];
A[p][j]=cn*apj-sn*aqj;
A[q][j]=sn*apj+cn*aqj;
}
// Givens旋转矩阵右乘A, 即更新A的p列和q列
for (int i=0; i
double aip=A[i][p];
double aiq=A[i][q];
A[i][p]=cn*aip-sn*aiq;
A[i][q]=sn*aip+cn*aiq;
}
// 更新特征向量存储矩阵V, V=J0×J1×J2...×Jit, 所以每次只更新V的p, q两列
for (int i2=0; i2
double vip=V[i2][p];
double viq=V[i2][q];
V[i2][p]=cn*vip-sn*viq;
V[i2][q]=sn*vip+cn*viq;
}
}
}
}
if (maxa
// 特征值向量
for (int j2=0; j2
eigsv[j2]=A[j2][j2];
//fprintf(fp2, "%f ",eigsv[j2]);
//fprintf(fp2, "%s ","##");
}
// 对特征值向量从大到小进行排序, 并调整特征向量顺序 (直接插入法)
double* tmp=new double[n];
for (int j=1; j
int i=j;
double a=eigsv[j];
for (int k=0; k
tmp[k]=V[k][j];
}
while(i>0 && eigsv[i-1]
for (int k=0; k
V[k][i]=V[k][i-1];
}
i--;
}
eigsv[i]=a;
for (int k2=0; k2
V[k2][i]=tmp[k2];
}
}
delete[] tmp;
cout<<"Hello World!"<<"非对角线元素已小于收敛标准，迭代结束"<
for(int ivc2=0;ivc2<120;ivc2++)
{
//cout<
}
fprintf(fp2, "%s","#\n");
//cout<
fclose(fp2);
for(int ieigsvc=0;ieigsvc<120;ieigsvc++){
fprintf(fp4, "%f%s",eigsv[ieigsvc],"\n");
//cout<
fclose(f
p2);
return success=1;
}
}
//---------------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------------
//jacbobi_loop(A,V,eigsv,0.0,100,10);
// printf("Hello World!\n"+c);
printf("Hello World!\n");
return success;
}