求矩阵特征值的vc代码

int cJcbEigen(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)
    {
        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);
}