24小时热门版块排行榜

返回列表

sleefd

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 26.5
帖子: 27
在线: 7.4小时
虫号: 1331238
注册: 2011-06-26
专业: 计算机软件

[求助] 求讲解一下，下面jacobi算法代码的具体实现思路

int jcb(double a[][200],double tzz[200],double v[][200],double eps,int jt,int rows)

//jacobi求实对称矩阵的特征值和特征向量
//a 为原矩阵，v为特征向量,tzz为特征值,eps为精度,jt为最大迭代次数
{
int i,j,p,q,l;
double fm,cn,sn,omega,x,y,d;
// double temp;
l=1;
for(i=0;i {
v=1.0;
for(j=0;j if(i!=j) v[j]=0;
}
while(true)
{
fm=0.0;
for(i=0;i<=rows-1;i++)
for(j=0;j<=i-1;j++)
{
d=fabs(a[j]);
if((i!=j)&&(d>fm)) //取最大的a[p][q]
{
fm=d;p=i;q=j;}
}
      if(fm {

for(i=0;i {
         tzz=a;
         //printf("%lf \n",a);
}
return(1);
break;
}
if(l>jt) //如果不能求解
{
for(i=0;i {
         tzz=a;
         //printf("%lf \n",a);
}
return(-1);
      break;
}
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]=y*omega+a[p][q]*(2*cn*cn-1);
a[q][p]=a[p][q];
for(j=0;j<=rows-1;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<=rows-1;i++)
if((i!=p)&&(i!=q))
{
fm=a[p];
a[p]=fm*cn+a[q]*sn;
a[q]=-fm*sn+a[q]*cn;
}
for(i=0;i<=rows-1;i++)
{
fm=v[p];
v[p]=fm*cn+v[q]*sn;
v[q]=-fm*sn+v[q]*cn;

}
}
}