First of all, thanks for creating such a great library, it is very easy to use and has almost everything I need.

I recently did some speed tests comparing Intel MKL PARDISO and Eigen for solving large sparse symmetric matrices arising from FEM for elasticity problems. The sizes of the matrices range from 30k to 300k. PARDISO is configured to use single thread only since Eigen does not support multi-threaded LLT (Eigen::SimplicialLLT). I found that the decomposition time of PARDISO is 10x to 30x times faster than Eigen, and the speed gap grows with the size of the matrix.

I know PARDISO is highly optimized by Intel people, but I did not expect it to be thus faster than Eigen if they are using the same algorithm. Is it because PARDISO uses more efficient algorithms in finding a better permutation and doing the decomposition?

Another interesting thing I found is that although Eigen is much slower in doing decomposition, its solve time always outperforms PARDISO, and can be 3x faster for some test cases.

I guess you are dealing with a 3D domain, then a 10x factor is expected because SimplicialLDLT only implements a simple algorithm which does not extract dense blocks (in contrast to supernodal/multifrontal approahes). In this case Eigen::SparseLU should be faster though it does not exploit the symmetry. ConjugateGradient should even be faster than PARDISO. See this table for basic recommendation:

http://eigen.tuxfamily.org/dox-devel/gr ... tml#title0

and see the row "poisson_SPD" for some timings on a Poisson problem discretised on a 2D tet mesh: http://eigen.tuxfamily.org/dox-devel/gr ... tml#title7