3

I am using finite differences for a square 100x100 domain (with neumann bcs on all sides) in c++ using Eigen's sparse matrix functionality, and built in solvers to compute x in Ax=b.

I have tried the following solvers, but am getting very different time results to what I expect from reading the documentation in http://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html which gives typical timescales for different solvers. In particular, the documentation suggests that conjugate gradients should be one of the fastest ways of solving this system, giving a timescale of 0.239 seconds for solving Poisson SPD with a size bigger than my system. By contrast this documentation suggests that SimplicialLLT should take roughly 3X as long.

When I run each of the solvers I am obtaining the following: - Conjugate gradients: 25 seconds - LLT: 0.35 seconds

I was wondering whether anyone could help me understand why there is two orders of magnitude between these two solvers, and in particular, why CG seems to be being beaten by LLT, in contrast to the literature? Also, if anyone else has an idea how I could substantially speed up the solvers by using different methods from other packages, then suggestions are welcome!

I am implementing the solvers by the following:

//Conjugate gradients
ConjugateGradient<SparseMatrix<double> > cg;
cg.compute(A);
MatrixXd vGDNF = cg.solve(b);

//SimplicialLLT
SimplicialLLT<SparseMatrix<double>> solver;
MatrixXd vGDNF = solver.compute(A).solve(b);

Here A is that of a finite difference Laplacian, and b is a vector of point sources of the field.

Thanks!

Ben

ben18785
  • 356
  • 1
  • 5
  • 15

1 Answers1

4

On this documentation, Poisson_SPD corresponds to a 3D problem which is much harder for direct methods as Cholesky because there are much more fill-in in the factors. As the same documentation page says in the first table, for 2D Laplacian/Poisson problems as yours, SimplicialLLT is the recommended solver.

ggael
  • 28,425
  • 2
  • 65
  • 71