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 
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;
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>
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
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
Scalar * b
Definition: benchVecAdd.cpp:17
size_t size() const
Definition: pytypes.h:2022
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
NumTraits< typename MatrixType::Scalar >::Real RealScalar
SkylineInplaceLU(MatrixType &matrix, int flags=0)
Inplace LU decomposition of a skyline matrix and associated features.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate offset
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DEVICE_FUNC CoeffReturnType value() const
Definition: DenseBase.h:526
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
MatrixType::Index Index
bool succeeded(void) const
m row(1)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
Array< int, Dynamic, 1 > v
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
MatrixType::Scalar Scalar
m col(1)
bool solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void setPrecision(RealScalar v)
RealScalar precision() const


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:50