Hello,
I've tried different methods of solving a sparse, symmetric system with eigen3 and I'm just surprised how bad the performance is, compared to matlab. To give you an example:
I have a 5149 x 5149 sparse, symmetric Matrix A and a (actual sparse but since the implemented solvers in eigen can just handle dense it's a dense) vector b. The solution vector x is also dense vector.
My first try was using implemented "SimplicialCholesky":

Code: Select all

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


(since I normally need the decomposition for more than one right side I used "compute()" and not directly "solve()", just an example here).
This code is very slow and these three lines need about 48 seconds to be computed.

So I tried cholmod:

Code: Select all

  Eigen::SparseLLT<Eigen::SparseMatrix<double>, Eigen::Cholmod> cholmoDec;
  cholmoDec.compute (A);
  x = cholmoDec.solve (b);


Here I was surprised about the good performance, in average it took 3.1 seconds to compute these lines.

I have also tried the latest svn and tried to use cg:

Code: Select all

  Eigen::ConjugateGradient<Eigen::SparseMatrix<double> > cg1;
  cg1.compute (A);
  x = cg1.solve (b);


Which took 35 seconds to complete, for smaller problems it was OK (still slower than cholmod but I realized lower memory usage).

I also tried umfpack:

Code: Select all

  Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::UmfPack> umfDec(A);
  umfDec.solve(b, &x);


Here it took 4 seconds, but a strange thing is that with bigger problems (I tried one with 13000 degrees of freedom) it just took about 6 seconds, cholmod needed about 30 seconds for this - and the problem was symmetric. Very strange in my opinion. Maybe someone can enlighten me?

But the really nasty thing is that matlab outperforms them all. With matlab and the ordinary backslash operator "\" it took in average 0.3 seconds to solve exactly this problem! 0.3 - unbelievable. I didn't get under 3 seconds with eigen3 and backends.
Can anybody give me advice here? Did I make something wrong or can I tune the performance somehow further?
I forgot to mention that I compiled my code using g++ (GCC) 4.6.1 with the options: "-Ofast -msse2 -fopenmp -DNDEBUG", I also tried "-O2" instead of "Ofast" with no luck.
So any hints?
Thanks
matse

Hi!

I was just about to ask the very same questions. Any updates?

Best!
Heiko

that's hard to help without more information. At least we would need the number of non zeros to have an estimate of sparse are your matrices. Also, which kind of problems are you trying to solve? Laplacian-like? bi-laplacian-like? 2D or 3D domain? circuit simulation? It would be nice if you could share your matrices in market format too.

anyway, 48s for a 5k x 5k problem is abnormally slow... your problems might be not sparse. I have example of 3D-Poisson problems without ~1M of unknown solved in a couple of seconds with CG, and 2D bi-harmonic problems with ~100k unknowns solved in about 0.1s with SimplicialLDLT.