Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
SkylineInplaceLU.h
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
13namespace Eigen {
14
24template<typename MatrixType>
26protected:
27 typedef typename MatrixType::Scalar Scalar;
28 typedef typename MatrixType::Index Index;
29
30 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
31
32public:
33
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
52 void setPrecision(RealScalar v) {
53 m_precision = v;
54 }
55
59 RealScalar precision() const {
60 return m_precision;
61 }
62
71 void setFlags(int f) {
72 m_flags = f;
73 }
74
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
89 void compute();
90 void computeRowMajor();
91
93 //inline const MatrixType& matrixL() const { return m_matrixL; }
94
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
103 inline bool succeeded(void) const {
104 return m_succeeded;
105 }
106
107protected:
108 RealScalar m_precision;
109 int m_flags;
110 mutable int m_status;
111 bool m_succeeded;
112 MatrixType& m_lu;
113};
114
118template<typename MatrixType>
119//template<typename _Scalar>
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
183template<typename MatrixType>
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
305template<typename MatrixType>
306template<typename BDerived, typename XDerived>
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_SKYLINEINPLACELU_H
Inplace LU decomposition of a skyline matrix and associated features.
Definition: SkylineInplaceLU.h:25
RealScalar precision() const
Definition: SkylineInplaceLU.h:59
void setPrecision(RealScalar v)
Definition: SkylineInplaceLU.h:52
bool solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
Definition: SkylineInplaceLU.h:307
void compute()
Definition: SkylineInplaceLU.h:120
int flags() const
Definition: SkylineInplaceLU.h:76
bool succeeded(void) const
Definition: SkylineInplaceLU.h:103
SkylineInplaceLU(MatrixType &matrix, int flags=0)
Definition: SkylineInplaceLU.h:36
void setFlags(int f)
Definition: SkylineInplaceLU.h:71
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index