8
$\begingroup$

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++.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Did you try Intel MKL (software.intel.com/en-us/mkl)? If I'm not wrong, it should be free for Linux. $\endgroup$
    – Tolga Birdal
    Jan 26, 2019 at 22:46
  • $\begingroup$ 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. $\endgroup$
    – KOF
    Jan 27, 2019 at 7:55
  • 1
    $\begingroup$ 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}$ :-) $\endgroup$
    – Wolfgang Bangerth
    Jan 28, 2019 at 16:37
  • $\begingroup$ @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$. $\endgroup$
    – Nox
    Jan 29, 2019 at 16:04
  • 1
    $\begingroup$ 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. $\endgroup$
    – rchilton1980
    Jan 29, 2019 at 16:07

1 Answer 1

Reset to default
10
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – KOF
    Feb 2, 2019 at 12:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nox
    Feb 2, 2019 at 18:27
  • 1
    $\begingroup$ github.com/flame/blis is based on goto-blas, with less hand-written assembler but still good performance. $\endgroup$
    – Mr. White
    Oct 25, 2019 at 11:54
  • 1
    $\begingroup$ Does the suggestion to use Eigen with OpenBLAS really apply to sparse matrices? AFAICT eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html says BLAS or LAPACK are used only for dense matrix products and dense matrix decompositions. $\endgroup$
    – Andrei Matveiakin
    Aug 30, 2021 at 14:23
  • $\begingroup$ @AndreiMatveiakin No, and that was never implied. Eigen3 allows one to use BLAS routines for BLAS operations which may, or may not, be used during sparse matrix operatinos (e.g. they would be if one were to use Krylov decomposition for Eigenvalues). As for BLIS: That is a very nice BLAS library, which actually uses a performance model for optimization. I, too, can highly recommend it! $\endgroup$
    – Nox
    Nov 8, 2021 at 21:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.