c语言求矩阵特征向量函数 c语言求矩阵的迹( 三 )


if (hd)
{
d=h;
}
m=j;
while ((m=n-1)(fabs(c[m])d))
{
m=m+1;
}
if (m!=j)
{
do
{
if (it==l)
{
printf("fail\n");
return(-1);
}
it=it+1;
g=b[j];
p=(b[j+1]-g)/(2.0*c[j]);
r=sqrt(p*p+1.0);
if (p=0.0)
{
b[j]=c[j]/(p+r);
}
else
{
b[j]=c[j]/(p-r);
}
h=g-b[j];
for (i=j+1; i=n-1; i++)
{
b[i]=b[i]-h;
}
f=f+h;
p=b[m];
e=1.0;
s=0.0;
for (i=m-1; i=j; i--)
{
g=e*c[i];
h=e*p;
if (fabs(p)=fabs(c[i]))
{
e=c[i]/p;
r=sqrt(e*e+1.0);
c[i+1]=s*p*r;
s=e/r;
e=1.0/r;
}
else
{
e=p/c[i];
r=sqrt(e*e+1.0);
c[i+1]=s*c[i]*r;
s=1.0/r;
e=e/r;
}
p=e*b[i]-s*g;
b[i+1]=h+s*(e*g+s*b[i]);
for (k=0; k=n-1; k++)
{
u=k*n+i+1;
v=u-1;
h=q[u];
q[u]=s*q[v]+e*h;
q[v]=e*q[v]-s*h;
}
}
c[j]=s*p;
b[j]=e*p;
}
while (fabs(c[j])d);
}
b[j]=b[j]+f;
}
for (i=0; i=n-1; i++)
{
k=i; p=b[i];
if (i+1=n-1)
{
j=i+1;
while ((j=n-1)(b[j]=p))
{
k=j;
p=b[j];
j=j+1;
}
}
if (k!=i)
{
b[k]=b[i];
b[i]=p;
for (j=0; j=n-1; j++)
{
u=j*n+i;
v=j*n+k;
p=q[u];
q[u]=q[v];
q[v]=p;
}
}
}
return(1);
}
//约化实矩阵为赫申伯格(Hessen berg)矩阵
//利用初等相似变换将n阶实矩阵约化为上H矩阵
//a-长度为n*n的数组,存放n阶实矩阵,返回时存放上H矩阵
//n-矩阵的阶数
void echbg(double a[],int n)
{ int i,j,k,u,v;
double d,t;
for (k=1; k=n-2; k++)
{
d=0.0;
for (j=k; j=n-1; j++)
{
u=j*n+k-1;
t=a[u];
if (fabs(t)fabs(d))
{
d=t;
i=j;
}
}
if (fabs(d)+1.0!=1.0)
{
if (i!=k)
{
for (j=k-1; j=n-1; j++)
{
u=i*n+j;
v=k*n+j;
t=a[u];
a[u]=a[v];
a[v]=t;
}
for (j=0; j=n-1; j++)
{
u=j*n+i;
v=j*n+k;
t=a[u];
a[u]=a[v];
a[v]=t;
}
}
for (i=k+1; i=n-1; i++)
{
u=i*n+k-1;
t=a[u]/d;
a[u]=0.0;
for (j=k; j=n-1; j++)
{
v=i*n+j;
a[v]=a[v]-t*a[k*n+j];
}
for (j=0; j=n-1; j++)
{
v=j*n+k;
a[v]=a[v]+t*a[j*n+i];
}
}
}
}
return;
}
//求赫申伯格(Hessen berg)矩阵的全部特征值
//利用带原点位移的双重步QR方法求上H矩阵的全部特征值
//返回值小于0表示超过迭代jt次仍未达到精度要求
//返回值大于0表示正常返回
//a-长度为n*n的数组,存放上H矩阵
//n-矩阵的阶数
//u-长度为n的数组,返回n个特征值的实部
//v-长度为n的数组,返回n个特征值的虚部
//eps-控制精度要求
//jt-整型变量,控制最大迭代次数
int edqr(double a[],int n,double u[],double v[],double eps,int jt)
{
int m,it,i,j,k,l,ii,jj,kk,ll;
double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
it=0;
m=n;
while (m!=0)
{
l=m-1;
while ((l0)(fabs(a[l*n+l-1])eps*(fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l]))))
{
l=l-1;
}
ii=(m-1)*n+m-1;
jj=(m-1)*n+m-2;
kk=(m-2)*n+m-1;
ll=(m-2)*n+m-2;
if (l==m-1)
{
u[m-1]=a[(m-1)*n+m-1];
v[m-1]=0.0;
m=m-1; it=0;
}
else if (l==m-2)
{
b=-(a[ii]+a[ll]);
c=a[ii]*a[ll]-a[jj]*a[kk];
w=b*b-4.0*c;
y=sqrt(fabs(w));
if (w0.0)
{
xy=1.0;

推荐阅读