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 
16 namespace internal {
17 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
18  : traits<_MatrixType>
19 {
20  typedef MatrixXpr XprKind;
22  typedef int StorageIndex;
24  enum {
25  Flags = BaseTraits::Flags & RowMajorBit,
26  CoeffReadCost = Dynamic
27  };
28 };
29 
30 template<typename T,typename Derived>
32 // {
33 // typedef Derived type;
34 // };
35 
36 template<typename T,typename Derived>
37 struct enable_if_ref<Ref<T>,Derived> {
38  typedef Derived type;
39 };
40 
41 } // end namespace internal
42 
76 template<typename _MatrixType> class PartialPivLU
77  : public SolverBase<PartialPivLU<_MatrixType> >
78 {
79  public:
80 
81  typedef _MatrixType MatrixType;
83  friend class SolverBase<PartialPivLU>;
84 
86  enum {
87  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
88  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
89  };
93 
100  PartialPivLU();
101 
108  explicit PartialPivLU(Index size);
109 
117  template<typename InputType>
118  explicit PartialPivLU(const EigenBase<InputType>& matrix);
119 
127  template<typename InputType>
129 
130  template<typename InputType>
132  m_lu = matrix.derived();
133  compute();
134  return *this;
135  }
136 
143  inline const MatrixType& matrixLU() const
144  {
145  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
146  return m_lu;
147  }
148 
151  inline const PermutationType& permutationP() const
152  {
153  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
154  return m_p;
155  }
156 
157  #ifdef EIGEN_PARSED_BY_DOXYGEN
158 
175  template<typename Rhs>
176  inline const Solve<PartialPivLU, Rhs>
177  solve(const MatrixBase<Rhs>& b) const;
178  #endif
179 
183  inline RealScalar rcond() const
184  {
185  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
186  return internal::rcond_estimate_helper(m_l1_norm, *this);
187  }
188 
196  inline const Inverse<PartialPivLU> inverse() const
197  {
198  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
199  return Inverse<PartialPivLU>(*this);
200  }
201 
215  Scalar determinant() const;
216 
217  MatrixType reconstructedMatrix() const;
218 
219  EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
220  EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
221 
222  #ifndef EIGEN_PARSED_BY_DOXYGEN
223  template<typename RhsType, typename DstType>
225  void _solve_impl(const RhsType &rhs, DstType &dst) const {
226  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
227  * So we proceed as follows:
228  * Step 1: compute c = Pb.
229  * Step 2: replace c by the solution x to Lx = c.
230  * Step 3: replace c by the solution x to Ux = c.
231  */
232 
233  // Step 1
234  dst = permutationP() * rhs;
235 
236  // Step 2
237  m_lu.template triangularView<UnitLower>().solveInPlace(dst);
238 
239  // Step 3
240  m_lu.template triangularView<Upper>().solveInPlace(dst);
241  }
242 
243  template<bool Conjugate, typename RhsType, typename DstType>
245  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
246  /* The decomposition PA = LU can be rewritten as A^T = U^T L^T P.
247  * So we proceed as follows:
248  * Step 1: compute c as the solution to L^T c = b
249  * Step 2: replace c by the solution x to U^T x = c.
250  * Step 3: update c = P^-1 c.
251  */
252 
253  eigen_assert(rhs.rows() == m_lu.cols());
254 
255  // Step 1
256  dst = m_lu.template triangularView<Upper>().transpose()
257  .template conjugateIf<Conjugate>().solve(rhs);
258  // Step 2
259  m_lu.template triangularView<UnitLower>().transpose()
260  .template conjugateIf<Conjugate>().solveInPlace(dst);
261  // Step 3
262  dst = permutationP().transpose() * dst;
263  }
264  #endif
265 
266  protected:
267 
269  {
271  }
272 
273  void compute();
274 
275  MatrixType m_lu;
276  PermutationType m_p;
277  TranspositionType m_rowsTranspositions;
279  signed char m_det_p;
281 };
282 
283 template<typename MatrixType>
285  : m_lu(),
286  m_p(),
287  m_rowsTranspositions(),
288  m_l1_norm(0),
289  m_det_p(0),
290  m_isInitialized(false)
291 {
292 }
293 
294 template<typename MatrixType>
296  : m_lu(size, size),
297  m_p(size),
298  m_rowsTranspositions(size),
299  m_l1_norm(0),
300  m_det_p(0),
301  m_isInitialized(false)
302 {
303 }
304 
305 template<typename MatrixType>
306 template<typename InputType>
308  : m_lu(matrix.rows(),matrix.cols()),
309  m_p(matrix.rows()),
310  m_rowsTranspositions(matrix.rows()),
311  m_l1_norm(0),
312  m_det_p(0),
313  m_isInitialized(false)
314 {
315  compute(matrix.derived());
316 }
317 
318 template<typename MatrixType>
319 template<typename InputType>
321  : m_lu(matrix.derived()),
322  m_p(matrix.rows()),
323  m_rowsTranspositions(matrix.rows()),
324  m_l1_norm(0),
325  m_det_p(0),
326  m_isInitialized(false)
327 {
328  compute();
329 }
330 
331 namespace internal {
332 
334 template<typename Scalar, int StorageOrder, typename PivIndex, int SizeAtCompileTime=Dynamic>
336 {
337  static const int UnBlockedBound = 16;
338  static const bool UnBlockedAtCompileTime = SizeAtCompileTime!=Dynamic && SizeAtCompileTime<=UnBlockedBound;
339  static const int ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic;
340  // Remaining rows and columns at compile-time:
341  static const int RRows = SizeAtCompileTime==2 ? 1 : Dynamic;
342  static const int RCols = SizeAtCompileTime==2 ? 1 : Dynamic;
347 
358  static Index unblocked_lu(MatrixTypeRef& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
359  {
360  typedef scalar_score_coeff_op<Scalar> Scoring;
361  typedef typename Scoring::result_type Score;
362  const Index rows = lu.rows();
363  const Index cols = lu.cols();
364  const Index size = (std::min)(rows,cols);
365  // For small compile-time matrices it is worth processing the last row separately:
366  // speedup: +100% for 2x2, +10% for others.
367  const Index endk = UnBlockedAtCompileTime ? size-1 : size;
368  nb_transpositions = 0;
369  Index first_zero_pivot = -1;
370  for(Index k = 0; k < endk; ++k)
371  {
372  int rrows = internal::convert_index<int>(rows-k-1);
373  int rcols = internal::convert_index<int>(cols-k-1);
374 
375  Index row_of_biggest_in_col;
376  Score biggest_in_corner
377  = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
378  row_of_biggest_in_col += k;
379 
380  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
381 
382  if(biggest_in_corner != Score(0))
383  {
384  if(k != row_of_biggest_in_col)
385  {
386  lu.row(k).swap(lu.row(row_of_biggest_in_col));
388  }
389 
390  lu.col(k).tail(fix<RRows>(rrows)) /= lu.coeff(k,k);
391  }
392  else if(first_zero_pivot==-1)
393  {
394  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
395  // and continue the factorization such we still have A = PLU
396  first_zero_pivot = k;
397  }
398 
399  if(k<rows-1)
400  lu.bottomRightCorner(fix<RRows>(rrows),fix<RCols>(rcols)).noalias() -= lu.col(k).tail(fix<RRows>(rrows)) * lu.row(k).tail(fix<RCols>(rcols));
401  }
402 
403  // special handling of the last entry
404  if(UnBlockedAtCompileTime)
405  {
406  Index k = endk;
407  row_transpositions[k] = PivIndex(k);
408  if (Scoring()(lu(k, k)) == Score(0) && first_zero_pivot == -1)
409  first_zero_pivot = k;
410  }
411 
412  return first_zero_pivot;
413  }
414 
430  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
431  {
432  MatrixTypeRef lu = MatrixType::Map(lu_data,rows, cols, OuterStride<>(luStride));
433 
434  const Index size = (std::min)(rows,cols);
435 
436  // if the matrix is too small, no blocking:
437  if(UnBlockedAtCompileTime || size<=UnBlockedBound)
438  {
439  return unblocked_lu(lu, row_transpositions, nb_transpositions);
440  }
441 
442  // automatically adjust the number of subdivisions to the size
443  // of the matrix so that there is enough sub blocks:
444  Index blockSize;
445  {
446  blockSize = size/8;
447  blockSize = (blockSize/16)*16;
448  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
449  }
450 
451  nb_transpositions = 0;
452  Index first_zero_pivot = -1;
453  for(Index k = 0; k < size; k+=blockSize)
454  {
455  Index bs = (std::min)(size-k,blockSize); // actual size of the block
456  Index trows = rows - k - bs; // trailing rows
457  Index tsize = size - k - bs; // trailing size
458 
459  // partition the matrix:
460  // A00 | A01 | A02
461  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
462  // A20 | A21 | A22
463  BlockType A_0 = lu.block(0,0,rows,k);
464  BlockType A_2 = lu.block(0,k+bs,rows,tsize);
465  BlockType A11 = lu.block(k,k,bs,bs);
466  BlockType A12 = lu.block(k,k+bs,bs,tsize);
467  BlockType A21 = lu.block(k+bs,k,trows,bs);
468  BlockType A22 = lu.block(k+bs,k+bs,trows,tsize);
469 
470  PivIndex nb_transpositions_in_panel;
471  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
472  // with a very small blocking size:
473  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
474  row_transpositions+k, nb_transpositions_in_panel, 16);
475  if(ret>=0 && first_zero_pivot==-1)
476  first_zero_pivot = k+ret;
477 
478  nb_transpositions += nb_transpositions_in_panel;
479  // update permutations and apply them to A_0
480  for(Index i=k; i<k+bs; ++i)
481  {
482  Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
483  A_0.row(i).swap(A_0.row(piv));
484  }
485 
486  if(trows)
487  {
488  // apply permutations to A_2
489  for(Index i=k;i<k+bs; ++i)
490  A_2.row(i).swap(A_2.row(row_transpositions[i]));
491 
492  // A12 = A11^-1 A12
493  A11.template triangularView<UnitLower>().solveInPlace(A12);
494 
495  A22.noalias() -= A21 * A12;
496  }
497  }
498  return first_zero_pivot;
499  }
500 };
501 
504 template<typename MatrixType, typename TranspositionType>
505 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
506 {
507  // Special-case of zero matrix.
508  if (lu.rows() == 0 || lu.cols() == 0) {
509  nb_transpositions = 0;
510  return;
511  }
512  eigen_assert(lu.cols() == row_transpositions.size());
513  eigen_assert(row_transpositions.size() < 2 || (&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
514 
516  < typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor,
517  typename TranspositionType::StorageIndex,
518  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime)>
519  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
520 }
521 
522 } // end namespace internal
523 
524 template<typename MatrixType>
526 {
527  check_template_parameters();
528 
529  // the row permutation is stored as int indices, so just to be sure:
530  eigen_assert(m_lu.rows()<NumTraits<int>::highest());
531 
532  if(m_lu.cols()>0)
533  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
534  else
535  m_l1_norm = RealScalar(0);
536 
537  eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
538  const Index size = m_lu.rows();
539 
540  m_rowsTranspositions.resize(size);
541 
543  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
544  m_det_p = (nb_transpositions%2) ? -1 : 1;
545 
546  m_p = m_rowsTranspositions;
547 
548  m_isInitialized = true;
549 }
550 
551 template<typename MatrixType>
553 {
554  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
555  return Scalar(m_det_p) * m_lu.diagonal().prod();
556 }
557 
561 template<typename MatrixType>
563 {
564  eigen_assert(m_isInitialized && "LU is not initialized.");
565  // LU
566  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
567  * m_lu.template triangularView<Upper>();
568 
569  // P^{-1}(LU)
570  res = m_p.inverse() * res;
571 
572  return res;
573 }
574 
575 /***** Implementation details *****************************************************/
576 
577 namespace internal {
578 
579 /***** Implementation of inverse() *****************************************************/
580 template<typename DstXprType, typename MatrixType>
581 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
582 {
585  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
586  {
587  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
588  }
589 };
590 } // end namespace internal
591 
592 /******** MatrixBase methods *******/
593 
600 template<typename Derived>
603 {
605 }
606 
615 template<typename Derived>
618 {
620 }
621 
622 } // end namespace Eigen
623 
624 #endif // EIGEN_PARTIALLU_H
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: Inverse.h:57
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: PartialPivLU.h:220
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:1264
#define max(a, b)
Definition: datatypes.h:20
MatrixType::RealScalar RealScalar
Definition: PartialPivLU.h:346
MatrixType::PlainObject PlainObject
Definition: PartialPivLU.h:92
PartialPivLU & compute(const EigenBase< InputType > &matrix)
Definition: PartialPivLU.h:131
Scalar * b
Definition: benchVecAdd.cpp:17
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:284
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:143
internal::traits< PartialPivLU< Matrix< Scalar, Dynamic, Dynamic > > >::Scalar Scalar
Definition: SolverBase.h:73
PermutationType m_p
Definition: PartialPivLU.h:276
#define min(a, b)
Definition: datatypes.h:19
int nb_transpositions
Definition: lapack/lu.cpp:30
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op< typename DstXprType::Scalar, typename LuType::Scalar > &)
Definition: PartialPivLU.h:585
void determinant(const MatrixType &m)
Definition: determinant.cpp:14
LU decomposition of a matrix with partial pivoting, and related features.
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Transpositions< RowsAtCompileTime, MaxRowsAtCompileTime > TranspositionType
Definition: PartialPivLU.h:91
signed char m_det_p
Definition: PartialPivLU.h:279
Decomposition::RealScalar rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Decomposition &dec)
Reciprocal condition number estimator.
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
const unsigned int RowMajorBit
Definition: Constants.h:66
#define EIGEN_SIZE_MIN_PREFER_FIXED(a, b)
Definition: Macros.h:1302
Matrix< Scalar, ActualSizeAtCompileTime, ActualSizeAtCompileTime, StorageOrder > MatrixType
Definition: PartialPivLU.h:343
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: PartialPivLU.h:219
Ref< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > > BlockType
Definition: PartialPivLU.h:345
Expression of the inverse of another expression.
Definition: Inverse.h:43
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEVICE_FUNC void _solve_impl(const RhsType &rhs, DstType &dst) const
Definition: PartialPivLU.h:225
void partial_lu_inplace(MatrixType &lu, TranspositionType &row_transpositions, typename TranspositionType::StorageIndex &nb_transpositions)
Definition: PartialPivLU.h:505
SolverBase< PartialPivLU > Base
Definition: PartialPivLU.h:82
#define EIGEN_NOEXCEPT
Definition: Macros.h:1418
static void check_template_parameters()
Definition: PartialPivLU.h:268
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
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:196
#define EIGEN_CONSTEXPR
Definition: Macros.h:787
cout<< "Here is the matrix m:"<< endl<< m<< endl;Eigen::FullPivLU< Matrix5x3 > lu(m)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
DenseIndex ret
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
PermutationMatrix< RowsAtCompileTime, MaxRowsAtCompileTime > PermutationType
Definition: PartialPivLU.h:90
const PermutationType & permutationP() const
Definition: PartialPivLU.h:151
TranspositionType m_rowsTranspositions
Definition: PartialPivLU.h:277
internal::nested_eval< T, 1 >::type eval(const T &xpr)
const int Dynamic
Definition: Constants.h:22
Pseudo expression representing a solving operation.
Definition: Solve.h:62
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:430
static Index unblocked_lu(MatrixTypeRef &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions)
Definition: PartialPivLU.h:358
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
The matrix class, also used for vectors and row-vectors.
Convenience specialization of Stride to specify only an outer stride See class Map for some examples...
Definition: Stride.h:106
RealScalar m_l1_norm
Definition: PartialPivLU.h:278
RealScalar rcond() const
Definition: PartialPivLU.h:183
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:68
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
_MatrixType MatrixType
Definition: PartialPivLU.h:81
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: Inverse.h:58
EIGEN_DEVICE_FUNC const XprTypeNestedCleaned & nestedExpression() const
Definition: Inverse.h:60
EIGEN_DEVICE_FUNC void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
Definition: PartialPivLU.h:245


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