cjcbi.c

这是用Jacob方法计算实对称矩阵特征根及特征向量的程序

//////////////////////////////////////////////////////////////////
//用雅可比方法求n阶对称矩阵A的特征值及相应特征向量的步骤如下:
//(1)令S＝In(In为单位矩阵)
//(2)在A中选取非对角线元家中绝对值最大者，设为apq
//(3)若fabs(apq)<ε则迭代过程结束。此时对角线元素aii(i＝0，1，…，M—1)即为特征值λi。矩阵S的第i列为与λi对应的特征向量。否则，继续下一步。
//(4)计算平面旋转矩阵的元素及其变换后的矩阵Al的元素。
//用jacobi方法计算实对称矩阵的特征值和特征向量
//a——双精度实型二维数组，体积为n*n。存放n阶实对称矩阵A；返回时，对角线上存放n个特征值。
//n-整型变量。实对称矩阵A的阶数。
//V--双精度实型二维数组，体积为n*n。返回特征向量。其中第i列为与λi(即返回的aii＝0，1，…，M—1)对应的特征向量。
 // eps--双精度实型变量。控制精度要求。
// jt——整型变量。控制最大迭代次数。
//本函数返回一个整型标志值。若返回的标志值小于0，则表示迭代了Jt次还达不到精度要求；若返回的标志值大于o，则表示正常返回。
///////////////////////////////////////////
  #include "math.h"
  int cjcbi(a,n,v,eps,jt)
  int n,jt;
  double a[],v[],eps;
  { 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=1; i<=n-1; i++)
        for (j=0; j<=i-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);
  }