jacobi方法求矩阵特征值和特征向量C语言代码

以下是Jacobi方法求解矩阵特征值和特征向量的C语言代码：

#include <stdio.h>
#include <math.h>
#define N 3

void matrix_print(double a[N][N])
{
    int i,j;

    for(i=0;i<N;i++){
        for(j=0;j<N;j++){
            printf("%lf ",a[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

void matrix_rot(double a[N][N], int p, int q, double theta)
{
    int i,j;
    double c = cos(theta);
    double s = sin(theta);
    double tau = (a[q][q]-a[p][p])/(2*a[p][q]);

    for(i=0;i<N;i++){
        if(i!=p && i!=q){
            double ap = a[i][p];
            double aq = a[i][q];
            a[i][p] = ap*c - aq*s;
            a[i][q] = ap*s + aq*c;
            a[p][i] = a[i][p];
            a[q][i] = a[i][q];
        }
    }
    a[p][p] = a[p][p]*c*c - 2*a[p][q]*c*s + a[q][q]*s*s;
    a[q][q] = a[p][p] + 2*a[p][q]*c*s + a[q][q]*c*c;
    a[p][q] = 0;
    a[q][p] = 0;
}

void jacobi(double a[N][N], double v[N][N])
{
    int i,j,p,q;
    double eps = 1e-8;
    double theta,t;
    double max;

    for(i=0;i<N;i++){
        for(j=0;j<N;j++){
            if(i==j) v[i][j] = 1.0;
            else v[i][j] = 0.0;
        }
    }
    while(1){
        max = -1.0;
        for(i=0;i<N;i++){
            for(j=0;j<i;j++){
                t = fabs(a[i][j]);
                if(t>max){
                    max = t;
                    p = i;
                    q = j;
                }
            }
        }
        if(max<eps) break;
        theta = 0.5*atan2(2*a[p][q],a[q][q]-a[p][p]);
        matrix_rot(a,p,q,theta);
        matrix_rot(v,p,q,theta);
    }
}

int main()
{
    double a[N][N] = {{4,-1,-1},{-1,3,-2},{-1,-2,3}};
    double v[N][N];
    jacobi(a,v);
    matrix_print(a);
    matrix_print(v);
    return 0;
}

以上代码中，`a`为待求特征值和特征向量的矩阵，`v`为输出的特征向量矩阵，`matrix_print`为矩阵输出函数，`matrix_rot`为旋转矩阵函数，`jacobi`为Jacobi方法主函数。在`jacobi`函数中，首先将特征向量矩阵初始化为单位矩阵，然后不断进行Jacobi旋转，以使矩阵对角化，最终得到特征值和特征向量。