namespace Eigen { /** \page TopicSparseSystems Solving Sparse Linear Systems 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 \ref TutorialSparse for a detailed introduction about sparse matrices in Eigen. In this page, we briefly present the main steps that are common to all the linear solvers in Eigen together with the main concepts behind them. Depending on the properties of the matrix, the desired accuracy, the end-user is able to tune these steps in order to improve the performance of its code. However, an impatient user does not need 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. \b Table \b of \b contents \n - \ref TheSparseCompute - \ref TheSparseSolve - \ref BenchmarkRoutine As summarized in \ref TutorialSparseDirectSolvers, there are many built-in solvers in Eigen as well as interface to external solvers libraries. All these solvers follow the same calling sequence. The basic steps are as follows : \code #include // ... SparseMatrix A; // fill A VectorXd b, x; // fill b // solve Ax = b SolverClassName > solver; solver.compute(A); if(solver.info()!=Succeeded) { // decomposition failed return; } x = solver.solve(b); if(solver.info()!=Succeeded) { // solving failed return; } \endcode \section TheSparseCompute The Compute Step In the compute() function, the matrix is generally factorized: LLT for self-adjoint matrices, LDLT for general hermitian matrices and LU for non hermitian 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 fill-in. 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 can note be used with other matrices. Eigen provides a limited set of methods to reorder the matrix in this step, either built-in (COLAMD, AMD) or external (METIS). These methods are set in template parameter list of the solver : \code DirectSolverClassName, OrderingMethod > solver; \endcode See \link Ordering_Modules the Ordering module \endlink 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. 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. \code IterativeSolverClassName, PreconditionerName > solver; \endcode The member function preconditioner() returns a read-write reference to the preconditioner to directly interact with it. For instance, with the ILUT preconditioner, the incomplete factors L and U are computed in this step. See \link Sparse_modules the Sparse module \endlink for the list of available preconditioners in Eigen. \section TheSparseSolve The Solve step The solve() function computes the solution of the linear systems with one or many right hand sides. \code X = solver.solve(B); \endcode 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. \code x1 = solver.solve(b1); // Get the second right hand side b2 x2 = solver.solve(b2); // ... \endcode 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 \link IterativeLinearSolvers_module Iterative solvers module \endlink. \section BenchmarkRoutine Most of the time, all you need is to know how much time it will take to qolve 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. First, it should be activated at the configuration step with the flag TEST_REAL_CASES. Then, in bench/spbench, you can compile the routine by typing \b make \e spbenchsolver. You can then run it with --help option to get the list of all available options. Basically, the matrices to test should be in \link http://math.nist.gov/MatrixMarket/formats.html MatrixMarket Coordinate format \endlink, and the routine returns the statistics from all available solvers in Eigen. The following table gives an example of XHTML statistics from several Eigen built-in 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.02545490.02156770.07018270.0001533880.01401070.01537090.01016010.009305020.0649689
Solve Time 0.003378350.0009518260.004843730.03748860.00464450.008477540.0005418130.0002936960.00485376
Total Time 0.02883330.02251950.07502650.0376420.01865520.02384840.01070190.009598710.0698227
Error(Iter) 1.299e-16 2.04207e-16 4.83393e-15 3.94856e-11 (80) 1.03861e-12 (3) 5.81088e-14 (6) 1.97578e-16 1.83927e-16 4.24115e-15
poisson_SPD 19788 308232 Compute Time 0.4250261.823780.6173670.0004789211.340011.334710.7964190.8575730.4730070.8148260.1847190.8615550.4705590.000458188
Solve Time 0.02800530.01944020.02687470.2494370.05484440.09269910.008502040.00531710.02589320.008746030.005781550.005303610.02489420.239093
Total Time 0.4530311.843220.6442410.2499161.394861.427410.8049210.8628910.49890.8235720.1905010.8668590.4954530.239551
Error(Iter) 4.67146e-16 1.068e-15 1.3397e-15 6.29233e-11 (201) 3.68527e-11 (6) 3.3168e-15 (16) 1.86376e-15 1.31518e-16 1.42593e-15 3.45361e-15 3.14575e-16 2.21723e-15 7.21058e-16 9.06435e-12 (261)
sherman2 1080 23094 Compute Time 0.006317540.0150520.0247514 -0.02144250.0217988
Solve Time 0.0004784240.0003379980.0010291 -0.002431520.00246152
Total Time 0.006795970.015390.0257805 -0.0238740.0242603
Error(Iter) 1.83099e-15 8.19351e-15 2.625e-14 1.3678e+69 (1080) 4.1911e-12 (7) 5.0299e-13 (12)
bcsstk01_SPD 48 400 Compute Time 0.0001690790.000107890.0005725381.425e-069.1612e-058.3985e-055.6489e-057.0913e-050.0004682515.7389e-058.0212e-055.8394e-050.0004630171.333e-06
Solve Time 1.2288e-051.1124e-050.0002863878.5896e-051.6381e-051.6984e-053.095e-064.115e-060.0003254383.504e-067.369e-063.454e-060.0002940956.0516e-05
Total Time 0.0001813670.0001190140.0008589258.7321e-050.0001079930.0001009695.9584e-057.5028e-050.0007936896.0893e-058.7581e-056.1848e-050.0007571126.1849e-05
Error(Iter) 1.03474e-16 2.23046e-16 2.01273e-16 4.87455e-07 (48) 1.03553e-16 (2) 3.55965e-16 (2) 2.48189e-16 1.88808e-16 1.97976e-16 2.37248e-16 1.82701e-16 2.71474e-16 2.11322e-16 3.547e-09 (48)
sherman1 1000 3750 Compute Time 0.002288050.002092310.005282689.846e-060.001635220.001621550.0007892590.0008044950.00438269
Solve Time 0.0002137889.7983e-050.0009388310.006298350.0003617640.000787944.3989e-052.5331e-050.000917166
Total Time 0.002501840.002190290.006221510.00630820.001996980.002409490.0008332480.0008298260.00529986
Error(Iter) 1.16839e-16 2.25968e-16 2.59116e-16 3.76779e-11 (248) 4.13343e-11 (4) 2.22347e-14 (10) 2.05861e-16 1.83555e-16 1.02917e-15
young1c 841 4089 Compute Time 0.002358430.002172280.005680751.2735e-050.002648660.00258236
Solve Time 0.0003295990.0001686340.000801180.05347380.001871930.00450211
Total Time 0.002688030.002340910.006481930.05348650.004520590.00708447
Error(Iter) 1.27029e-16 2.81321e-16 5.0492e-15 8.0507e-11 (706) 3.00447e-12 (8) 1.46532e-12 (16)
mhd1280b 1280 22778 Compute Time 0.002348980.002070790.005709182.5976e-050.003025630.002980360.001445250.0009199220.00426444
Solve Time 0.001033920.0002119110.001050.01104320.0006282870.003920890.0001383036.2446e-050.00097564
Total Time 0.00338290.00228270.006759180.01106920.003653920.006901240.001583550.0009823680.00524008
Error(Iter) 1.32953e-16 3.08646e-16 6.734e-16 8.83132e-11 (40) 1.51153e-16 (1) 6.08556e-16 (8) 1.89264e-16 1.97477e-16 6.68126e-09
crashbasis 160000 1750416 Compute Time 3.20195.789215.75730.003835153.10063.09921
Solve Time 0.2619150.1062250.4021411.490890.248880.443673
Total Time 3.463815.8954216.15941.494733.349483.54288
Error(Iter) 1.76348e-16 4.58395e-16 1.67982e-14 8.64144e-11 (61) 8.5996e-12 (2) 6.04042e-14 (5)
*/ }