Solving Systems of Linear Equations using Eigen

Question

I'm currently working on a fluid simulation in C++, and part of the algorithm is to solve a sparse system of linear equations. People recommended using the library Eigen for this. I decided to test it out using this short program that I wrote:

#include <Eigen/SparseCholesky>

#include <vector>
#include <iostream>

int main() {
    std::vector<Eigen::Triplet<double>> triplets;
    triplets.push_back(Eigen::Triplet<double>(0, 0, 1));
    triplets.push_back(Eigen::Triplet<double>(0, 1, -2));
    triplets.push_back(Eigen::Triplet<double>(1, 0, 3));
    triplets.push_back(Eigen::Triplet<double>(1, 1, -2));

    Eigen::SparseMatrix<double> A(2, 2);
    A.setFromTriplets(triplets.begin(), triplets.end());

    Eigen::VectorXd b(2);
    b[0] = -2;
    b[1] = 2;

    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
    Eigen::VectorXd x = chol.solve(b);

    std::cout << x[0] << ' ' << x[1] << std::endl;

    system("pause");
}

It gives it these two equations:

x - 2y = -2

3x - 2y = 2

The correct solution is:

x = 2

y = 2

But the problem is that when the program runs, it outputs: 0.181818 -0.727273

Which is totally wrong! I have been debugging this for hours, but it's a very short program and I'm following the tutorial on the Eigen website exactly. Does anybody know what is causing this issue?

P.S. I know that the classes I'm using are for sparse matrices, but the only difference between those and the normal Matrix classes is the way the elements are stored.

@n.m. I thought that might be it, but even when I transpose it I still get the wrong answer :/ — paper man
@n.m. What are you saying? These answers still don’t work even when I transpose the coefficient matrix — paper man
Hmm no sorry I messed up my math. This is a solution not for a transposed matrix, but for (1,3,3,-2). — n. 1.8e9-where's-my-share m.
Which is only logical, as simplicial Cholesky only works for self-adjoint matrices! It is also deprecated in Eigen. — n. 1.8e9-where's-my-share m.

ggael ggael · Accepted Answer · 2018-09-04T07:15:26

SimplicialCholesky is for symmetric positive definite (SPD) matrices, the matrix you assembled is not even symmetric. By default it only reads the entries in the lower triangular part ignoring the others, so it solved:

x + 3y = -2
3x -2y = 2

As you noticed, for non-symmetric square problems you need to use a direct solver based on LU or BICGSTAB in the world of iterative solvers. This is all summarized in the doc.

Solving Systems of Linear Equations using Eigen

2 Answers