C++ Jacobi Algorithm Class for Eigenvalue and Eigenvector Calculation

#include #include

class Jacobi { public: // 输入对称矩阵A的维度和元素值 Jacobi(int n, const std::vector& A) : n(n), A(A), V(nn, 0.0) { for (int i = 0; i < n; i++) { V[in+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 A; // 对称矩阵A std::vector 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

int main() { // 定义对称矩阵A int n = 3; std::vector A{nn}; A[0n+0] = 4.0; A[0n+1] = 1.0; A[0n+2] = 1.0; A[1n+0] = 1.0; A[1n+1] = 3.0; A[1n+2] = 1.0; A[2n+0] = 1.0; A[2n+1] = 1.0; A[2n+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;