PartialPivLU.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) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_PARTIALLU_H
12 #define EIGEN_PARTIALLU_H
13 
14 namespace Eigen {
15 
47 template<typename _MatrixType> class PartialPivLU
48 {
49  public:
50 
51  typedef _MatrixType MatrixType;
52  enum {
53  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55  Options = MatrixType::Options,
56  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
58  };
59  typedef typename MatrixType::Scalar Scalar;
62  typedef typename MatrixType::Index Index;
65 
66 
73  PartialPivLU();
74 
81  PartialPivLU(Index size);
82 
90  PartialPivLU(const MatrixType& matrix);
91 
92  PartialPivLU& compute(const MatrixType& matrix);
93 
100  inline const MatrixType& matrixLU() const
101  {
102  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
103  return m_lu;
104  }
105 
108  inline const PermutationType& permutationP() const
109  {
110  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
111  return m_p;
112  }
113 
131  template<typename Rhs>
133  solve(const MatrixBase<Rhs>& b) const
134  {
135  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
136  return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
137  }
138 
147  {
148  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
150  (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
151  }
152 
167 
168  MatrixType reconstructedMatrix() const;
169 
170  inline Index rows() const { return m_lu.rows(); }
171  inline Index cols() const { return m_lu.cols(); }
172 
173  protected:
174  MatrixType m_lu;
175  PermutationType m_p;
176  TranspositionType m_rowsTranspositions;
177  Index m_det_p;
179 };
180 
181 template<typename MatrixType>
183  : m_lu(),
184  m_p(),
186  m_det_p(0),
187  m_isInitialized(false)
188 {
189 }
190 
191 template<typename MatrixType>
193  : m_lu(size, size),
194  m_p(size),
195  m_rowsTranspositions(size),
196  m_det_p(0),
197  m_isInitialized(false)
198 {
199 }
200 
201 template<typename MatrixType>
203  : m_lu(matrix.rows(), matrix.rows()),
204  m_p(matrix.rows()),
205  m_rowsTranspositions(matrix.rows()),
206  m_det_p(0),
207  m_isInitialized(false)
208 {
209  compute(matrix);
210 }
211 
212 namespace internal {
213 
215 template<typename Scalar, int StorageOrder, typename PivIndex>
217 {
218  // FIXME add a stride to Map, so that the following mapping becomes easier,
219  // another option would be to create an expression being able to automatically
220  // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
221  // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
222  // and Block.
226  typedef typename MatrixType::RealScalar RealScalar;
227  typedef typename MatrixType::Index Index;
228 
239  static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
240  {
241  const Index rows = lu.rows();
242  const Index cols = lu.cols();
243  const Index size = (std::min)(rows,cols);
244  nb_transpositions = 0;
245  Index first_zero_pivot = -1;
246  for(Index k = 0; k < size; ++k)
247  {
248  Index rrows = rows-k-1;
249  Index rcols = cols-k-1;
250 
251  Index row_of_biggest_in_col;
252  RealScalar biggest_in_corner
253  = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
254  row_of_biggest_in_col += k;
255 
256  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
257 
258  if(biggest_in_corner != RealScalar(0))
259  {
260  if(k != row_of_biggest_in_col)
261  {
262  lu.row(k).swap(lu.row(row_of_biggest_in_col));
263  ++nb_transpositions;
264  }
265 
266  // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
267  // overflow but not the actual quotient?
268  lu.col(k).tail(rrows) /= lu.coeff(k,k);
269  }
270  else if(first_zero_pivot==-1)
271  {
272  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
273  // and continue the factorization such we still have A = PLU
274  first_zero_pivot = k;
275  }
276 
277  if(k<rows-1)
278  lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
279  }
280  return first_zero_pivot;
281  }
282 
298  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
299  {
300  MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
301  MatrixType lu(lu1,0,0,rows,cols);
302 
303  const Index size = (std::min)(rows,cols);
304 
305  // if the matrix is too small, no blocking:
306  if(size<=16)
307  {
308  return unblocked_lu(lu, row_transpositions, nb_transpositions);
309  }
310 
311  // automatically adjust the number of subdivisions to the size
312  // of the matrix so that there is enough sub blocks:
313  Index blockSize;
314  {
315  blockSize = size/8;
316  blockSize = (blockSize/16)*16;
317  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
318  }
319 
320  nb_transpositions = 0;
321  Index first_zero_pivot = -1;
322  for(Index k = 0; k < size; k+=blockSize)
323  {
324  Index bs = (std::min)(size-k,blockSize); // actual size of the block
325  Index trows = rows - k - bs; // trailing rows
326  Index tsize = size - k - bs; // trailing size
327 
328  // partition the matrix:
329  // A00 | A01 | A02
330  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
331  // A20 | A21 | A22
332  BlockType A_0(lu,0,0,rows,k);
333  BlockType A_2(lu,0,k+bs,rows,tsize);
334  BlockType A11(lu,k,k,bs,bs);
335  BlockType A12(lu,k,k+bs,bs,tsize);
336  BlockType A21(lu,k+bs,k,trows,bs);
337  BlockType A22(lu,k+bs,k+bs,trows,tsize);
338 
339  PivIndex nb_transpositions_in_panel;
340  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
341  // with a very small blocking size:
342  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
343  row_transpositions+k, nb_transpositions_in_panel, 16);
344  if(ret>=0 && first_zero_pivot==-1)
345  first_zero_pivot = k+ret;
346 
347  nb_transpositions += nb_transpositions_in_panel;
348  // update permutations and apply them to A_0
349  for(Index i=k; i<k+bs; ++i)
350  {
351  Index piv = (row_transpositions[i] += k);
352  A_0.row(i).swap(A_0.row(piv));
353  }
354 
355  if(trows)
356  {
357  // apply permutations to A_2
358  for(Index i=k;i<k+bs; ++i)
359  A_2.row(i).swap(A_2.row(row_transpositions[i]));
360 
361  // A12 = A11^-1 A12
362  A11.template triangularView<UnitLower>().solveInPlace(A12);
363 
364  A22.noalias() -= A21 * A12;
365  }
366  }
367  return first_zero_pivot;
368  }
369 };
370 
373 template<typename MatrixType, typename TranspositionType>
374 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions)
375 {
376  eigen_assert(lu.cols() == row_transpositions.size());
377  eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
378 
380  <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index>
381  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
382 }
383 
384 } // end namespace internal
385 
386 template<typename MatrixType>
388 {
389  // the row permutation is stored as int indices, so just to be sure:
390  eigen_assert(matrix.rows()<NumTraits<int>::highest());
391 
392  m_lu = matrix;
393 
394  eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
395  const Index size = matrix.rows();
396 
398 
399  typename TranspositionType::Index nb_transpositions;
400  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
401  m_det_p = (nb_transpositions%2) ? -1 : 1;
402 
404 
405  m_isInitialized = true;
406  return *this;
407 }
408 
409 template<typename MatrixType>
411 {
412  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
413  return Scalar(m_det_p) * m_lu.diagonal().prod();
414 }
415 
419 template<typename MatrixType>
421 {
422  eigen_assert(m_isInitialized && "LU is not initialized.");
423  // LU
424  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
425  * m_lu.template triangularView<Upper>();
426 
427  // P^{-1}(LU)
428  res = m_p.inverse() * res;
429 
430  return res;
431 }
432 
433 /***** Implementation of solve() *****************************************************/
434 
435 namespace internal {
436 
437 template<typename _MatrixType, typename Rhs>
438 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
439  : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
440 {
442 
443  template<typename Dest> void evalTo(Dest& dst) const
444  {
445  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
446  * So we proceed as follows:
447  * Step 1: compute c = Pb.
448  * Step 2: replace c by the solution x to Lx = c.
449  * Step 3: replace c by the solution x to Ux = c.
450  */
451 
452  eigen_assert(rhs().rows() == dec().matrixLU().rows());
453 
454  // Step 1
455  dst = dec().permutationP() * rhs();
456 
457  // Step 2
458  dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
459 
460  // Step 3
461  dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
462  }
463 };
464 
465 } // end namespace internal
466 
467 /******** MatrixBase methods *******/
468 
475 template<typename Derived>
478 {
479  return PartialPivLU<PlainObject>(eval());
480 }
481 
482 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
483 
491 template<typename Derived>
494 {
495  return PartialPivLU<PlainObject>(eval());
496 }
497 #endif
498 
499 } // end namespace Eigen
500 
501 #endif // EIGEN_PARTIALLU_H
internal::traits< Derived >::Scalar Scalar
MatrixType::Scalar Scalar
Definition: PartialPivLU.h:59
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:182
PermutationType m_p
Definition: PartialPivLU.h:175
Block< MatrixType, Dynamic, Dynamic > BlockType
Definition: PartialPivLU.h:225
Map< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > > MapLU
Definition: PartialPivLU.h:223
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
internal::traits< Derived >::Index Index
PartialPivLU & compute(const MatrixType &matrix)
Definition: PartialPivLU.h:387
LU decomposition of a matrix with partial pivoting, and related features.
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
void partial_lu_inplace(MatrixType &lu, TranspositionType &row_transpositions, typename TranspositionType::Index &nb_transpositions)
Definition: PartialPivLU.h:374
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
MatrixType::Index Index
Definition: PartialPivLU.h:62
const internal::solve_retval< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:133
Transpositions< RowsAtCompileTime, MaxRowsAtCompileTime > TranspositionType
Definition: PartialPivLU.h:64
MatrixType::RealScalar RealScalar
Definition: PartialPivLU.h:226
const unsigned int RowMajorBit
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:477
void resize(int newSize)
const PermutationType & permutationP() const
Definition: PartialPivLU.h:108
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
const internal::solve_retval< PartialPivLU, typename MatrixType::IdentityReturnType > inverse() const
Definition: PartialPivLU.h:146
Transpose< PermutationBase > inverse() const
void rhs(const real_t *x, real_t *f)
PermutationMatrix< RowsAtCompileTime, MaxRowsAtCompileTime > PermutationType
Definition: PartialPivLU.h:63
Expression of a fixed-size or dynamic-size block.
Definition: Core/Block.h:102
Index cols() const
Definition: PartialPivLU.h:171
static Index unblocked_lu(MatrixType &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions)
Definition: PartialPivLU.h:239
Block< MapLU, Dynamic, Dynamic > MatrixType
Definition: PartialPivLU.h:224
internal::traits< MatrixType >::Scalar determinant() const
Definition: PartialPivLU.h:410
TranspositionType m_rowsTranspositions
Definition: PartialPivLU.h:176
NumTraits< typename MatrixType::Scalar >::Real RealScalar
Definition: PartialPivLU.h:60
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:420
Index rows() const
Definition: PartialPivLU.h:170
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:100
#define eigen_assert(x)
internal::traits< MatrixType >::StorageKind StorageKind
Definition: PartialPivLU.h:61
static Index blocked_lu(Index rows, Index cols, Scalar *lu_data, Index luStride, PivIndex *row_transpositions, PivIndex &nb_transpositions, Index maxBlockSize=256)
Definition: PartialPivLU.h:298
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
_MatrixType MatrixType
Definition: PartialPivLU.h:51


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