写一个Jacobi算法的c++类库

以下是一个简单的Jacobi算法C++类库，用于计算对称矩阵的特征值和特征向量。

#include <cmath>
#include <vector>

class Jacobi {
public:
    // 输入对称矩阵A的维度和元素值
    Jacobi(int n, const std::vector<double>& A) : n(n), A(A), V(n*n, 0.0) {
        for (int i = 0; i < n; i++) {
            V[i*n+i] = 1.0;
        }
    }

    // 进行Jacobi迭代计算特征值和特征向量
    void solve(std::vector<double>& eigenvalues, std::vector<double>& eigenvectors) {
        // 计算初始矩阵的最大非对角元素
        double maxoff = maxoffdiag();
        // 迭代直到所有非对角元素都小于给定的阈值
        while (maxoff > 1.0e-10) {
            // 找到最大非对角元素的位置
            int p = 0, q = 0;
            findmaxoffdiag(p, q);
            // 计算旋转角度
            double theta = calc_theta(p, q);
            // 更新矩阵A和V
            update_matrix(p, q, theta);
            // 更新最大非对角元素
            maxoff = maxoffdiag();
        }
        // 将对角线元素存入eigenvalues中
        eigenvalues.resize(n);
        for (int i = 0; i < n; i++) {
            eigenvalues[i] = A[i*n+i];
        }
        // 将V矩阵的列向量存入eigenvectors中
        eigenvectors.resize(n*n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                eigenvectors[j*n+i] = V[i*n+j];
            }
        }
    }

private:
    int n;  // 矩阵维度
    std::vector<double> A;  // 对称矩阵A
    std::vector<double> V;  // 特征向量矩阵V

    // 计算矩阵A的最大非对角元素
    double maxoffdiag() {
        double maxoff = 0.0;
        for (int i = 0; i < n; i++) {
            for (int j = i+1; j < n; j++) {
                double aij = std::abs(A[i*n+j]);
                if (aij > maxoff) {
                    maxoff = aij;
                }
            }
        }
        return maxoff;
    }

    // 找到矩阵A的最大非对角元素的位置
    void findmaxoffdiag(int& p, int& q) {
        double maxoff = 0.0;
        for (int i = 0; i < n; i++) {
            for (int j = i+1; j < n; j++) {
                double aij = std::abs(A[i*n+j]);
                if (aij > maxoff) {
                    maxoff = aij;
                    p = i;
                    q = j;
                }
            }
        }
    }

    // 计算旋转角度
    double calc_theta(int p, int q) {
        double app = A[p*n+p];
        double aqq = A[q*n+q];
        double apq = A[p*n+q];
        return 0.5 * std::atan2(2*apq, aqq-app);
    }

    // 更新矩阵A和V
    void update_matrix(int p, int q, double theta) {
        double c = std::cos(theta);
        double s = std::sin(theta);
        for (int k = 0; k < n; k++) {
            if (k != p && k != q) {
                double apk = A[p*n+k];
                double aqk = A[q*n+k];
                A[p*n+k] = apk*c - aqk*s;
                A[q*n+k] = aqk*c + apk*s;
                A[k*n+p] = A[p*n+k];
                A[k*n+q] = A[q*n+k];
            }
        }
        double app = A[p*n+p];
        double aqq = A[q*n+q];
        double apq = A[p*n+q];
        A[p*n+p] = app*c*c - 2*apq*c*s + aqq*s*s;
        A[q*n+q] = aqq*c*c + 2*apq*c*s + app*s*s;
        A[p*n+q] = (aqq-app)*c*s + apq*(c*c-s*s);
        A[q*n+p] = A[p*n+q];
        for (int k = 0; k < n; k++) {
            double vkp = V[k*n+p];
            double vkq = V[k*n+q];
            V[k*n+p] = vkp*c - vkq*s;
            V[k*n+q] = vkq*c + vkp*s;
        }
    }
};

使用示例：

#include <iostream>

int main() {
    // 定义对称矩阵A
    int n = 3;
    std::vector<double> A{n*n};
    A[0*n+0] = 4.0; A[0*n+1] = 1.0; A[0*n+2] = 1.0;
    A[1*n+0] = 1.0; A[1*n+1] = 3.0; A[1*n+2] = 1.0;
    A[2*n+0] = 1.0; A[2*n+1] = 1.0; A[2*n+2] = 5.0;

    // 计算特征值和特征向量
    Jacobi jacobi(n, A);
    std::vector<double> eigenvalues;
    std::vector<double> eigenvectors;
    jacobi.solve(eigenvalues, eigenvectors);

    // 输出结果
    std::cout << "Eigenvalues: ";
    for (double lambda : eigenvalues) {
        std::cout << lambda << " ";
    }
    std::cout << std::endl;

    std::cout << "Eigenvectors: " << std::endl;
    for (int i = 0; i < n; i++) {
        std::cout << "[";
        for (int j = 0; j < n; j++) {
            std::cout << eigenvectors[j*n+i] << " ";
        }
        std::cout << "]" << std::endl;
    }

    return 0;
}
``