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)
m_flags(flags)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() */
setPrecision(RealScalar v)52 void setPrecision(RealScalar v) {
53 m_precision = v;
54 }
55
56 /** \returns the current precision.
57 *
58 * \sa setPrecision() */
precision()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() */
setFlags(int f)71 void setFlags(int f) {
72 m_flags = f;
73 }
74
75 /** \returns the current flags */
flags()76 int flags() const {
77 return m_flags;
78 }
79
setOrderingMethod(int m)80 void setOrderingMethod(int m) {
81 m_flags = m;
82 }
83
orderingMethod()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 */
succeeded(void)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>
compute()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>
computeRowMajor()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>
solve(const MatrixBase<BDerived> & b,MatrixBase<XDerived> * x,const int transposed)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
353