SkylineInplaceLU.h
Go to the documentation of this file.
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 
24 template<typename MatrixType>
26 protected:
27  typedef typename MatrixType::Scalar Scalar;
28  typedef typename MatrixType::Index Index;
29 
31 
32 public:
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>
100  const int transposed = 0) const;
101 
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;
112  MatrixType& m_lu;
113 };
114 
118 template<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 
183 template<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 
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
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
NumTraits< typename MatrixType::Scalar >::Real RealScalar
SkylineInplaceLU(MatrixType &matrix, int flags=0)
Inplace LU decomposition of a skyline matrix and associated features.
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
MatrixType::Index Index
bool succeeded(void) const
bool solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
CoeffReturnType value() const
Definition: DenseBase.h:424
RealScalar precision() const
MatrixType::Scalar Scalar
#define v
RowXpr row(Index i)
Definition: BlockMethods.h:725
ColXpr col(Index i)
Definition: BlockMethods.h:708
#define eigen_assert(x)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void setPrecision(RealScalar v)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:06