Fast solve a sparse positive definite linear system with Eigen3

Question

I need to solve a linear system with a sparse symmetric and positive definite matrix (2630x2630) millions of times. I have ploted the matrix in Mathematica an it's shown bellow.

I have chosen the Eigen3 lib with the LLT decomposition to solve the linear system, which compared to others methods like LU is much faster. The system solution took 0.385894 seconds in a intel 10700 with 4.8 GHz processor. Code:

#include <Eigen/Core>
using namespace Eigen;
using namespace std;


int main()
{
    ifstream datamatrix("matrix.txt");
    ifstream datavec("vector.txt");
    int m=2630;
    MatrixXd A(m,m);
    VectorXd b(m);
    for (int i = 0; i < m; i++)
    {
     for (int j = 0; j < m; j++)
     {
         datamatrix >> A(i,j);
     }
     datavec >> b(i);
    }
    chrono::steady_clock sc;
    auto start = sc.now();
    VectorXd x = A.llt().solve(b);
    auto end = sc.now();
    // measure time span between start & end
    auto time_span = static_cast<chrono::duration<double>>(end - start);
    cout << "Operation took: " << time_span.count() << " seconds.";
}

Is it possible to accelerate this with Eigen3 or MKL?

The matrix and vector files can be downloaded here:

Matrix: https://www.dropbox.com/s/k521t91cd8u7t5h/matrix.txt?dl=0

Vector: https://www.dropbox.com/s/ldajnzl2qj3a7zh/vector.txt?dl=0

EDIT

I found that the majority of the time spent was in filling the sparse matrix. With the code bellow using Triplet to fill the sparse matrix it is much faster! obs. MatDoub is the full-dense matrix from numerical recipes code - nr3.h

    void slopeproject::SolveEigenSparse(MatDoub A, MatDoub b, MatDoub& x)
    {
    
        typedef Eigen::Triplet<double> T;
        std::vector<T> tripletList;
        int sz=A.nrows();
    
//approximated number of non-zero entries in the matrix
        tripletList.reserve(sz*50);
       // tripletList.reserve(80000);
    
        x.assign(sz, 1, 0.);
        SparseMatrix<double> AA(sz, sz);
        VectorXd bbb(sz);
        for (int i = 0; i < sz; i++)
        {
            for (int j = 0; j < sz; j++)
            {
//checking if the value of the dense matrix is zero. If not it is appended to the tripletlist

                if(fabs(A[i][j])>1.e-12)
                {
                    tripletList.push_back(T(i,j,A[i][j]));
                }
            }
            bbb(i) = b[i][0];
        }
    //transfer from the tripletlist to the sparse matrix very fast

        AA.setFromTriplets(tripletList.begin(), tripletList.end());
    //optional. I dont know wath this function do.
        AA.makeCompressed();
        SimplicialLLT< SparseMatrix<double> > solver;
//solve the system
        VectorXd xx = solver.compute(AA).solve(bbb);
        for(int i=0;i<sz;i++)x[i][0]=xx(i);
    }

Answer 1

You have a sparse matrix, but you're representing it in Eigen as a dense matrix. The matrix file that you have is also dense, it would be more convenient to use if it was stored in sparse form, the Market format for example.

If I change the matrix to a sparse one, and use

SimplicialLLT< SparseMatrix<double> > solver;
VectorXd x = solver.compute(A).solve(b);

Then on my PC (which should be slower than yours, it has only an old 4770K) the actual factor/solve only takes 0.01 seconds.

By the way even the dense solve should be faster than what you saw, so I guess that you did not enable SIMD in your compiler settings.