特征稀疏解算器错误结果

2022-07-19 00:00:00 sparse-matrix linear-algebra c++ eigen

我正在尝试用C++的本征库来解决一个稀疏线性系统Ax=B,但是下面这个简单的例子似乎给出了一个不正确的解决方案:

#include <Eigen/SparseCholesky>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <iostream>
#include <vector>

using namespace std;
using namespace Eigen;

int main(){

    SimplicialLDLT<SparseMatrix<double>> solver;
    SparseMatrix<double> A(9,9);
    typedef Triplet<double> T;
    vector<T> triplet;
    VectorXd B(9);

    for(int i=0; i<4; i++){
        triplet.push_back(T(i,i,1));
        triplet.push_back(T(i+5,i+5,1));
    }

    triplet.push_back(T(4,1,-1));
    triplet.push_back(T(4,3,-1));
    triplet.push_back(T(4,5,-1));
    triplet.push_back(T(4,7,-1));
    triplet.push_back(T(4,4,4));

    A.setFromTriplets(triplet.begin(),triplet.end());
    B << 0,0,0,0,0.387049,0,0,0,0;

    solver.compute(A);
    VectorXd x = solver.solve(B);

    cout << "A
" << A << "
";
    cout << "B
" << B << "
";
    cout << "x
" << x << "
";

    return 0;
}

我看不到任何错误,算法返回&0&q;表示成功,但我得到的解决方案是

x = 0 0.193524 0 0.193524 0.193524 0 0 0 0

这显然不是此系统的解决方案,正确的解决方案是

x = 0 0 0 0 0.0967621 0 0 0 0

解决方案

Here's documentationSimplicialLDLT求解器:

此类提供了一个不含稀疏矩阵平方根的LDL^T Cholesky分解稀疏矩阵是自伴且是正定的。

当矩阵在元素中存储实数时,自伴==对称。你的矩阵显然不是对称的。而且,并不是每个对称矩阵都是正定的see examples。

简而言之,您选择的求解器仅适用于非常窄的矩阵类。正如您已经发现的,SparseLU求解器适用于您的输入数据。

ConjugateGradient求解器也不工作,它不要求矩阵是正定的,但it does要求它是自伴的。

相关文章