（急求）关于求解小型三对角矩阵特征值和特征向量的问题

daxia

我正在做流固耦合问题，已经把大型非对称矩阵降阶为小型三对角矩阵，下一步就是求解三对角矩阵的特征值和特征向量的问题。这个问题已经困扰我好长时间了，希望各位能帮小弟一把，不胜感激！！！

luckyfox

网络上或者期刊网上应该有类似的文章吧？
对于对称的三对角阵，这里有一个从别的地方得来的源程序。不过本人没有做过这方面的内容。

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

liuhaibo

这是清华出的一本书《C语言常用算法集》上的源程序。请问哪位大哥能把它注释一下不？本人在原理上不能将其联系起来。不胜感谢！！！如果能帮助本人解决此问题，本人愿献上一些最近的研究成果。谢谢！liuhaibo20061985@126.com