source: pacpussensors/trunk/Vislab/lib3dv-1.2.0/lib3dv/eigen/doc/SparseLinearSystems.dox

Last change on this file was 136, checked in by ldecherf, 8 years ago

Doc

File size: 17.0 KB
Line 
1namespace Eigen {
2/** \eigenManualPage TopicSparseSystems Solving Sparse Linear Systems
3In 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. 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 end-user 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.
4
5\eigenAutoToc
6
7\section TutorialSparseDirectSolvers Sparse solvers
8
9%Eigen currently provides a limited set of built-in solvers, as well as wrappers to external solver libraries.
10They are summarized in the following table:
11
12<table class="manual">
13<tr><th>Class</th><th>Module</th><th>Solver kind</th><th>Matrix kind</th><th>Features related to performance</th>
14 <th>Dependencies,License</th><th class="width20em"><p>Notes</p></th></tr>
15<tr><td>SimplicialLLT </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
16 <td>built-in, LGPL</td>
17 <td>SimplicialLDLT is often preferable</td></tr>
18<tr><td>SimplicialLDLT </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LDLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
19 <td>built-in, LGPL</td>
20 <td>Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)</td></tr>
21<tr><td>ConjugateGradient</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Classic iterative CG</td><td>SPD</td><td>Preconditionning</td>
22 <td>built-in, MPL2</td>
23 <td>Recommended for large symmetric problems (e.g., 3D Poisson eq.)</td></tr>
24<tr><td>BiCGSTAB</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Iterative stabilized bi-conjugate gradient</td><td>Square</td><td>Preconditionning</td>
25 <td>built-in, MPL2</td>
26 <td>To speedup the convergence, try it with the \ref IncompleteLUT preconditioner.</td></tr>
27<tr><td>SparseLU</td> <td>\link SparseLU_Module SparseLU \endlink </td> <td>LU factorization </td>
28 <td>Square </td><td>Fill-in reducing, Leverage fast dense algebra</td>
29 <td> built-in, MPL2</td> <td>optimized for small and large problems with irregular patterns </td></tr>
30<tr><td>SparseQR</td> <td>\link SparseQR_Module SparseQR \endlink</td> <td> QR factorization</td>
31 <td>Any, rectangular</td><td> Fill-in reducing</td>
32 <td>built-in, MPL2</td><td>recommended for least-square problems, has a basic rank-revealing feature</td></tr>
33<tr> <th colspan="7"> Wrappers to external solvers </th></tr>
34<tr><td>PastixLLT \n PastixLDLT \n PastixLU</td><td>\link PaStiXSupport_Module PaStiXSupport \endlink</td><td>Direct LLt, LDLt, LU factorizations</td><td>SPD \n SPD \n Square</td><td>Fill-in reducing, Leverage fast dense algebra, Multithreading</td>
35 <td>Requires the <a href="http://pastix.gforge.inria.fr">PaStiX</a> package, \b CeCILL-C </td>
36 <td>optimized for tough problems and symmetric patterns</td></tr>
37<tr><td>CholmodSupernodalLLT</td><td>\link CholmodSupport_Module CholmodSupport \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing, Leverage fast dense algebra</td>
38 <td>Requires the <a href="http://www.suitesparse.com">SuiteSparse</a> package, \b GPL </td>
39 <td></td></tr>
40<tr><td>UmfPackLU</td><td>\link UmfPackSupport_Module UmfPackSupport \endlink</td><td>Direct LU factorization</td><td>Square</td><td>Fill-in reducing, Leverage fast dense algebra</td>
41 <td>Requires the <a href="http://www.suitesparse.com">SuiteSparse</a> package, \b GPL </td>
42 <td></td></tr>
43<tr><td>SuperLU</td><td>\link SuperLUSupport_Module SuperLUSupport \endlink</td><td>Direct LU factorization</td><td>Square</td><td>Fill-in reducing, Leverage fast dense algebra</td>
44 <td>Requires the <a href="http://crd-legacy.lbl.gov/~xiaoye/SuperLU/">SuperLU</a> library, (BSD-like)</td>
45 <td></td></tr>
46<tr><td>SPQR</td><td>\link SPQRSupport_Module SPQRSupport \endlink </td> <td> QR factorization </td>
47 <td> Any, rectangular</td><td>fill-in reducing, multithreaded, fast dense algebra</td>
48 <td> requires the <a href="http://www.suitesparse.com">SuiteSparse</a> package, \b GPL </td><td>recommended for linear least-squares problems, has a rank-revealing feature</tr>
49</table>
50
51Here \c SPD means symmetric positive definite.
52
53All these solvers follow the same general concept.
54Here is a typical and general example:
55\code
56#include <Eigen/RequiredModuleName>
57// ...
58SparseMatrix<double> A;
59// fill A
60VectorXd b, x;
61// fill b
62// solve Ax = b
63SolverClassName<SparseMatrix<double> > solver;
64solver.compute(A);
65if(solver.info()!=Success) {
66 // decomposition failed
67 return;
68}
69x = solver.solve(b);
70if(solver.info()!=Success) {
71 // solving failed
72 return;
73}
74// solve for another right hand side:
75x1 = solver.solve(b1);
76\endcode
77
78For \c SPD solvers, a second optional template argument allows to specify which triangular part have to be used, e.g.:
79
80\code
81#include <Eigen/IterativeLinearSolvers>
82
83ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
84x = solver.compute(A).solve(b);
85\endcode
86In 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.
87
88In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:
89\code
90SolverClassName<SparseMatrix<double> > solver;
91solver.analyzePattern(A); // for this step the numerical values of A are not used
92solver.factorize(A);
93x1 = solver.solve(b1);
94x2 = solver.solve(b2);
95...
96A = ...; // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
97solver.factorize(A);
98x1 = solver.solve(b1);
99x2 = solver.solve(b2);
100...
101\endcode
102The compute() method is equivalent to calling both analyzePattern() and factorize().
103
104Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on.
105More details are availble in the documentations of the respective classes.
106
107\section TheSparseCompute The Compute Step
108In the compute() function, the matrix is generally factorized: LLT for self-adjoint 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().
109
110The 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 should not be used with other matrices.
111
112Eigen 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 :
113\code
114DirectSolverClassName<SparseMatrix<double>, OrderingMethod<IndexType> > solver;
115\endcode
116
117See the \link OrderingMethods_Module OrderingMethods module \endlink for the list of available methods and the associated options.
118
119In 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.
120
121For 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.
122\code
123IterativeSolverClassName<SparseMatrix<double>, PreconditionerName<SparseMatrix<double> > solver;
124\endcode
125The member function preconditioner() returns a read-write reference to the preconditioner
126 to directly interact with it. See the \link IterativeLinearSolvers_Module Iterative solvers module \endlink and the documentation of each class for the list of available methods.
127
128\section TheSparseSolve The Solve step
129The solve() function computes the solution of the linear systems with one or many right hand sides.
130\code
131X = solver.solve(B);
132\endcode
133Here, 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.
134\code
135x1 = solver.solve(b1);
136// Get the second right hand side b2
137x2 = solver.solve(b2);
138// ...
139\endcode
140For 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 \b setTolerance(). For all the available functions, please, refer to the documentation of the \link IterativeLinearSolvers_Module Iterative solvers module \endlink.
141
142\section BenchmarkRoutine
143Most 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. In the build directory, navigate to bench/spbench and compile the routine by typing \b make \e spbenchsolver. Run it with --help option to get the list of all available options. Basically, the matrices to test should be in <a href="http://math.nist.gov/MatrixMarket/formats.html">MatrixMarket Coordinate format</a>, and the routine returns the statistics from all available solvers in Eigen.
144
145The following table gives an example of XML statistics from several Eigen built-in and external solvers.
146<TABLE border="1">
147 <TR><TH>Matrix <TH> N <TH> NNZ <TH> <TH > UMFPACK <TH > SUPERLU <TH > PASTIX LU <TH >BiCGSTAB <TH > BiCGSTAB+ILUT <TH >GMRES+ILUT<TH > LDLT <TH> CHOLMOD LDLT <TH > PASTIX LDLT <TH > LLT <TH > CHOLMOD SP LLT <TH > CHOLMOD LLT <TH > PASTIX LLT <TH> CG</TR>
148<TR><TH rowspan="4">vector_graphics <TD rowspan="4"> 12855 <TD rowspan="4"> 72069 <TH>Compute Time <TD>0.0254549<TD>0.0215677<TD>0.0701827<TD>0.000153388<TD>0.0140107<TD>0.0153709<TD>0.0101601<TD style="background-color:red">0.00930502<TD>0.0649689
149<TR><TH>Solve Time <TD>0.00337835<TD>0.000951826<TD>0.00484373<TD>0.0374886<TD>0.0046445<TD>0.00847754<TD>0.000541813<TD style="background-color:red">0.000293696<TD>0.00485376
150<TR><TH>Total Time <TD>0.0288333<TD>0.0225195<TD>0.0750265<TD>0.037642<TD>0.0186552<TD>0.0238484<TD>0.0107019<TD style="background-color:red">0.00959871<TD>0.0698227
151<TR><TH>Error(Iter) <TD> 1.299e-16 <TD> 2.04207e-16 <TD> 4.83393e-15 <TD> 3.94856e-11 (80) <TD> 1.03861e-12 (3) <TD> 5.81088e-14 (6) <TD> 1.97578e-16 <TD> 1.83927e-16 <TD> 4.24115e-15
152<TR><TH rowspan="4">poisson_SPD <TD rowspan="4"> 19788 <TD rowspan="4"> 308232 <TH>Compute Time <TD>0.425026<TD>1.82378<TD>0.617367<TD>0.000478921<TD>1.34001<TD>1.33471<TD>0.796419<TD>0.857573<TD>0.473007<TD>0.814826<TD style="background-color:red">0.184719<TD>0.861555<TD>0.470559<TD>0.000458188
153<TR><TH>Solve Time <TD>0.0280053<TD>0.0194402<TD>0.0268747<TD>0.249437<TD>0.0548444<TD>0.0926991<TD>0.00850204<TD>0.0053171<TD>0.0258932<TD>0.00874603<TD style="background-color:red">0.00578155<TD>0.00530361<TD>0.0248942<TD>0.239093
154<TR><TH>Total Time <TD>0.453031<TD>1.84322<TD>0.644241<TD>0.249916<TD>1.39486<TD>1.42741<TD>0.804921<TD>0.862891<TD>0.4989<TD>0.823572<TD style="background-color:red">0.190501<TD>0.866859<TD>0.495453<TD>0.239551
155<TR><TH>Error(Iter) <TD> 4.67146e-16 <TD> 1.068e-15 <TD> 1.3397e-15 <TD> 6.29233e-11 (201) <TD> 3.68527e-11 (6) <TD> 3.3168e-15 (16) <TD> 1.86376e-15 <TD> 1.31518e-16 <TD> 1.42593e-15 <TD> 3.45361e-15 <TD> 3.14575e-16 <TD> 2.21723e-15 <TD> 7.21058e-16 <TD> 9.06435e-12 (261)
156<TR><TH rowspan="4">sherman2 <TD rowspan="4"> 1080 <TD rowspan="4"> 23094 <TH>Compute Time <TD style="background-color:red">0.00631754<TD>0.015052<TD>0.0247514 <TD> -<TD>0.0214425<TD>0.0217988
157<TR><TH>Solve Time <TD style="background-color:red">0.000478424<TD>0.000337998<TD>0.0010291 <TD> -<TD>0.00243152<TD>0.00246152
158<TR><TH>Total Time <TD style="background-color:red">0.00679597<TD>0.01539<TD>0.0257805 <TD> -<TD>0.023874<TD>0.0242603
159<TR><TH>Error(Iter) <TD> 1.83099e-15 <TD> 8.19351e-15 <TD> 2.625e-14 <TD> 1.3678e+69 (1080) <TD> 4.1911e-12 (7) <TD> 5.0299e-13 (12)
160<TR><TH rowspan="4">bcsstk01_SPD <TD rowspan="4"> 48 <TD rowspan="4"> 400 <TH>Compute Time <TD>0.000169079<TD>0.00010789<TD>0.000572538<TD>1.425e-06<TD>9.1612e-05<TD>8.3985e-05<TD style="background-color:red">5.6489e-05<TD>7.0913e-05<TD>0.000468251<TD>5.7389e-05<TD>8.0212e-05<TD>5.8394e-05<TD>0.000463017<TD>1.333e-06
161<TR><TH>Solve Time <TD>1.2288e-05<TD>1.1124e-05<TD>0.000286387<TD>8.5896e-05<TD>1.6381e-05<TD>1.6984e-05<TD style="background-color:red">3.095e-06<TD>4.115e-06<TD>0.000325438<TD>3.504e-06<TD>7.369e-06<TD>3.454e-06<TD>0.000294095<TD>6.0516e-05
162<TR><TH>Total Time <TD>0.000181367<TD>0.000119014<TD>0.000858925<TD>8.7321e-05<TD>0.000107993<TD>0.000100969<TD style="background-color:red">5.9584e-05<TD>7.5028e-05<TD>0.000793689<TD>6.0893e-05<TD>8.7581e-05<TD>6.1848e-05<TD>0.000757112<TD>6.1849e-05
163<TR><TH>Error(Iter) <TD> 1.03474e-16 <TD> 2.23046e-16 <TD> 2.01273e-16 <TD> 4.87455e-07 (48) <TD> 1.03553e-16 (2) <TD> 3.55965e-16 (2) <TD> 2.48189e-16 <TD> 1.88808e-16 <TD> 1.97976e-16 <TD> 2.37248e-16 <TD> 1.82701e-16 <TD> 2.71474e-16 <TD> 2.11322e-16 <TD> 3.547e-09 (48)
164<TR><TH rowspan="4">sherman1 <TD rowspan="4"> 1000 <TD rowspan="4"> 3750 <TH>Compute Time <TD>0.00228805<TD>0.00209231<TD>0.00528268<TD>9.846e-06<TD>0.00163522<TD>0.00162155<TD>0.000789259<TD style="background-color:red">0.000804495<TD>0.00438269
165<TR><TH>Solve Time <TD>0.000213788<TD>9.7983e-05<TD>0.000938831<TD>0.00629835<TD>0.000361764<TD>0.00078794<TD>4.3989e-05<TD style="background-color:red">2.5331e-05<TD>0.000917166
166<TR><TH>Total Time <TD>0.00250184<TD>0.00219029<TD>0.00622151<TD>0.0063082<TD>0.00199698<TD>0.00240949<TD>0.000833248<TD style="background-color:red">0.000829826<TD>0.00529986
167<TR><TH>Error(Iter) <TD> 1.16839e-16 <TD> 2.25968e-16 <TD> 2.59116e-16 <TD> 3.76779e-11 (248) <TD> 4.13343e-11 (4) <TD> 2.22347e-14 (10) <TD> 2.05861e-16 <TD> 1.83555e-16 <TD> 1.02917e-15
168<TR><TH rowspan="4">young1c <TD rowspan="4"> 841 <TD rowspan="4"> 4089 <TH>Compute Time <TD>0.00235843<TD style="background-color:red">0.00217228<TD>0.00568075<TD>1.2735e-05<TD>0.00264866<TD>0.00258236
169<TR><TH>Solve Time <TD>0.000329599<TD style="background-color:red">0.000168634<TD>0.00080118<TD>0.0534738<TD>0.00187193<TD>0.00450211
170<TR><TH>Total Time <TD>0.00268803<TD style="background-color:red">0.00234091<TD>0.00648193<TD>0.0534865<TD>0.00452059<TD>0.00708447
171<TR><TH>Error(Iter) <TD> 1.27029e-16 <TD> 2.81321e-16 <TD> 5.0492e-15 <TD> 8.0507e-11 (706) <TD> 3.00447e-12 (8) <TD> 1.46532e-12 (16)
172<TR><TH rowspan="4">mhd1280b <TD rowspan="4"> 1280 <TD rowspan="4"> 22778 <TH>Compute Time <TD>0.00234898<TD>0.00207079<TD>0.00570918<TD>2.5976e-05<TD>0.00302563<TD>0.00298036<TD>0.00144525<TD style="background-color:red">0.000919922<TD>0.00426444
173<TR><TH>Solve Time <TD>0.00103392<TD>0.000211911<TD>0.00105<TD>0.0110432<TD>0.000628287<TD>0.00392089<TD>0.000138303<TD style="background-color:red">6.2446e-05<TD>0.00097564
174<TR><TH>Total Time <TD>0.0033829<TD>0.0022827<TD>0.00675918<TD>0.0110692<TD>0.00365392<TD>0.00690124<TD>0.00158355<TD style="background-color:red">0.000982368<TD>0.00524008
175<TR><TH>Error(Iter) <TD> 1.32953e-16 <TD> 3.08646e-16 <TD> 6.734e-16 <TD> 8.83132e-11 (40) <TD> 1.51153e-16 (1) <TD> 6.08556e-16 (8) <TD> 1.89264e-16 <TD> 1.97477e-16 <TD> 6.68126e-09
176<TR><TH rowspan="4">crashbasis <TD rowspan="4"> 160000 <TD rowspan="4"> 1750416 <TH>Compute Time <TD>3.2019<TD>5.7892<TD>15.7573<TD style="background-color:red">0.00383515<TD>3.1006<TD>3.09921
177<TR><TH>Solve Time <TD>0.261915<TD>0.106225<TD>0.402141<TD style="background-color:red">1.49089<TD>0.24888<TD>0.443673
178<TR><TH>Total Time <TD>3.46381<TD>5.89542<TD>16.1594<TD style="background-color:red">1.49473<TD>3.34948<TD>3.54288
179<TR><TH>Error(Iter) <TD> 1.76348e-16 <TD> 4.58395e-16 <TD> 1.67982e-14 <TD> 8.64144e-11 (61) <TD> 8.5996e-12 (2) <TD> 6.04042e-14 (5)
180
181</TABLE>
182*/
183}
Note: See TracBrowser for help on using the repository browser.