//實(shí)對(duì)稱(chēng)陣特征值及特征向量雅可比法
#include <iostream.h>
#include <math.h>

int  main()
{
	cout.precision(2);
	cout.setf(ios::fixed);
	float a[3][3]={
		{2.0,-1.0,0.0},
		{-1.0,2.0,-1.0},
		{0.0,-1.0,2.0}
	};
	float v[3][3];
	float eps=0.001;//精度
    int	jt=10;//迭代次數(shù)



    int j, p, q, l(1), stRank;
    float fm,cn,sn,omega,x,y,d;
	
	//if(MatrixSymmetry(a)!=true)	//不是對(duì)稱(chēng)陣
	//	return(0);					

//	stRank = a.GetColNum();		// 矩陣階數(shù)
      stRank=3;
    for(int i=0; i<stRank; i++)
    {
		v[i][i]=1.0;
        for(j=0; j<stRank; j++)
          if(i!=j) v[i][j]=0.0;
    }
    while(1)
    { 
		fm=0.0;
        for(i=1; i<stRank; i++)
        for(j=0; j<i; j++)
        {
			d=fabs(a[i][j]);
            if((i!=j)&&(d>fm))
            {
				fm=d; 
				p=i;
				q=j;
			}
        }
        if(fm<eps)  
		{
			cout<<"特征向量為："<<endl;
			for(int st=0; st<stRank; st++)
			{
        for(int tt=0; tt<stRank; tt++)
		{cout<<v[st][tt]<<"       ";
         //cout<<a[st][tt]<<"   ";
		}
		cout<<endl;
		
			}
			cout<<endl;
			cout<<"對(duì)角線為三個(gè)特征值："<<endl;
			for(int s1t=0; s1t<stRank; s1t++)
			{
        for(int t1t=0; t1t<stRank; t1t++)
		{
         cout<<a[s1t][t1t]<<"   ";
		}
		cout<<endl;
			}
			return 1;
		}
        if(l>jt)  return(-1);
        l=l+1;
        x=-a[p][q];
		y=(a[q][q]-a[p][p])/2.0;
        omega=x/sqrt(x*x+y*y);
        if(y<0.0) omega=-omega;
        sn=1.0+sqrt(1.0-omega*omega);
        sn=omega/sqrt(2.0*sn);
        cn=float (sqrt(1.0-sn*sn));
        fm=a[p][p];
        a[p][p]=fm*cn*cn+a[q][q]*sn*sn+a[p][q]*omega;
        a[q][q]=fm*sn*sn+a[q][q]*cn*cn-a[p][q]*omega;
        a[p][q]=0.0; 
		a[q][p]=0.0;
        for(j=0; j<stRank; j++)
        if((j!=p)&&(j!=q))
        {
            fm=a[p][j];
            a[p][j]=fm*cn+a[q][j]*sn;
            a[q][j]=-fm*sn+a[q][j]*cn;
        }
        for(i=0; i<stRank; i++)
          if((i!=p)&&(i!=q))
          { 
              fm=a[i][p];
              a[i][p]=fm*cn+a[i][q]*sn;
              a[i][q]=-fm*sn+a[i][q]*cn;
          }
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        for(i=0; i<stRank; i++)
        {
            fm=v[i][p];
            v[i][p]=fm*cn+v[i][q]*sn;
            v[i][q]=-fm*sn+v[i][q]*cn;
        }
    }
	
}