guoguan.c

jacobi过关法求解矩阵特征值主函数

#include"math.h"
void jcbm(int n,double a[],double v[],double eps)
{ int i,j,p,q,u,w,t,s;
double ff,fm,cn,sn,omega,x,y,d;
for(i=0;i<=n-1;i++)
{ v[i*n+i]=1.0;
for(j=0;j<=n-1;j++)
if(i!=j) v[i*n+j]=0.0;}
ff=0.0;

for(i=1;i<=n-1;i++)
for(j=0;j<=i-1;j++)
{
	d=a[i*n+j];
    ff=ff+d*d;}
    ff=sqrt(2.0*ff);
loop0:
    ff=ff/(1.0*n);
loop1:
for(i=1;i<=n-1;i++)
for(j=0;j<i-1;j++)
{
  d=fabs(a[i*n+j]);
if(d>ff)
{  p=i;q=j;
goto loop;
}
}
if(ff<eps) return;
goto loop0;
loop:
u=p*n+q;w=p*n+p;t=q*n+p;s=q*n+q;
x=-a[u];y=(a[s]-a[w])/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[w];
a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;
a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;
a[u]=0.0;a[t]=0.0;
for(j=0;j<=n-1;j++)
if((j!=p)&&(j!=q))
{u=p*n+j;
w=q*n+j;
fm=a[u];
a[u]=fm*cn+a[w]*sn;
a[w]=-fm*sn+a[w]*cn;
}
for(i=0;i<=n-1;i++)
if((i!=p)&&(i!=q))
{u=i*n+p;
 w=i*n+q;
fm=a[u];
a[u]=fm*cn+a[w]*sn;
a[w]=-fm*sn+a[w]*cn;
}
for(i=0;i<=n-1;i++)
{u=i*n+p;
 w=i*n+q;
fm=v[u];
v[u]=fm*cn+v[w]*sn;
v[w]=-fm*sn+v[w]*cn;}
goto loop1;
}