CODEWITHSUNDEEP
c++sparse-matrixlinear-algebraeigen
Kyraz998
Kyraz998

Reputation: 23

Eigen sparse solver wrong results

I am trying to solve a sparse linear system Ax=B with the Eigen library in C++, however the following trivial example seems to give an incorrect solution:

#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\n" << A << "\n";
    cout << "B\n" << B << "\n";
    cout << "x\n" << x << "\n";

    return 0;
}

I don't see any errors, the algorithm returns "0" meaning "Success", however the solution I get is

x = 0 0.193524 0 0.193524 0.193524 0 0 0 0

which is obviously not the solution to this system, the correct one is

x = 0 0 0 0 0.0967621 0 0 0 0

Upvotes: 2

Views: 1326

Answers (1)

Soonts
Soonts

Reputation: 21936

Here's documentation for SimplicialLDLT solver:

This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are selfadjoint and positive definite.

When the matrix stores real numbers in the elements, self-adjoint == symmetrical. Your matrix is clearly not symmetrical. Also, not every symmetrical matrix is positive-definite, see examples.

In short, the solver you have chosen is only applicable to very narrow class of matrices. As you have already discovered, SparseLU solver works for your input data.

ConjugateGradient solver will not work either, it doesn’t require the matrix to be positive-definite but it does require it to be self-adjoint.

Upvotes: 4

Related Questions