Eigen
3.3.90 (mercurial changeset 8071cda5714d)

In Eigen, there are several methods available to solve linear systems when the coefficient matrix is sparse. Because of the special representation of this class of matrices, special care should be taken in order to get a good performance. See Sparse matrix manipulations for a detailed introduction about sparse matrices in Eigen. This page lists the sparse solvers available in Eigen. The main steps that are common to all these linear solvers are introduced as well. Depending on the properties of the matrix, the desired accuracy, the enduser is able to tune those steps in order to improve the performance of its code. Note that it is not required to know deeply what's hiding behind these steps: the last section presents a benchmark routine that can be easily used to get an insight on the performance of all the available solvers.
Eigen currently provides a wide set of builtin solvers, as well as wrappers to external solver libraries. They are summarized in the following tables:
Class  Solver kind  Matrix kind  Features related to performance  License 
Notes 

SimplicialLLT #include<Eigen/SparseCholesky>  Direct LLt factorization  SPD  Fillin reducing  LGPL  SimplicialLDLT is often preferable 
SimplicialLDLT #include<Eigen/SparseCholesky>  Direct LDLt factorization  SPD  Fillin reducing  LGPL  Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.) 
SparseLU #include<Eigen/SparseLU>  LU factorization  Square  Fillin reducing, Leverage fast dense algebra  MPL2  optimized for small and large problems with irregular patterns 
SparseQR #include<Eigen/SparseQR>  QR factorization  Any, rectangular  Fillin reducing  MPL2  recommended for leastsquare problems, has a basic rankrevealing feature 
Class  Solver kind  Matrix kind  Supported preconditioners, [default]  License 
Notes 

ConjugateGradient #include<Eigen/IterativeLinearSolvers>  Classic iterative CG  SPD  IdentityPreconditioner, [DiagonalPreconditioner], IncompleteCholesky  MPL2  Recommended for large symmetric problems (e.g., 3D Poisson eq.) 
LeastSquaresConjugateGradient #include<Eigen/IterativeLinearSolvers>  CG for rectangular leastsquare problem  Rectangular  IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]  MPL2  Solve for min A'Axb^2 without forming A'A 
BiCGSTAB #include<Eigen/IterativeLinearSolvers>  Iterative stabilized biconjugate gradient  Square  IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT  MPL2  To speedup the convergence, try it with the IncompleteLUT preconditioner. 
Class  Module  Solver kind  Matrix kind  Features related to performance  Dependencies,License 
Notes 

PastixLLT PastixLDLT PastixLU  PaStiXSupport  Direct LLt, LDLt, LU factorizations  SPD SPD Square  Fillin reducing, Leverage fast dense algebra, Multithreading  Requires the PaStiX package, CeCILLC  optimized for tough problems and symmetric patterns 
CholmodSupernodalLLT  CholmodSupport  Direct LLt factorization  SPD  Fillin reducing, Leverage fast dense algebra  Requires the SuiteSparse package, GPL  
UmfPackLU  UmfPackSupport  Direct LU factorization  Square  Fillin reducing, Leverage fast dense algebra  Requires the SuiteSparse package, GPL  
KLU  KLUSupport  Direct LU factorization  Square  Fillin reducing, suitted for circuit simulation  Requires the SuiteSparse package, GPL  
SuperLU  SuperLUSupport  Direct LU factorization  Square  Fillin reducing, Leverage fast dense algebra  Requires the SuperLU library, (BSDlike)  
SPQR  SPQRSupport  QR factorization  Any, rectangular  fillin reducing, multithreaded, fast dense algebra  requires the SuiteSparse package, GPL  recommended for linear leastsquares problems, has a rankrevealing feature 
PardisoLLT PardisoLDLT PardisoLU  PardisoSupport  Direct LLt, LDLt, LU factorizations  SPD SPD Square  Fillin reducing, Leverage fast dense algebra, Multithreading  Requires the Intel MKL package, Proprietary  optimized for tough problems patterns, see also using MKL with Eigen 
Here SPD
means symmetric positive definite.
All these solvers follow the same general concept. Here is a typical and general example:
For SPD
solvers, a second optional template argument allows to specify which triangular part have to be used, e.g.:
In the above example, only the upper triangular part of the input matrix A is considered for solving. The opposite triangle might either be empty or contain arbitrary values.
In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:
The compute() method is equivalent to calling both analyzePattern() and factorize().
Each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on. More details are available in the documentations of the respective classes.
Finally, most of the iterative solvers, can also be used in a matrixfree context, see the following example .
In the compute() function, the matrix is generally factorized: LLT for selfadjoint matrices, LDLT for general hermitian matrices, LU for non hermitian matrices and QR for rectangular matrices. These are the results of using direct solvers. For this class of solvers precisely, the compute step is further subdivided into analyzePattern() and factorize().
The goal of analyzePattern() is to reorder the nonzero elements of the matrix, such that the factorization step creates less fillin. This step exploits only the structure of the matrix. Hence, the results of this step can be used for other linear systems where the matrix has the same structure. Note however that sometimes, some external solvers (like SuperLU) require that the values of the matrix are set in this step, for instance to equilibrate the rows and columns of the matrix. In this situation, the results of this step should not be used with other matrices.
Eigen provides a limited set of methods to reorder the matrix in this step, either builtin (COLAMD, AMD) or external (METIS). These methods are set in template parameter list of the solver :
See the OrderingMethods module for the list of available methods and the associated options.
In factorize(), the factors of the coefficient matrix are computed. This step should be called each time the values of the matrix change. However, the structural pattern of the matrix should not change between multiple calls.
For iterative solvers, the compute step is used to eventually setup a preconditioner. For instance, with the ILUT preconditioner, the incomplete factors L and U are computed in this step. Remember that, basically, the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear system where the coefficient matrix has more clustered eigenvalues. For real problems, an iterative solver should always be used with a preconditioner. In Eigen, a preconditioner is selected by simply adding it as a template parameter to the iterative solver object.
The member function preconditioner() returns a readwrite reference to the preconditioner to directly interact with it. See the Iterative solvers module and the documentation of each class for the list of available methods.
The solve() function computes the solution of the linear systems with one or many right hand sides.
Here, B can be a vector or a matrix where the columns form the different right hand sides. The solve() function can be called several times as well, for instance when all the right hand sides are not available at once.
For direct methods, the solution are computed at the machine precision. Sometimes, the solution need not be too accurate. In this case, the iterative methods are more suitable and the desired accuracy can be set before the solve step using setTolerance(). For all the available functions, please, refer to the documentation of the Iterative solvers module .
Most of the time, all you need is to know how much time it will take to solve your system, and hopefully, what is the most suitable solver. In Eigen, we provide a benchmark routine that can be used for this purpose. It is very easy to use. In the build directory, navigate to bench/spbench and compile the routine by typing make spbenchsolver. Run it with –help option to get the list of all available options. Basically, the matrices to test should be in MatrixMarket Coordinate format, and the routine returns the statistics from all available solvers in Eigen.
To export your matrices and righthandside vectors in the matrixmarket format, you can the the unsupported SparseExtra module:
The following table gives an example of XML statistics from several Eigen builtin and external solvers.
Matrix  N  NNZ  UMFPACK  SUPERLU  PASTIX LU  BiCGSTAB  BiCGSTAB+ILUT  GMRES+ILUT  LDLT  CHOLMOD LDLT  PASTIX LDLT  LLT  CHOLMOD SP LLT  CHOLMOD LLT  PASTIX LLT  CG  

vector_graphics  12855  72069  Compute Time  0.0254549  0.0215677  0.0701827  0.000153388  0.0140107  0.0153709  0.0101601  0.00930502  0.0649689  
Solve Time  0.00337835  0.000951826  0.00484373  0.0374886  0.0046445  0.00847754  0.000541813  0.000293696  0.00485376  
Total Time  0.0288333  0.0225195  0.0750265  0.037642  0.0186552  0.0238484  0.0107019  0.00959871  0.0698227  
Error(Iter)  1.299e16  2.04207e16  4.83393e15  3.94856e11 (80)  1.03861e12 (3)  5.81088e14 (6)  1.97578e16  1.83927e16  4.24115e15  
poisson_SPD  19788  308232  Compute Time  0.425026  1.82378  0.617367  0.000478921  1.34001  1.33471  0.796419  0.857573  0.473007  0.814826  0.184719  0.861555  0.470559  0.000458188 
Solve Time  0.0280053  0.0194402  0.0268747  0.249437  0.0548444  0.0926991  0.00850204  0.0053171  0.0258932  0.00874603  0.00578155  0.00530361  0.0248942  0.239093  
Total Time  0.453031  1.84322  0.644241  0.249916  1.39486  1.42741  0.804921  0.862891  0.4989  0.823572  0.190501  0.866859  0.495453  0.239551  
Error(Iter)  4.67146e16  1.068e15  1.3397e15  6.29233e11 (201)  3.68527e11 (6)  3.3168e15 (16)  1.86376e15  1.31518e16  1.42593e15  3.45361e15  3.14575e16  2.21723e15  7.21058e16  9.06435e12 (261)  
sherman2  1080  23094  Compute Time  0.00631754  0.015052  0.0247514    0.0214425  0.0217988  
Solve Time  0.000478424  0.000337998  0.0010291    0.00243152  0.00246152  
Total Time  0.00679597  0.01539  0.0257805    0.023874  0.0242603  
Error(Iter)  1.83099e15  8.19351e15  2.625e14  1.3678e+69 (1080)  4.1911e12 (7)  5.0299e13 (12)  
bcsstk01_SPD  48  400  Compute Time  0.000169079  0.00010789  0.000572538  1.425e06  9.1612e05  8.3985e05  5.6489e05  7.0913e05  0.000468251  5.7389e05  8.0212e05  5.8394e05  0.000463017  1.333e06 
Solve Time  1.2288e05  1.1124e05  0.000286387  8.5896e05  1.6381e05  1.6984e05  3.095e06  4.115e06  0.000325438  3.504e06  7.369e06  3.454e06  0.000294095  6.0516e05  
Total Time  0.000181367  0.000119014  0.000858925  8.7321e05  0.000107993  0.000100969  5.9584e05  7.5028e05  0.000793689  6.0893e05  8.7581e05  6.1848e05  0.000757112  6.1849e05  
Error(Iter)  1.03474e16  2.23046e16  2.01273e16  4.87455e07 (48)  1.03553e16 (2)  3.55965e16 (2)  2.48189e16  1.88808e16  1.97976e16  2.37248e16  1.82701e16  2.71474e16  2.11322e16  3.547e09 (48)  
sherman1  1000  3750  Compute Time  0.00228805  0.00209231  0.00528268  9.846e06  0.00163522  0.00162155  0.000789259  0.000804495  0.00438269  
Solve Time  0.000213788  9.7983e05  0.000938831  0.00629835  0.000361764  0.00078794  4.3989e05  2.5331e05  0.000917166  
Total Time  0.00250184  0.00219029  0.00622151  0.0063082  0.00199698  0.00240949  0.000833248  0.000829826  0.00529986  
Error(Iter)  1.16839e16  2.25968e16  2.59116e16  3.76779e11 (248)  4.13343e11 (4)  2.22347e14 (10)  2.05861e16  1.83555e16  1.02917e15  
young1c  841  4089  Compute Time  0.00235843  0.00217228  0.00568075  1.2735e05  0.00264866  0.00258236  
Solve Time  0.000329599  0.000168634  0.00080118  0.0534738  0.00187193  0.00450211  
Total Time  0.00268803  0.00234091  0.00648193  0.0534865  0.00452059  0.00708447  
Error(Iter)  1.27029e16  2.81321e16  5.0492e15  8.0507e11 (706)  3.00447e12 (8)  1.46532e12 (16)  
mhd1280b  1280  22778  Compute Time  0.00234898  0.00207079  0.00570918  2.5976e05  0.00302563  0.00298036  0.00144525  0.000919922  0.00426444  
Solve Time  0.00103392  0.000211911  0.00105  0.0110432  0.000628287  0.00392089  0.000138303  6.2446e05  0.00097564  
Total Time  0.0033829  0.0022827  0.00675918  0.0110692  0.00365392  0.00690124  0.00158355  0.000982368  0.00524008  
Error(Iter)  1.32953e16  3.08646e16  6.734e16  8.83132e11 (40)  1.51153e16 (1)  6.08556e16 (8)  1.89264e16  1.97477e16  6.68126e09  
crashbasis  160000  1750416  Compute Time  3.2019  5.7892  15.7573  0.00383515  3.1006  3.09921  
Solve Time  0.261915  0.106225  0.402141  1.49089  0.24888  0.443673  
Total Time  3.46381  5.89542  16.1594  1.49473  3.34948  3.54288  
Error(Iter)  1.76348e16  4.58395e16  1.67982e14  8.64144e11 (61)  8.5996e12 (2)  6.04042e14 (5) 