1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008 Guillaume Saupin <guillaume.saupin@cea.fr>
|
---|
5 | //
|
---|
6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
9 |
|
---|
10 | #ifndef EIGEN_SKYLINEINPLACELU_H
|
---|
11 | #define EIGEN_SKYLINEINPLACELU_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \ingroup Skyline_Module
|
---|
16 | *
|
---|
17 | * \class SkylineInplaceLU
|
---|
18 | *
|
---|
19 | * \brief Inplace LU decomposition of a skyline matrix and associated features
|
---|
20 | *
|
---|
21 | * \param MatrixType the type of the matrix of which we are computing the LU factorization
|
---|
22 | *
|
---|
23 | */
|
---|
24 | template<typename MatrixType>
|
---|
25 | class SkylineInplaceLU {
|
---|
26 | protected:
|
---|
27 | typedef typename MatrixType::Scalar Scalar;
|
---|
28 | typedef typename MatrixType::Index Index;
|
---|
29 |
|
---|
30 | typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
|
---|
31 |
|
---|
32 | public:
|
---|
33 |
|
---|
34 | /** Creates a LU object and compute the respective factorization of \a matrix using
|
---|
35 | * flags \a flags. */
|
---|
36 | SkylineInplaceLU(MatrixType& matrix, int flags = 0)
|
---|
37 | : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0), m_lu(matrix) {
|
---|
38 | m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
|
---|
39 | m_lu.IsRowMajor ? computeRowMajor() : compute();
|
---|
40 | }
|
---|
41 |
|
---|
42 | /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
|
---|
43 | *
|
---|
44 | * Setting a value greater than zero speeds up computation, and yields to an imcomplete
|
---|
45 | * factorization with fewer non zero coefficients. Such approximate factors are especially
|
---|
46 | * useful to initialize an iterative solver.
|
---|
47 | *
|
---|
48 | * Note that the exact meaning of this parameter might depends on the actual
|
---|
49 | * backend. Moreover, not all backends support this feature.
|
---|
50 | *
|
---|
51 | * \sa precision() */
|
---|
52 | void setPrecision(RealScalar v) {
|
---|
53 | m_precision = v;
|
---|
54 | }
|
---|
55 |
|
---|
56 | /** \returns the current precision.
|
---|
57 | *
|
---|
58 | * \sa setPrecision() */
|
---|
59 | RealScalar precision() const {
|
---|
60 | return m_precision;
|
---|
61 | }
|
---|
62 |
|
---|
63 | /** Sets the flags. Possible values are:
|
---|
64 | * - CompleteFactorization
|
---|
65 | * - IncompleteFactorization
|
---|
66 | * - MemoryEfficient
|
---|
67 | * - one of the ordering methods
|
---|
68 | * - etc...
|
---|
69 | *
|
---|
70 | * \sa flags() */
|
---|
71 | void setFlags(int f) {
|
---|
72 | m_flags = f;
|
---|
73 | }
|
---|
74 |
|
---|
75 | /** \returns the current flags */
|
---|
76 | int flags() const {
|
---|
77 | return m_flags;
|
---|
78 | }
|
---|
79 |
|
---|
80 | void setOrderingMethod(int m) {
|
---|
81 | m_flags = m;
|
---|
82 | }
|
---|
83 |
|
---|
84 | int orderingMethod() const {
|
---|
85 | return m_flags;
|
---|
86 | }
|
---|
87 |
|
---|
88 | /** Computes/re-computes the LU factorization */
|
---|
89 | void compute();
|
---|
90 | void computeRowMajor();
|
---|
91 |
|
---|
92 | /** \returns the lower triangular matrix L */
|
---|
93 | //inline const MatrixType& matrixL() const { return m_matrixL; }
|
---|
94 |
|
---|
95 | /** \returns the upper triangular matrix U */
|
---|
96 | //inline const MatrixType& matrixU() const { return m_matrixU; }
|
---|
97 |
|
---|
98 | template<typename BDerived, typename XDerived>
|
---|
99 | bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
|
---|
100 | const int transposed = 0) const;
|
---|
101 |
|
---|
102 | /** \returns true if the factorization succeeded */
|
---|
103 | inline bool succeeded(void) const {
|
---|
104 | return m_succeeded;
|
---|
105 | }
|
---|
106 |
|
---|
107 | protected:
|
---|
108 | RealScalar m_precision;
|
---|
109 | int m_flags;
|
---|
110 | mutable int m_status;
|
---|
111 | bool m_succeeded;
|
---|
112 | MatrixType& m_lu;
|
---|
113 | };
|
---|
114 |
|
---|
115 | /** Computes / recomputes the in place LU decomposition of the SkylineInplaceLU.
|
---|
116 | * using the default algorithm.
|
---|
117 | */
|
---|
118 | template<typename MatrixType>
|
---|
119 | //template<typename _Scalar>
|
---|
120 | void SkylineInplaceLU<MatrixType>::compute() {
|
---|
121 | const size_t rows = m_lu.rows();
|
---|
122 | const size_t cols = m_lu.cols();
|
---|
123 |
|
---|
124 | eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
|
---|
125 | eigen_assert(!m_lu.IsRowMajor && "LU decomposition does not work with rowMajor Storage");
|
---|
126 |
|
---|
127 | for (Index row = 0; row < rows; row++) {
|
---|
128 | const double pivot = m_lu.coeffDiag(row);
|
---|
129 |
|
---|
130 | //Lower matrix Columns update
|
---|
131 | const Index& col = row;
|
---|
132 | for (typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
|
---|
133 | lIt.valueRef() /= pivot;
|
---|
134 | }
|
---|
135 |
|
---|
136 | //Upper matrix update -> contiguous memory access
|
---|
137 | typename MatrixType::InnerLowerIterator lIt(m_lu, col);
|
---|
138 | for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
|
---|
139 | typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
|
---|
140 | typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
|
---|
141 | const double coef = lIt.value();
|
---|
142 |
|
---|
143 | uItPivot += (rrow - row - 1);
|
---|
144 |
|
---|
145 | //update upper part -> contiguous memory access
|
---|
146 | for (++uItPivot; uIt && uItPivot;) {
|
---|
147 | uIt.valueRef() -= uItPivot.value() * coef;
|
---|
148 |
|
---|
149 | ++uIt;
|
---|
150 | ++uItPivot;
|
---|
151 | }
|
---|
152 | ++lIt;
|
---|
153 | }
|
---|
154 |
|
---|
155 | //Upper matrix update -> non contiguous memory access
|
---|
156 | typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
|
---|
157 | for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
|
---|
158 | typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
|
---|
159 | const double coef = lIt3.value();
|
---|
160 |
|
---|
161 | //update lower part -> non contiguous memory access
|
---|
162 | for (Index i = 0; i < rrow - row - 1; i++) {
|
---|
163 | m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
|
---|
164 | ++uItPivot;
|
---|
165 | }
|
---|
166 | ++lIt3;
|
---|
167 | }
|
---|
168 | //update diag -> contiguous
|
---|
169 | typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
|
---|
170 | for (Index rrow = row + 1; rrow < m_lu.rows(); rrow++) {
|
---|
171 |
|
---|
172 | typename MatrixType::InnerUpperIterator uItPivot(m_lu, row);
|
---|
173 | typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
|
---|
174 | const double coef = lIt2.value();
|
---|
175 |
|
---|
176 | uItPivot += (rrow - row - 1);
|
---|
177 | m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
|
---|
178 | ++lIt2;
|
---|
179 | }
|
---|
180 | }
|
---|
181 | }
|
---|
182 |
|
---|
183 | template<typename MatrixType>
|
---|
184 | void SkylineInplaceLU<MatrixType>::computeRowMajor() {
|
---|
185 | const size_t rows = m_lu.rows();
|
---|
186 | const size_t cols = m_lu.cols();
|
---|
187 |
|
---|
188 | eigen_assert(rows == cols && "We do not (yet) support rectangular LU.");
|
---|
189 | eigen_assert(m_lu.IsRowMajor && "You're trying to apply rowMajor decomposition on a ColMajor matrix !");
|
---|
190 |
|
---|
191 | for (Index row = 0; row < rows; row++) {
|
---|
192 | typename MatrixType::InnerLowerIterator llIt(m_lu, row);
|
---|
193 |
|
---|
194 |
|
---|
195 | for (Index col = llIt.col(); col < row; col++) {
|
---|
196 | if (m_lu.coeffExistLower(row, col)) {
|
---|
197 | const double diag = m_lu.coeffDiag(col);
|
---|
198 |
|
---|
199 | typename MatrixType::InnerLowerIterator lIt(m_lu, row);
|
---|
200 | typename MatrixType::InnerUpperIterator uIt(m_lu, col);
|
---|
201 |
|
---|
202 |
|
---|
203 | const Index offset = lIt.col() - uIt.row();
|
---|
204 |
|
---|
205 |
|
---|
206 | Index stop = offset > 0 ? col - lIt.col() : col - uIt.row();
|
---|
207 |
|
---|
208 | //#define VECTORIZE
|
---|
209 | #ifdef VECTORIZE
|
---|
210 | Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
|
---|
211 | Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
|
---|
212 |
|
---|
213 |
|
---|
214 | Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
|
---|
215 | #else
|
---|
216 | if (offset > 0) //Skip zero value of lIt
|
---|
217 | uIt += offset;
|
---|
218 | else //Skip zero values of uIt
|
---|
219 | lIt += -offset;
|
---|
220 | Scalar newCoeff = m_lu.coeffLower(row, col);
|
---|
221 |
|
---|
222 | for (Index k = 0; k < stop; ++k) {
|
---|
223 | const Scalar tmp = newCoeff;
|
---|
224 | newCoeff = tmp - lIt.value() * uIt.value();
|
---|
225 | ++lIt;
|
---|
226 | ++uIt;
|
---|
227 | }
|
---|
228 | #endif
|
---|
229 |
|
---|
230 | m_lu.coeffRefLower(row, col) = newCoeff / diag;
|
---|
231 | }
|
---|
232 | }
|
---|
233 |
|
---|
234 | //Upper matrix update
|
---|
235 | const Index col = row;
|
---|
236 | typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
|
---|
237 | for (Index rrow = uuIt.row(); rrow < col; rrow++) {
|
---|
238 |
|
---|
239 | typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
|
---|
240 | typename MatrixType::InnerUpperIterator uIt(m_lu, col);
|
---|
241 | const Index offset = lIt.col() - uIt.row();
|
---|
242 |
|
---|
243 | Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
|
---|
244 |
|
---|
245 | #ifdef VECTORIZE
|
---|
246 | Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
|
---|
247 | Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
|
---|
248 |
|
---|
249 | Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
|
---|
250 | #else
|
---|
251 | if (offset > 0) //Skip zero value of lIt
|
---|
252 | uIt += offset;
|
---|
253 | else //Skip zero values of uIt
|
---|
254 | lIt += -offset;
|
---|
255 | Scalar newCoeff = m_lu.coeffUpper(rrow, col);
|
---|
256 | for (Index k = 0; k < stop; ++k) {
|
---|
257 | const Scalar tmp = newCoeff;
|
---|
258 | newCoeff = tmp - lIt.value() * uIt.value();
|
---|
259 |
|
---|
260 | ++lIt;
|
---|
261 | ++uIt;
|
---|
262 | }
|
---|
263 | #endif
|
---|
264 | m_lu.coeffRefUpper(rrow, col) = newCoeff;
|
---|
265 | }
|
---|
266 |
|
---|
267 |
|
---|
268 | //Diag matrix update
|
---|
269 | typename MatrixType::InnerLowerIterator lIt(m_lu, row);
|
---|
270 | typename MatrixType::InnerUpperIterator uIt(m_lu, row);
|
---|
271 |
|
---|
272 | const Index offset = lIt.col() - uIt.row();
|
---|
273 |
|
---|
274 |
|
---|
275 | Index stop = offset > 0 ? lIt.size() : uIt.size();
|
---|
276 | #ifdef VECTORIZE
|
---|
277 | Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
|
---|
278 | Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
|
---|
279 | Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
|
---|
280 | #else
|
---|
281 | if (offset > 0) //Skip zero value of lIt
|
---|
282 | uIt += offset;
|
---|
283 | else //Skip zero values of uIt
|
---|
284 | lIt += -offset;
|
---|
285 | Scalar newCoeff = m_lu.coeffDiag(row);
|
---|
286 | for (Index k = 0; k < stop; ++k) {
|
---|
287 | const Scalar tmp = newCoeff;
|
---|
288 | newCoeff = tmp - lIt.value() * uIt.value();
|
---|
289 | ++lIt;
|
---|
290 | ++uIt;
|
---|
291 | }
|
---|
292 | #endif
|
---|
293 | m_lu.coeffRefDiag(row) = newCoeff;
|
---|
294 | }
|
---|
295 | }
|
---|
296 |
|
---|
297 | /** Computes *x = U^-1 L^-1 b
|
---|
298 | *
|
---|
299 | * If \a transpose is set to SvTranspose or SvAdjoint, the solution
|
---|
300 | * of the transposed/adjoint system is computed instead.
|
---|
301 | *
|
---|
302 | * Not all backends implement the solution of the transposed or
|
---|
303 | * adjoint system.
|
---|
304 | */
|
---|
305 | template<typename MatrixType>
|
---|
306 | template<typename BDerived, typename XDerived>
|
---|
307 | bool SkylineInplaceLU<MatrixType>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed) const {
|
---|
308 | const size_t rows = m_lu.rows();
|
---|
309 | const size_t cols = m_lu.cols();
|
---|
310 |
|
---|
311 |
|
---|
312 | for (Index row = 0; row < rows; row++) {
|
---|
313 | x->coeffRef(row) = b.coeff(row);
|
---|
314 | Scalar newVal = x->coeff(row);
|
---|
315 | typename MatrixType::InnerLowerIterator lIt(m_lu, row);
|
---|
316 |
|
---|
317 | Index col = lIt.col();
|
---|
318 | while (lIt.col() < row) {
|
---|
319 |
|
---|
320 | newVal -= x->coeff(col++) * lIt.value();
|
---|
321 | ++lIt;
|
---|
322 | }
|
---|
323 |
|
---|
324 | x->coeffRef(row) = newVal;
|
---|
325 | }
|
---|
326 |
|
---|
327 |
|
---|
328 | for (Index col = rows - 1; col > 0; col--) {
|
---|
329 | x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
|
---|
330 |
|
---|
331 | const Scalar x_col = x->coeff(col);
|
---|
332 |
|
---|
333 | typename MatrixType::InnerUpperIterator uIt(m_lu, col);
|
---|
334 | uIt += uIt.size()-1;
|
---|
335 |
|
---|
336 |
|
---|
337 | while (uIt) {
|
---|
338 | x->coeffRef(uIt.row()) -= x_col * uIt.value();
|
---|
339 | //TODO : introduce --operator
|
---|
340 | uIt += -1;
|
---|
341 | }
|
---|
342 |
|
---|
343 |
|
---|
344 | }
|
---|
345 | x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
|
---|
346 |
|
---|
347 | return true;
|
---|
348 | }
|
---|
349 |
|
---|
350 | } // end namespace Eigen
|
---|
351 |
|
---|
352 | #endif // EIGEN_SKYLINELU_H
|
---|