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 
53  m_precision = v;
54  }
55 
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:
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
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
col
m col(1)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::SkylineInplaceLU::m_status
int m_status
Definition: SkylineInplaceLU.h:110
b
Scalar * b
Definition: benchVecAdd.cpp:17
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
gtsam::diag
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
x
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
Definition: gnuplot_common_settings.hh:12
Eigen::SkylineInplaceLU::Scalar
MatrixType::Scalar Scalar
Definition: SkylineInplaceLU.h:27
Eigen::SkylineInplaceLU::setOrderingMethod
void setOrderingMethod(int m)
Definition: SkylineInplaceLU.h:80
Eigen::SkylineInplaceLU::RealScalar
NumTraits< typename MatrixType::Scalar >::Real RealScalar
Definition: SkylineInplaceLU.h:30
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
Eigen::SkylineInplaceLU::m_lu
MatrixType & m_lu
Definition: SkylineInplaceLU.h:112
Eigen::SkylineInplaceLU::computeRowMajor
void computeRowMajor()
Definition: SkylineInplaceLU.h:184
Eigen::SkylineInplaceLU::setPrecision
void setPrecision(RealScalar v)
Definition: SkylineInplaceLU.h:52
Eigen::SkylineInplaceLU::precision
RealScalar precision() const
Definition: SkylineInplaceLU.h:59
Eigen::SkylineInplaceLU::orderingMethod
int orderingMethod() const
Definition: SkylineInplaceLU.h:84
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::SkylineInplaceLU::solve
bool solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
Definition: SkylineInplaceLU.h:307
Eigen::SkylineInplaceLU::flags
int flags() const
Definition: SkylineInplaceLU.h:76
Eigen::SkylineInplaceLU::m_flags
int m_flags
Definition: SkylineInplaceLU.h:109
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
anyset::size
size_t size() const
Definition: pytypes.h:2220
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
offset
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
Definition: gnuplot_common_settings.hh:64
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
Eigen::SkylineInplaceLU::compute
void compute()
Definition: SkylineInplaceLU.h:120
Eigen::SkylineInplaceLU::setFlags
void setFlags(int f)
Definition: SkylineInplaceLU.h:71
Eigen::SkylineInplaceLU::m_succeeded
bool m_succeeded
Definition: SkylineInplaceLU.h:111
row
m row(1)
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
Eigen::SkylineInplaceLU::succeeded
bool succeeded(void) const
Definition: SkylineInplaceLU.h:103
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::SkylineInplaceLU
Inplace LU decomposition of a skyline matrix and associated features.
Definition: SkylineInplaceLU.h:25
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::SkylineInplaceLU::Index
MatrixType::Index Index
Definition: SkylineInplaceLU.h:28
Eigen::SkylineInplaceLU::SkylineInplaceLU
SkylineInplaceLU(MatrixType &matrix, int flags=0)
Definition: SkylineInplaceLU.h:36
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::SkylineInplaceLU::m_precision
RealScalar m_precision
Definition: SkylineInplaceLU.h:108
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:28