用C++编写matlab中eig(A,B)的函数，输出特征向量矩阵和特征值矩阵

以下是用C语言编写的类似于matlab中eig(A,B)函数的代码，实现了求解一般化特征值问题（Ax = λBx）的功能，并输出特征向量矩阵和特征值矩阵：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 3        // 矩阵A和B的维度
#define EPS 1e-10  // 判断精度

void mat_print(double A[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            printf("%.3f ", A[i][j]);
        }
        printf("\n");
    }
}

void mat_transpose(double A[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=i+1; j<n; j++) {
            double tmp = A[i][j];
            A[i][j] = A[j][i];
            A[j][i] = tmp;
        }
    }
}

void mat_sub(double A[][N], double B[][N], double C[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

void mat_mul(double A[][N], double B[][N], double C[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            C[i][j] = 0;
            for (int k=0; k<n; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

void mat_copy(double A[][N], double B[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            B[i][j] = A[i][j];
        }
    }
}

void mat_eye(double A[][N], int n) {
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            if (i == j) {
                A[i][j] = 1;
            } else {
                A[i][j] = 0;
            }
        }
    }
}

double mat_norm(double A[][N], int n) {
    double res = 0.0;
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            res += A[i][j] * A[i][j];
        }
    }
    return sqrt(res);
}

void qr_decomp(double A[][N], double Q[][N], double R[][N], int n) {
    double V[N][N], H[N][N], I[N][N];
    mat_eye(Q, n);
    mat_copy(A, R, n);
    for (int k=0; k<n-1; k++) {
        mat_eye(V, n);
        double normx = 0.0;
        for (int i=k; i<n; i++) {
            normx += R[i][k] * R[i][k];
        }
        double x1 = R[k][k] + (R[k][k] > 0 ? 1 : -1) * sqrt(normx);
        V[k][k] = x1;
        for (int i=k+1; i<n; i++) {
            V[i][k] = R[i][k];
        }
        mat_mul(V, V, H, n);
        double normh = mat_norm(H, n);
        mat_transpose(V, n);
        mat_copy(Q, I, n);
        double alpha = 2.0 / (normh * normh);
        double beta = -alpha * alpha;
        for (int i=0; i<n; i++) {
            for (int j=0; j<n; j++) {
                Q[i][j] = I[i][j] + alpha * H[i][j] + beta * mat_mul(Q, H, Q, n)[i][j];
            }
        }
        mat_mul(H, R, R, n);
    }
}

void eig(double A[][N], double B[][N], double V[][N], double D[][N], int n) {
    double C[N][N], Q[N][N], R[N][N];
    mat_mul(B, A, C, n);
    qr_decomp(C, Q, R, n);
    mat_mul(R, Q, C, n);
    int iter = 0;
    while (mat_norm(C, n) > EPS && iter < 1000) {
        iter++;
        mat_mul(B, C, D, n);
        qr_decomp(D, Q, R, n);
        mat_mul(R, Q, D, n);
        mat_mul(C, D, C, n);
    }
    for (int i=0; i<n; i++) {
        double lambda = C[i][i];
        double x[N], y[N];
        for (int j=0; j<n; j++) {
            x[j] = C[j][i];
            y[j] = B[j][j];
        }
        for (int j=0; j<n; j++) {
            V[j][i] = x[j] / sqrt(x[j] * x[j] + y[j] * y[j]);
            D[i][i] = lambda;
        }
    }
}

int main() {
    double A[N][N] = {{1, 2, 3}, {2, 5, 6}, {3, 6, 9}};
    double B[N][N] = {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}};
    double V[N][N], D[N][N];
    eig(A, B, V, D, N);
    printf("特征向量矩阵：\n");
    mat_print(V, N);
    printf("特征值矩阵：\n");
    mat_print(D, N);
    return 0;
}

这里给出的是一个比较简单的实现，可能不够高效或者不够健壮，但可以作为参考。需要注意的是，对于一般化特征值问题，矩阵B必须是对称正定的，否则算法可能无法收敛或者得到正确的结果。