# Which C++ linear algebra library is probably the fastest on solving huge sparse [square matrix] linear system?

I am developing a 2D CFD solver for fluid-particle interaction. To solve Navier-Stokes equations on a grid of size $$10000\times 10000$$ cells (or >1 million cells), a large linear system $$Ax=b$$ with $$A$$ being the $$10000\times10000$$ sparse coefficient matrix needs to be solved efficiently in each time-step.

What I am looking for is a high-performance/parallel C++ linear algebra library to solve this large sparse linear system. An iterative approach such as biconjugate gradient stabilized method is preferred.

There are a lot of existing libraries out there like: Eigen3, PETSc, Trilinos, MLT4, GNU GSL, Armadillo, LAPACK++, and the list goes on.

Among the well-known libraries, which one should I choose for my project, in terms of high-performance (better with OpenMP/MPI support), and easy-to-use in c++.

• Did you try Intel MKL (software.intel.com/en-us/mkl)? If I'm not wrong, it should be free for Linux. – Tolga Birdal Jan 26 at 22:46
• I had experience on Eigen3, but I am not sure the performance on large sparse linear system pointed out by some posts. I will try PETSc and Eigen3 for a comparison. – KOF Jan 27 at 7:55
• I'll note that for most practitioners, 10,000 x 10,000 is actually a pretty tiny sparse matrix these days, not a huge one. Huge these days would be, say, $10^{10}\times 10^{10}$ :-) – Wolfgang Bangerth Jan 28 at 16:37
• @WolfgangBangerth I concur. 10k x 10k is something you can pack as a dense matrix on a consumer-grade GPU. Even 500k x 500k is still moderate. Largest I've used in research (linear algebra methods) was $10^7\times 10^7$. – Nox Jan 29 at 16:04
• Although the OP is a little confusing, I think it might be discussing a 2D grid with 10K elements on a side, so more in the neighborhood of (10K)^2 = 10^8 unknowns. – rchilton1980 Jan 29 at 16:07

Eigen 3 is a nice C++ template library some of whose routines are parallelized. c.f. Eigen documentation The parallelization is OMP only, so if you intend to parallelise using MPI (and OMP) it is probably not suitable for your purpose. The nice feature of Eigen is that you can swap in a high performance BLAS library (like MKL or OpenBLAS) for some routines by simply using #define EIGEN_USE_BLAS (and other macros).

Similarly Armadillo allows for node-level parallelism only. In my experience it is better to use Eigen since it is easier to interface with the raw C++ arrays in Eigen, which facilitates use of other libraries (e.g. ARPACK++).

In my experience I would advise against using GSL for linear algebra. I have found its performance to be lacking and the usability to be worse than that of Eigen.

If you plan to execute linear solvers (e.g. BiCGstab) on multiple nodes I would advise you to use Trilinos. I have used it in my research codes with fairly little delving into its documentation due to the good examples available on the Trilinos homepage. Furthermore its performance is decent and can be fine-tuned by including good BLAS/LAPACK libraries during the compilation.

Similar should hold for PETSc, although I have never actively used its LA routines. In my experience PETSc is dependency hell, if you want to use performance-optimized (e.g. optimized for the CPU architecture you're using) versions of the libraries it requires. The performance should be fairly decent, too, I think, since PETSc relies on common LA libraries (BLAS, LAPACK, SCALAPACK etc.).

Long story short: For interoperability and good performance on a single node (using OpenMP) I advise to use Eigen with OpenBLAS. If you want to use multiple nodes via MPI and let the library figure out how to solve a system using multiple nodes then use Trilinos.

• thanks very much for the detailed information. I will try Eigen 3 first to see if it can solve the sparse linear system (c.a 1 million cells) within 5 seconds per time-step. – KOF Feb 2 at 12:37
• I would've liked to give a few empiric results, but our cluster is currently full. If you want to use Eigen it is imperative that you compile with all necessary optimization flags enabled (see Eigen documentation). When using OpenBLAS as BLAS provider it must be optimized, this is less of a Problem when using Intel MKL as BLAS provider for Eigen. – Nox Feb 2 at 18:27