c++ eigen MatrixXd崩溃

Question

c++ Eigen库实现SOR法于与直接法解不同维度线性方程组时间对比, MatrixXd 来定义A,x,b 求解Ax=b.但是MatrixXd A(SIZE, SIZE)在SIZE=20以上时程序就崩溃,没有具体报错.是因为MatrixXd double内存太大吗?还是怎么回事?

#include <iostream>
#include <cmath>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Dense>
#include <vector> 
#include "matplotlibcpp.h"
#include <locale>  
#include <chrono>


using namespace std;
using namespace Eigen;
namespace plt = matplotlibcpp;
void plot(vector<double> ii, vector<double> ww, vector<double> jj);
void plot(vector<vector<double>> vctime);

int main(int argc, char **argv) {

    vector<vector<double>> vctime;
    // 系数矩阵大小 1-30
    for(int SIZE = 1; SIZE <= 30; SIZE ++){
        int i, j;
        MatrixXd A(SIZE, SIZE);
        MatrixXd b(SIZE, 1);
        std::srand(std::time(0)); //保证每次运行随机生成的A,b于上次不同
        A = A.setRandom();    
        b = b.setRandom();
        for (i = 0; i < A.cols(); i++){ 
            A(i, i) = A(i, i) + 3;    
        }
        cout << "\nA = :\n " << A << endl;
        cout << "b = : " << b.transpose() << endl;
        if( A.determinant() == 0){
            cout << "A 不可逆" << endl;
            return 0;
        }
        
        // 判断杜利特尔分解条件是否成立
        double k = 1;
        for (i = 1; i <= A.cols(); i++ ){
            k = k * A.block(0,0,i,i).determinant();
        }

        if (k == 0){
            cout << "T 矩阵顺序主子式都不等零 不成立!" << endl;
            return 0;
        }
        
        // 杜利特尔分解
        // 修改U,L不会改变T
        MatrixXd U(SIZE, SIZE) ;
        U = A.triangularView<Upper>(); 
        MatrixXd L(SIZE, SIZE)  ;
        L = A.triangularView<Lower>(); 
        MatrixXd D(SIZE, SIZE)  ;
        D = D.setZero();
        U = - U;
        L = - L;
        for (i = 0; i < U.cols(); i++){
            L(i, i) = 0;
            U(i, i) = 0;
            D(i, i) = A(i, i);
        }
        /*cout << "上三角矩阵 U :\n" << U << endl ;
        cout << "下三角矩阵 L :\n" << L << endl ;
        cout << "对角矩阵 D :\n" << D << endl ;*/
        
        // 计算谱半径 元素是std::complex<double, double> 复数!
        MatrixXd I(SIZE, SIZE);
        I = I.setIdentity();
        MatrixXd J(SIZE, SIZE); 
        J = I - D.inverse() * A;
        VectorXcd _eigenvalues = J.eigenvalues();  
        // 复数特征值取模长
        VectorXd eigenvalues = _eigenvalues.cwiseAbs();
        //cout << "谱半径 : "<< eigenvalues.maxCoeff() << "雅各比矩阵 J :\n" << J << endl;
        
        double w;
        double error = 1e-5;
        vector<double> cost;
        vector<double> ii;
        vector<double> ww;
        vector<double> tt;
        vector<double> jj;
        vector<MatrixXd> y;
        
        //SOR法
        double step = 0.05;
        for (w = 0 + step ; w < 1; w = w + step){
            MatrixXd B(SIZE, SIZE);
            B = (D - w * L).inverse() * ((1 - w) * D + w * U);
            MatrixXd f(SIZE, 1);
            f = w * (D - w * L).inverse() * b;
            MatrixXd x(SIZE, 1);
            x = x.setOnes();// 初始化为1
            
            VectorXcd _eigenvalues = B.eigenvalues();  
            // 复数特征值取模长
            VectorXd eigenvalues = _eigenvalues.cwiseAbs();
            // 谱半径
            double max_eigenvalues = eigenvalues.maxCoeff();
            
            for (i = 0; i < 100; i++){
                chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
                
                MatrixXd temp(SIZE, 1);
                temp = x;
                x = B * x + f;
                double _cost = abs((x - temp).norm());

                if (_cost <= error){
                    /*cout << "ρ(B) = " << max_eigenvalues << endl; 
                    cout << "i = " << i << endl;
                    cout << "x = " << x.transpose() << endl;*/
                    chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
                    chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
                    tt.push_back(time_used.count());
                    cost.push_back(_cost); 
                    ii.push_back(i); 
                    ww.push_back(w);
                    y.push_back(x);
                    jj.push_back(max_eigenvalues);
                    break;
                    }
            }
        }
        // 最佳 ω 和索引
        double i_min = *std::min_element(ii.begin(), ii.end());
        int index = std::distance(ii.begin(), std::min_element(ii.begin(), ii.end())); 
        cout << "最佳松弛因子 ω = " << ww.at(index) <<" 最小迭代次数 i = " << i_min << " \nSOR迭代时间 t = " << tt.at(index) 
        << "s SOR迭代结果 x = \n" << y.at(index).transpose() << endl;
        
        // 直接法
        chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
        MatrixXd dx(SIZE, 1);
        dx = A.inverse() * b;
        chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
        chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
        cout << "直接法用时 t = " << time_used.count() << "s 直接法结果 x =\n" << dx.transpose(); 
        
        vctime.push_back({double(SIZE), tt.at(index), time_used.count()});
    }    
    //plot(ii, ww, jj);
    plot(vctime);
}
void plot(vector<double> ii, vector<double> ww, vector<double> jj){
    // ω, ρ
    plt::subplot(3, 1, 1);
    plt::plot(ww,ii);
    plt::xlabel("$\\omega$");
    plt::ylabel("iteration(s)");
    
    plt::subplot(3, 1, 2);
    plt::plot(ww,jj);
    plt::xlabel("$\\omega$");
    plt::ylabel("$\\rho(B_\\omega)$");
    
    //convergence rate
    vector<double> cr;
    for (auto a : jj){
        cr.push_back(- log(a));
    }
    
    plt::subplot(3, 1, 3);
    plt::plot(ww,cr);
    plt::xlabel("$\\omega$");
    plt::ylabel("convergence rate");
    
    plt::show();
}
void plot(vector<vector<double>> vctime){
    // ω, ρ
    vector<double>  NN;
    vector<double>  SORtime;
    vector<double>  Directtime;
    for(int i = 0; i < vctime.size(); i++){
        NN.push_back(vctime.at(i).at(0));
        SORtime.push_back(vctime.at(i).at(1));
        Directtime.push_back(vctime.at(i).at(2));
    }
    plt::named_plot("SOR", NN, SORtime, "-r");
    plt::named_plot("Direct", NN, Directtime, "-b");
    plt::xlabel("$\\omega$");
    plt::ylabel("iteration time ");
    plt::legend();
    plt::show();
}

Answer 1

【以下回答由 GPT 生成】

问题分析：

根据您提供的代码，当矩阵大小SIZE大于20时，程序会崩溃。根据您的描述，崩溃可能是由于MatrixXd的内存过大导致的，或者可能由于其他原因引起。为了解决这个问题，我将采取以下步骤：

1. 检查当前程序崩溃的具体原因。
2. 优化矩阵的内存管理。
3. 检查其他可能导致程序崩溃的原因。

步骤1：检查当前程序崩溃的具体原因：

为了检查程序崩溃的具体原因，我们可以使用try-catch块来捕获异常并获取更多的错误信息。

try {
    // 进行Ax=b的求解操作
} catch (const std::exception& e) {
    std::cerr << "Exception caught: " << e.what() << std::endl;
    return 1;
}

通过添加这个try-catch块，任何引发的异常都将被捕获，并将异常信息输出到标准错误流。这样我们就可以知道具体导致程序崩溃的原因。

步骤2：优化矩阵的内存管理：

为了解决内存过大导致的崩溃问题，可以考虑使用动态分配内存。在这种情况下，可以使用new运算符创建动态分配的矩阵和向量。

MatrixXd* A = new MatrixXd(SIZE, SIZE);
MatrixXd* b = new MatrixXd(SIZE, 1);

然后，在进行完线性方程求解操作后，记得释放动态分配的内存。

delete A;
delete b;

这样可以避免在每次循环迭代中重复分配和释放内存，从而提高程序的性能。

步骤3：检查其他可能导致程序崩溃的原因：

如果在优化了内存管理之后程序仍然崩溃，那么可能存在其他原因导致崩溃。以下是一些可能导致程序崩溃的其他原因：

• 内存溢出：如果矩阵A和向量b的大小太大，可能会导致内存溢出。可以尝试减小矩阵的大小或者增加程序的内存限制。
• 数值问题：如果矩阵A的条件数很高，可能会导致数值稳定性问题，从而导致程序崩溃。可以考虑使用数值稳定性更好的求解方法，如QR分解或SVD分解。
• 编译器问题：在某些情况下，程序崩溃可能是由于编译器的bug或不兼容性引起的。可以尝试使用其他的C++编译器或更新编译器版本来解决这个问题。

这些是一些可能导致程序崩溃的其他原因，可能需要进一步的调试和优化来解决这些问题。

根据以上几点，我已经给出了一些优化建议，请您根据我的建议尝试改进代码。如果还有问题，请告诉我，我将尽力帮助解决。

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^