Large sparse linear system solving with SimplicialCholesky Eigen

Question

I am actually trying to solve large sparse linear systems using C++ lib Eigen. Sparse matrices are taken from this page. Each system as this structure: Ax = b where A is the sparse matrix (n x n), b is computed as A*xe with xe vector of dimension n containing only zeros. After computing x I need to compute the relative error between xe and x. I have written some code but I can't understand why the relative error is so high (1.49853e+08) at the end of the computation.

#include <iostream>
#include <Eigen/Dense>
#include <unsupported/Eigen/SparseExtra>
#include<Eigen/SparseCholesky>
#include <sys/time.h>
#include <sys/resource.h>


using namespace std;
using namespace Eigen;


int main()
{

    SparseMatrix<double> mat;
    loadMarket(mat, "/Users/anto/Downloads/ex15/ex15.mtx");

	VectorXd xe = VectorXd::Constant(mat.rows(), 1);
	VectorXd b = mat*xe;

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

    double relative_error = (x-xe).norm()/(xe).norm(); 
    cout << relative_error << endl;
    
}

The matrix ex15 can be downloaded from this page. It is a symmetric, positive definite matrix. Can anyone help me to solve the problem? Thank you in advance for your help.

Have you tested this on a much simpler matric with a known, expected result? — Scott Hunter
Hi Scott, ex15 matrix is a very simple sparse matrix. I use as benchmark Matlab's solution for the same system. — Anto

ggael ggael · Accepted Answer · 2019-04-09T19:54:37

According to this page, ex15 is not full rank. You should check that each step went well:

SimplicialLDLT<Eigen::SparseMatrix<double> > chol(mat);
if(chol.info()!=Eigen::Success)
  return;
VectorXd x = chol.solve(b); 
if(chol.info()!=Eigen::Success)
  return;

and then check that you got one solution (if it's not full rank and that at least one solution exists, then there exist a whole subspace of solutions):

cout << (mat*x-b).norm()/b.norm() << "\n";

Large sparse linear system solving with SimplicialCholesky Eigen

1 Answers