数值计算程序大放送-特征值和特征向量

dormouse_buaa

<>选自&lt;&lt;徐世良数值计算程序集(C)&gt;&gt;</P>
<>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。</P>
<>今天都给贴出来了</P>
<P>#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&lt;=n-1; i++)
{
  for (j=0; j&lt;=n-1; j++)
  {
   u=i*n+j; q=a;
  }
}
    for (i=n-1; i&gt;=1; i--)
{
  h=0.0;
        if (i&gt;1)
  {
   for (k=0; k&lt;=i-1; k++)
            {
    u=i*n+k;
    h=h+q*q;
   }
  }
  
        if (h+1.0==1.0)
  {
   c[i-1]=0.0;
            if (i==1)
   {
    c[i-1]=q[i*n+i-1];
   }
            b=0.0;
  }
        else
  {
   c[i-1]=sqrt(h);
            u=i*n+i-1;
            if (q&gt;0.0)
   {
    c[i-1]=-c[i-1];
   }
            h=h-q*c[i-1];
            q=q-c[i-1];
            f=0.0;
            for (j=0; j&lt;=i-1; j++)
   {
    q[j*n+i]=q[i*n+j]/h;
                g=0.0;
                for (k=0; k&lt;=j; k++)
    {
     g=g+q[j*n+k]*q[i*n+k];
    }
                if (j+1&lt;=i-1)
    {
     for (k=j+1; k&lt;=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&lt;=i-1; j++)
   {
    f=q[i*n+j];
                g=c[j-1]-h2*f;
                c[j-1]=g;
                for (k=0; k&lt;=j; k++)
    {
     u=j*n+k;
                    q=q-f*c[k-1]-g*q[i*n+k];
    }
   }
            b=h;
  }
}
    b[0]=0.0;
    for (i=0; i&lt;=n-1; i++)
{
  if ((b!=0.0)&amp;&amp;(i-1&gt;=0))
  {
   for (j=0; j&lt;=i-1; j++)
            {
    g=0.0;
    for (k=0; k&lt;=i-1; k++)
    {
     g=g+q[i*n+k]*q[k*n+j];
    }
    for (k=0; k&lt;=i-1; k++)
                {
     u=k*n+j;
     q=q-g*q[k*n+i];
                }
            }
  }
        u=i*n+i;
        b=q;
  q=1.0;
        if (i-1&gt;=0)
  {
   for (j=0; j&lt;=i-1; j++)
            {
    q[i*n+j]=0.0;
    q[j*n+i]=0.0;
   }
  }
}
    return;
}
</P>

dormouse_buaa

//求实对称三对角对称矩阵的全部特征值及特征向量
//利用变型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&lt;=n-1; j++)
{
  it=0;
        h=eps*(fabs(b[j])+fabs(c[j]));
        if (h&gt;d)
  {
   d=h;
  }
        m=j;
        while ((m&lt;=n-1)&amp;&amp;(fabs(c[m])&gt;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&gt;=0.0)
    {
     b[j]=c[j]/(p+r);
    }
                else
    {
     b[j]=c[j]/(p-r);
    }
                h=g-b[j];
                for (i=j+1; i&lt;=n-1; i++)
    {
     b=b-h;
    }
                f=f+h;
    p=b[m];
    e=1.0;
    s=0.0;
                for (i=m-1; i&gt;=j; i--)
    {
     g=e*c;
     h=e*p;
                    if (fabs(p)&gt;=fabs(c))
     {
      e=c/p;
      r=sqrt(e*e+1.0);
                        c[i+1]=s*p*r;
      s=e/r;
      e=1.0/r;
     }
                    else
     {
      e=p/c;
      r=sqrt(e*e+1.0);
                        c[i+1]=s*c*r;
                        s=1.0/r;
      e=e/r;
     }
                    p=e*b-s*g;
                    b[i+1]=h+s*(e*g+s*b);
                    for (k=0; k&lt;=n-1; k++)
     {
      u=k*n+i+1;
      v=u-1;
                        h=q;
      q=s*q[v]+e*h;
                        q[v]=e*q[v]-s*h;
     }
    }
                c[j]=s*p;
    b[j]=e*p;
   }
            while (fabs(c[j])&gt;d);
  }
        b[j]=b[j]+f;
}
    for (i=0; i&lt;=n-1; i++)
{
  k=i; p=b;
        if (i+1&lt;=n-1)
  {
   j=i+1;
            while ((j&lt;=n-1)&amp;&amp;(b[j]&lt;=p))
   {
    k=j;
    p=b[j];
    j=j+1;
   }
  }
        if (k!=i)
  {
   b[k]=b;
   b=p;
            for (j=0; j&lt;=n-1; j++)
   {
    u=j*n+i;
    v=j*n+k;
                p=q;
    q=q[v];
    q[v]=p;
   }
  }
}
    return(1);
}

dormouse_buaa

//约化实矩阵为赫申伯格(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&lt;=n-2; k++)
{
  d=0.0;
        for (j=k; j&lt;=n-1; j++)
  {
   u=j*n+k-1;
   t=a;
            if (fabs(t)&gt;fabs(d))
   {
    d=t;
    i=j;
   }
  }
        if (fabs(d)+1.0!=1.0)
  {
   if (i!=k)
   {
    for (j=k-1; j&lt;=n-1; j++)
    {
     u=i*n+j;
     v=k*n+j;
                    t=a;
     a=a[v];
     a[v]=t;
    }
                for (j=0; j&lt;=n-1; j++)
    {
     u=j*n+i;
     v=j*n+k;
                    t=a;
     a=a[v];
     a[v]=t;
    }
   }
            for (i=k+1; i&lt;=n-1; i++)
   {
    u=i*n+k-1;
    t=a/d;
    a=0.0;
                for (j=k; j&lt;=n-1; j++)
    {
     v=i*n+j;
                    a[v]=a[v]-t*a[k*n+j];
    }
                for (j=0; j&lt;=n-1; j++)
    {
     v=j*n+k;
                    a[v]=a[v]+t*a[j*n+i];
    }
   }
  }
}
    return;
}

dormouse_buaa

//求赫申伯格(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&gt;0)&amp;&amp;(fabs(a[l*n+l-1])&gt;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&gt;0.0)
   {
    xy=1.0;
                if (b&lt;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&gt;=jt)
   {
    printf("fail\n");
                return(-1);
   }
            it=it+1;
            for (j=l+2; j&lt;=m-1; j++)
   {
    a[j*n+j-2]=0.0;
   }
            for (j=l+3; j&lt;=m-1; j++)
   {
    a[j*n+j-3]=0.0;
   }
            for (k=l; k&lt;=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&lt;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&lt;=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&gt;=m-1)
     {
      j=m-1;
     }
                    for (i=l; i&lt;=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);
}

dormouse_buaa

<>//求实对称矩阵的特征值及特征向量的雅格比法
//利用雅格比(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&lt;=n-1; i++)
{
  v[i*n+i]=1.0;
        for (j=0; j&lt;=n-1; j++)
  {
   if (i!=j)
   {
    v[i*n+j]=0.0;
   }
  }
}
    while (1==1)
{
  fm=0.0;
        for (i=0; i&lt;=n-1; i++)
  {
   for (j=0; j&lt;=n-1; j++)
   {
    d=fabs(a[i*n+j]);
    if ((i!=j)&amp;&amp;(d&gt;fm))
    {
     fm=d;
     p=i;
     q=j;
    }
   }
  }
        if (fm&lt;eps)  
  {
   return(1);
  }
        if (l&gt;jt)  
  {
   return(-1);
  }
        l=l+1;
        u=p*n+q;
  w=p*n+p;
  t=q*n+p;
  s=q*n+q;
        x=-a;
  y=(a-a[w])/2.0;
        omega=x/sqrt(x*x+y*y);
        if (y&lt;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*sn*sn+a*omega;
        a=fm*sn*sn+a*cn*cn-a*omega;
        a=0.0;
  a[t]=0.0;
        for (j=0; j&lt;=n-1; j++)
  {
   if ((j!=p)&amp;&amp;(j!=q))
   {
    u=p*n+j;
    w=q*n+j;
    fm=a;
    a=fm*cn+a[w]*sn;
    a[w]=-fm*sn+a[w]*cn;
   }
  }
        for (i=0; i&lt;=n-1; i++)
  {
   if ((i!=p)&amp;&amp;(i!=q))
            {
    u=i*n+p;
    w=i*n+q;
    fm=a;
    a=fm*cn+a[w]*sn;
    a[w]=-fm*sn+a[w]*cn;
            }
  }
        for (i=0; i&lt;=n-1; i++)
  {
   u=i*n+p;
   w=i*n+q;
            fm=v;
            v=fm*cn+v[w]*sn;
            v[w]=-fm*sn+v[w]*cn;
  }
}
return(1);
}</P>

dormouse_buaa

//////////////////////////////////////////////////////////////
//约化对称矩阵为三对角对称矩阵
//利用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);
//////////////////////////////////////////////////////////////