数值计算程序特征值和特征向量 [转]

//数值计算程序-特征值和特征向量 

////////////////////////////////////////////////////////////// 
//约化对称矩阵为三对角对称矩阵 
//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 
//a-长度为n*n的数组，存放n阶实对称矩阵 
//n-矩阵的阶数 
//q-长度为n*n的数组，返回时存放Householder变换矩阵 
//b-长度为n的数组，返回时存放三对角阵的主对角线元素 
//c-长度为n的数组，返回时前n-1个元素存放次对角线元素 
void eastrq(double a[],int n,double q[],double b[],double c[]); 
////////////////////////////////////////////////////////////// 
//求实对称三对角对称矩阵的全部特征值及特征向量 
//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 
//n-矩阵的阶数 
//b-长度为n的数组，返回时存放三对角阵的主对角线元素 
//c-长度为n的数组，返回时前n-1个元素存放次对角线元素 
//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组 
// 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组 
//a-长度为n*n的数组，存放n阶实对称矩阵 
int ebstq(int n,double b[],double c[],double q[],double eps,int l); 
////////////////////////////////////////////////////////////// 
//约化实矩阵为赫申伯格(Hessen berg)矩阵 
//利用初等相似变换将n阶实矩阵约化为上H矩阵 
//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵 
//n-矩阵的阶数 
void echbg(double a[],int n); 
////////////////////////////////////////////////////////////// 
//求赫申伯格(Hessen berg)矩阵的全部特征值 
//返回值小于0表示超过迭代jt次仍未达到精度要求 
//返回值大于0表示正常返回 
//利用带原点位移的双重步QR方法求上H矩阵的全部特征值 
//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); 
////////////////////////////////////////////////////////////// 
//求实对称矩阵的特征值及特征向量的雅格比法 
//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 
//返回值小于0表示超过迭代jt次仍未达到精度要求 
//返回值大于0表示正常返回 
//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 
//n-矩阵的阶数 
//u-长度为n*n的数组，返回特征向量(按列存储) 
//eps-控制精度要求 
//jt-整型变量，控制最大迭代次数 
int eejcb(double a[],int n,double v[],double eps,int jt); 
////////////////////////////////////////////////////////////// 



选自<<徐世良数值计算程序集(C)>> 

每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。 

今天都给贴出来了 

#include "stdio.h" 
#include "math.h" 
//约化对称矩阵为三对角对称矩阵 
//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 
//a-长度为n*n的数组，存放n阶实对称矩阵 
//n-矩阵的阶数 
//q-长度为n*n的数组，返回时存放Householder变换矩阵 
//b-长度为n的数组，返回时存放三对角阵的主对角线元素 
//c-长度为n的数组，返回时前n-1个元素存放次对角线元素 
void eastrq(double a[],int n,double q[],double b[],double c[]) 
{ 
int i,j,k,u,v; 
double h,f,g,h2; 
for (i=0; i<=n-1; i++) 
{ 
for (j=0; j<=n-1; j++) 
{ 
u=i*n+j; q[u]=a[u]; 
} 
} 
for (i=n-1; i>=1; i--) 
{ 
h=0.0; 
if (i>1) 
{ 
for (k=0; k<=i-1; k++) 
{ 
u=i*n+k; 
h=h+q[u]*q[u]; 
} 
} 

if (h+1.0==1.0) 
{ 
c[i-1]=0.0; 
if (i==1) 
{ 
c[i-1]=q[i*n+i-1]; 
} 
b[i]=0.0; 
} 
else 
{ 
c[i-1]=sqrt(h); 
u=i*n+i-1; 
if (q[u]>0.0) 
{ 
c[i-1]=-c[i-1]; 
} 
h=h-q[u]*c[i-1]; 
q[u]=q[u]-c[i-1]; 
f=0.0; 
for (j=0; j<=i-1; j++) 
{ 
q[j*n+i]=q[i*n+j]/h; 
g=0.0; 
for (k=0; k<=j; k++) 
{ 
g=g+q[j*n+k]*q[i*n+k]; 
} 
if (j+1<=i-1) 
{ 
for (k=j+1; k<=i-1; k++) 
{ 
g=g+q[k*n+j]*q[i*n+k]; 
} 
} 
c[j-1]=g/h; 
f=f+g*q[j*n+i]; 
} 
h2=f/(h+h); 
for (j=0; j<=i-1; j++) 
{ 
f=q[i*n+j]; 
g=c[j-1]-h2*f; 
c[j-1]=g; 
for (k=0; k<=j; k++) 
{ 
u=j*n+k; 
q[u]=q[u]-f*c[k-1]-g*q[i*n+k]; 
} 
} 
b[i]=h; 
} 
} 
b[0]=0.0; 
for (i=0; i<=n-1; i++) 
{ 
if ((b[i]!=0.0)&&(i-1>=0)) 
{ 
for (j=0; j<=i-1; j++) 
{ 
g=0.0; 
for (k=0; k<=i-1; k++) 
{ 
g=g+q[i*n+k]*q[k*n+j]; 
} 
for (k=0; k<=i-1; k++) 
{ 
u=k*n+j; 
q[u]=q[u]-g*q[k*n+i]; 
} 
} 
} 
u=i*n+i; 
b[i]=q[u]; 
q[u]=1.0; 
if (i-1>=0) 
{ 
for (j=0; j<=i-1; j++) 
{ 
q[i*n+j]=0.0; 
q[j*n+i]=0.0; 
} 
} 
} 
return; 




//求实对称三对角对称矩阵的全部特征值及特征向量 
//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量 
//n-矩阵的阶数 
//b-长度为n的数组，返回时存放三对角阵的主对角线元素 
//c-长度为n的数组，返回时前n-1个元素存放次对角线元素 
//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组 
// 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组 
//a-长度为n*n的数组，存放n阶实对称矩阵 
int ebstq(int n,double b[],double c[],double q[],double eps,int l) 
{ 
int i,j,k,m,it,u,v; 
double d,f,h,g,p,r,e,s; 
c[n-1]=0.0; 
d=0.0; 
f=0.0; 
for (j=0; j<=n-1; j++) 
{ 
it=0; 
h=eps*(fabs(b[j])+fabs(c[j])); 
if (h>d) 
{ 
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 ((l>0)&&(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 (w>0.0) 
{ 
xy=1.0; 
if (b<0.0) 
{ 
xy=-1.0; 
} 
u[m-1]=(-b-xy*y)/2.0; 
u[m-2]=c/u[m-1]; 
v[m-1]=0.0; v[m-2]=0.0; 
} 
else 
{ 
u[m-1]=-b/2.0; 
u[m-2]=u[m-1]; 
v[m-1]=y/2.0; 
v[m-2]=-v[m-1]; 
} 
m=m-2; 
it=0; 
} 
else 
{ 
if (it>=jt) 
{ 
printf("fail\n"); 
return(-1); 
} 
it=it+1; 
for (j=l+2; j<=m-1; j++) 
{ 
a[j*n+j-2]=0.0; 
} 
for (j=l+3; j<=m-1; j++) 
{ 
a[j*n+j-3]=0.0; 
} 
for (k=l; k<=m-2; k++) 
{ 
if (k!=l) 
{ 
p=a[k*n+k-1]; 
q=a[(k+1)*n+k-1]; 
r=0.0; 
if (k!=m-2) 
{ 
r=a[(k+2)*n+k-1]; 
} 
} 
else 
{ 
x=a[ii]+a[ll]; 
y=a[ll]*a[ii]-a[kk]*a[jj]; 
ii=l*n+l; 
jj=l*n+l+1; 
kk=(l+1)*n+l; 
ll=(l+1)*n+l+1; 
p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y; 
q=a[kk]*(a[ii]+a[ll]-x); 
r=a[kk]*a[(l+2)*n+l+1]; 
} 
if ((fabs(p)+fabs(q)+fabs(r))!=0.0) 
{ 
xy=1.0; 
if (p<0.0) 
{ 
xy=-1.0; 
} 
s=xy*sqrt(p*p+q*q+r*r); 
if (k!=l) 
{ 
a[k*n+k-1]=-s; 
} 
e=-q/s; 
f=-r/s; 
x=-p/s; 
y=-x-f*r/(p+s); 
g=e*r/(p+s); 
z=-x-e*q/(p+s); 
for (j=k; j<=m-1; j++) 
{ 
ii=k*n+j; 
jj=(k+1)*n+j; 
p=x*a[ii]+e*a[jj]; 
q=e*a[ii]+y*a[jj]; 
r=f*a[ii]+g*a[jj]; 
if (k!=m-2) 
{ 
kk=(k+2)*n+j; 
p=p+f*a[kk]; 
q=q+g*a[kk]; 
r=r+z*a[kk]; 
a[kk]=r; 
} 
a[jj]=q; 
a[ii]=p; 
} 
j=k+3; 
if (j>=m-1) 
{ 
j=m-1; 
} 
for (i=l; i<=j; i++) 
{ 
ii=i*n+k; 
jj=i*n+k+1; 
p=x*a[ii]+e*a[jj]; 
q=e*a[ii]+y*a[jj]; 
r=f*a[ii]+g*a[jj]; 
if (k!=m-2) 
{ 
kk=i*n+k+2; 
p=p+f*a[kk]; 
q=q+g*a[kk]; 
r=r+z*a[kk]; 
a[kk]=r; 
} 
a[jj]=q; 
a[ii]=p; 
} 
} 
} 
} 
} 
return(1); 
} 



//求实对称矩阵的特征值及特征向量的雅格比法 
//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量 
//返回值小于0表示超过迭代jt次仍未达到精度要求 
//返回值大于0表示正常返回 
//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值 
//n-矩阵的阶数 
//u-长度为n*n的数组，返回特征向量(按列存储) 
//eps-控制精度要求 
//jt-整型变量，控制最大迭代次数 
int eejcb(double a[],int n,double v[],double eps,int jt) 
{ 
int i,j,p,q,u,w,t,s,l; 
double fm,cn,sn,omega,x,y,d; 
l=1; 
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; 
} 
} 
} 
while (1==1) 
{ 
fm=0.0; 
for (i=0; i<=n-1; i++) 
{ 
for (j=0; j<=n-1; j++) 
{ 
d=fabs(a[i*n+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; 
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; 
} 
} 
return(1); 
}