//////////////////////////////////////////////////////////////
//利用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个元素存放次对角线元素
//数值计算程序-特征值和特征向量//////////////////////////////////////////////////////////////
//约化对称矩阵为三对角对称矩阵
//利用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=0i<=n-1i++)
{
for (j=0j<=n-1j++)
{
u=i*n+jq[u]=a[u]
}
}
for (i=n-1i>=1i--)
{
h=0.0
if (i>1)
{
for (k=0k<=i-1k++)
{
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=0j<=i-1j++)
{
q[j*n+i]=q[i*n+j]/h
g=0.0
for (k=0k<=jk++)
{
g=g+q[j*n+k]*q[i*n+k]
}
if (j+1<=i-1)
{
for (k=j+1k<=i-1k++)
{
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=0j<=i-1j++)
{
f=q[i*n+j]
g=c[j-1]-h2*f
c[j-1]=g
for (k=0k<=jk++)
{
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=0i<=n-1i++)
{
if ((b[i]!=0.0)&&(i-1>=0))
{
for (j=0j<=i-1j++)
{
g=0.0
for (k=0k<=i-1k++)
{
g=g+q[i*n+k]*q[k*n+j]
}
for (k=0k<=i-1k++)
{
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=0j<=i-1j++)
{
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=0j<=n-1j++)
{
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+1i<=n-1i++)
{
b[i]=b[i]-h
}
f=f+h
p=b[m]
e=1.0
s=0.0
for (i=m-1i>=ji--)
{
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=0k<=n-1k++)
{
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=0i<=n-1i++)
{
k=ip=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=0j<=n-1j++)
{
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=1k<=n-2k++)
{
d=0.0
for (j=kj<=n-1j++)
{
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-1j<=n-1j++)
{
u=i*n+j
v=k*n+j
t=a[u]
a[u]=a[v]
a[v]=t
}
for (j=0j<=n-1j++)
{
u=j*n+i
v=j*n+k
t=a[u]
a[u]=a[v]
a[v]=t
}
}
for (i=k+1i<=n-1i++)
{
u=i*n+k-1
t=a[u]/d
a[u]=0.0
for (j=kj<=n-1j++)
{
v=i*n+j
a[v]=a[v]-t*a[k*n+j]
}
for (j=0j<=n-1j++)
{
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-1it=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.0v[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+2j<=m-1j++)
{
a[j*n+j-2]=0.0
}
for (j=l+3j<=m-1j++)
{
a[j*n+j-3]=0.0
}
for (k=lk<=m-2k++)
{
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=kj<=m-1j++)
{
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=li<=ji++)
{
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=0i<=n-1i++)
{
v[i*n+i]=1.0
for (j=0j<=n-1j++)
{
if (i!=j)
{
v[i*n+j]=0.0
}
}
}
while (1==1)
{
fm=0.0
for (i=0i<=n-1i++)
{
for (j=0j<=n-1j++)
{
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=0j<=n-1j++)
{
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=0i<=n-1i++)
{
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=0i<=n-1i++)
{
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)
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）