雅克比法求取矩阵特征值和特征向量（c语言）

#include "stdio.h" #include "math.h" int jcbi(double a[][5],int n,double v[][5],double eps,int jt) { int i,j,p,q,u,w,l; double fm,cn,sn,omega,x,y,d; l=1; for (i=0; i<=n-1; i++) { for (j=0; j<=n-1; j++) if (i!=j) v[i][j]=0.0; else v[i][j] = 1.0; } while (1==1) { fm=0.0; for (i=1; i<=n-1; i++) for (j=0; j<=i-1; j++) { d=fabs(a[i][j]); if ((i!=j)&&(d>fm)) { fm=d; p=i; q=j;} } if (fm<eps) 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=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<=n-1; j++) if ((j!=p)&&(j!=q)) { u=p*n+j; w=q*n+j; fm=a[p][q]; a[p][q]=fm*cn+a[p][p]*sn; a[p][p]=-fm*sn+a[p][p]*cn; } for (i=0; i<=n-1; i++) if ((i!=p)&&(i!=q)) { u=i*n+p; w=i*n+q; fm=a[p][q]; a[p][q]=fm*cn+a[p][p]*sn; a[p][p]=-fm*sn+a[p][p]*cn; } for (i=0; i<=n-1; i++) { fm=v[i][p]; v[i][p]=fm*cn+v[i][q]*sn; v[i][q]=-fm*sn+v[i][q]*cn; } } return(1); } void main() { int i,j,k; double eps,b[5],temp; static double v[5][5]; static double a[5][5]={ {10.0,1.0,2.0,3.0,4.0}, {1.0,9.0,-1.0,2.0,-3.0}, {2.0,-1.0,7.0,3.0,-5.0}, {3.0,2.0,3.0,12.0,-1.0}, {4.0,-3.0,-5.0,-1.0,15.0}}; eps=0.000001; jcbi(a,5,v,eps,100); for (i=0; i<=4; i++) { printf("%13.7e\n",a[i][i]); b[i]=a[i][i]; } printf("\n\n"); for (i=0; i<=4; i++) { for (j=0; j<=4; j++) printf("%12.6e ",v[i][j]); printf("\n"); } printf("\n"); for(i=0;i<5;i++) for(j=0;j<5-i;j++) if(b[j]<b[j+1]) {temp=b[j]; b[j]=b[j+1]; b[j+1]=temp; } for(i=0;i<5;i++) printf("%12.6e\n ",b[i]); scanf("%d",&k); }