Jacobi法求实对称矩阵的特征值与特征向量(c语言,不要复制的,能注释清楚更好)

(1) 编制的程序能求解任意解实对称矩阵;
(2) 把正交变换的操作编制成独立的功能函数;
(3) 输出每一步的p、q、正交矩阵、正交变换的最终结果、全部特征值及对应的特征向量。

第1个回答  2013-11-05
void Jacobi(int n,float (*a)[7][7],float (*s)[8][8]) //n为矩阵阶数,a为输入矩阵,s为输出矩阵{ int i,j,i1,l,iq,iql,ip; float g,s1,s2,s3,v1,v2,v3,u,st,ct; for(i=0;i<=n;i++) { for(j=0;j<=i;j++) { if((i-j)==0) (*s)[i][j]=1; else { (*s)[i][j]=0.0; (*s)[j][i]=0.0; } } } g=0; for(i=1;i<n;i++) { i1=i-1; for(j=0;j<=i1;j++) g=g+2*(*a)[i][j]*(*a)[i][j]; } s1=sqrt(g); s2=min/(n)*s1; s3=s1; do { s3=s3/(n); do { l=0; for(iq=1;iq<n;iq++) { iql=iq-1; for(ip=0;ip<=iql;ip++) { if(fabs((*a)[ip][iq])>=s3) { l=1; v1=(*a)[ip][ip]; v2=(*a)[ip][iq]; v3=(*a)[iq][iq]; u=0.5*(v1-v3); if(u==0) g=1; if(fabs(u)>=1e-10) g=-(u/fabs(u)*l)*v2/sqrt(v2*v2+u*u); st=g/sqrt(2*(l+sqrt(l-g*g))); ct=sqrt(1-st*st); for(i=0;i<n;i++) { g=(*a)[i][ip]*ct-(*a)[i][iq]*st; (*a)[i][iq]=(*a)[i][ip]*st+(*a)[i][iq]*ct; (*a)[i][ip]=g; g=(*s)[i][ip]*ct-(*s)[i][iq]*st; (*s)[i][iq]=(*s)[i][ip]*st+(*s)[i][iq]*ct; (*s)[i][ip]=g; } for(i=0;i<n;i++) { (*a)[ip][i]=(*a)[i][ip]; (*a)[iq][i]=(*a)[i][iq]; } g=2*v2*st*ct; (*a)[ip][ip]=v1*ct*ct+v3*st*st-g; (*a)[iq][iq]=v1*st*st+v3*ct*ct+g; (*a)[ip][iq]=(v1-v3)*st*ct+v2*(ct*ct-st*st); (*a)[iq][ip]=(*a)[ip][iq]; } } } }while(l==1); }while(s3>s2); }
相似回答