Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
PartialPivLU.h
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 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 template<typename MatrixType_> struct traits<PartialPivLU<MatrixType_> >
20  : traits<MatrixType_>
21 {
22  typedef MatrixXpr XprKind;
23  typedef SolverStorage StorageKind;
24  typedef int StorageIndex;
25  typedef traits<MatrixType_> BaseTraits;
26  enum {
27  Flags = BaseTraits::Flags & RowMajorBit,
28  CoeffReadCost = Dynamic
29  };
30 };
31 
32 template<typename T,typename Derived>
33 struct enable_if_ref;
34 // {
35 // typedef Derived type;
36 // };
37 
38 template<typename T,typename Derived>
39 struct enable_if_ref<Ref<T>,Derived> {
40  typedef Derived type;
41 };
42 
43 } // end namespace internal
44 
78 template<typename MatrixType_> class PartialPivLU
79  : public SolverBase<PartialPivLU<MatrixType_> >
80 {
81  public:
82 
83  typedef MatrixType_ MatrixType;
85  friend class SolverBase<PartialPivLU>;
86 
87  EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
88  enum {
89  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
90  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
91  };
94  typedef typename MatrixType::PlainObject PlainObject;
95 
102  PartialPivLU();
103 
110  explicit PartialPivLU(Index size);
111 
119  template<typename InputType>
120  explicit PartialPivLU(const EigenBase<InputType>& matrix);
121 
129  template<typename InputType>
130  explicit PartialPivLU(EigenBase<InputType>& matrix);
131 
132  template<typename InputType>
133  PartialPivLU& compute(const EigenBase<InputType>& matrix) {
134  m_lu = matrix.derived();
135  compute();
136  return *this;
137  }
138 
145  inline const MatrixType& matrixLU() const
146  {
147  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
148  return m_lu;
149  }
150 
153  inline const PermutationType& permutationP() const
154  {
155  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
156  return m_p;
157  }
158 
159  #ifdef EIGEN_PARSED_BY_DOXYGEN
177  template<typename Rhs>
178  inline const Solve<PartialPivLU, Rhs>
179  solve(const MatrixBase<Rhs>& b) const;
180  #endif
181 
185  inline RealScalar rcond() const
186  {
187  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
188  return internal::rcond_estimate_helper(m_l1_norm, *this);
189  }
190 
198  inline const Inverse<PartialPivLU> inverse() const
199  {
200  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
201  return Inverse<PartialPivLU>(*this);
202  }
203 
217  Scalar determinant() const;
218 
219  MatrixType reconstructedMatrix() const;
220 
221  EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
222  EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
223 
224  #ifndef EIGEN_PARSED_BY_DOXYGEN
225  template<typename RhsType, typename DstType>
226  EIGEN_DEVICE_FUNC
227  void _solve_impl(const RhsType &rhs, DstType &dst) const {
228  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
229  * So we proceed as follows:
230  * Step 1: compute c = Pb.
231  * Step 2: replace c by the solution x to Lx = c.
232  * Step 3: replace c by the solution x to Ux = c.
233  */
234 
235  // Step 1
236  dst = permutationP() * rhs;
237 
238  // Step 2
239  m_lu.template triangularView<UnitLower>().solveInPlace(dst);
240 
241  // Step 3
242  m_lu.template triangularView<Upper>().solveInPlace(dst);
243  }
244 
245  template<bool Conjugate, typename RhsType, typename DstType>
246  EIGEN_DEVICE_FUNC
247  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
248  /* The decomposition PA = LU can be rewritten as A^T = U^T L^T P.
249  * So we proceed as follows:
250  * Step 1: compute c as the solution to L^T c = b
251  * Step 2: replace c by the solution x to U^T x = c.
252  * Step 3: update c = P^-1 c.
253  */
254 
255  eigen_assert(rhs.rows() == m_lu.cols());
256 
257  // Step 1
258  dst = m_lu.template triangularView<Upper>().transpose()
259  .template conjugateIf<Conjugate>().solve(rhs);
260  // Step 2
261  m_lu.template triangularView<UnitLower>().transpose()
262  .template conjugateIf<Conjugate>().solveInPlace(dst);
263  // Step 3
264  dst = permutationP().transpose() * dst;
265  }
266  #endif
267 
268  protected:
269 
270  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
271 
272  void compute();
273 
274  MatrixType m_lu;
275  PermutationType m_p;
276  TranspositionType m_rowsTranspositions;
277  RealScalar m_l1_norm;
278  signed char m_det_p;
279  bool m_isInitialized;
280 };
281 
282 template<typename MatrixType>
284  : m_lu(),
285  m_p(),
286  m_rowsTranspositions(),
287  m_l1_norm(0),
288  m_det_p(0),
289  m_isInitialized(false)
290 {
291 }
292 
293 template<typename MatrixType>
295  : m_lu(size, size),
296  m_p(size),
297  m_rowsTranspositions(size),
298  m_l1_norm(0),
299  m_det_p(0),
300  m_isInitialized(false)
301 {
302 }
303 
304 template<typename MatrixType>
305 template<typename InputType>
307  : m_lu(matrix.rows(),matrix.cols()),
308  m_p(matrix.rows()),
309  m_rowsTranspositions(matrix.rows()),
310  m_l1_norm(0),
311  m_det_p(0),
312  m_isInitialized(false)
313 {
314  compute(matrix.derived());
315 }
316 
317 template<typename MatrixType>
318 template<typename InputType>
320  : m_lu(matrix.derived()),
321  m_p(matrix.rows()),
322  m_rowsTranspositions(matrix.rows()),
323  m_l1_norm(0),
324  m_det_p(0),
325  m_isInitialized(false)
326 {
327  compute();
328 }
329 
330 namespace internal {
331 
333 template<typename Scalar, int StorageOrder, typename PivIndex, int SizeAtCompileTime=Dynamic>
334 struct partial_lu_impl
335 {
336  static constexpr int UnBlockedBound = 16;
337  static constexpr bool UnBlockedAtCompileTime = SizeAtCompileTime!=Dynamic && SizeAtCompileTime<=UnBlockedBound;
338  static constexpr int ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic;
339  // Remaining rows and columns at compile-time:
340  static constexpr int RRows = SizeAtCompileTime==2 ? 1 : Dynamic;
341  static constexpr int RCols = SizeAtCompileTime==2 ? 1 : Dynamic;
343  typedef Ref<MatrixType> MatrixTypeRef;
345  typedef typename MatrixType::RealScalar RealScalar;
346 
357  static Index unblocked_lu(MatrixTypeRef& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
358  {
359  typedef scalar_score_coeff_op<Scalar> Scoring;
360  typedef typename Scoring::result_type Score;
361  const Index rows = lu.rows();
362  const Index cols = lu.cols();
363  const Index size = (std::min)(rows,cols);
364  // For small compile-time matrices it is worth processing the last row separately:
365  // speedup: +100% for 2x2, +10% for others.
366  const Index endk = UnBlockedAtCompileTime ? size-1 : size;
367  nb_transpositions = 0;
368  Index first_zero_pivot = -1;
369  for(Index k = 0; k < endk; ++k)
370  {
371  int rrows = internal::convert_index<int>(rows-k-1);
372  int rcols = internal::convert_index<int>(cols-k-1);
373 
374  Index row_of_biggest_in_col;
375  Score biggest_in_corner
376  = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
377  row_of_biggest_in_col += k;
378 
379  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
380 
381  if(!numext::is_exactly_zero(biggest_in_corner))
382  {
383  if(k != row_of_biggest_in_col)
384  {
385  lu.row(k).swap(lu.row(row_of_biggest_in_col));
386  ++nb_transpositions;
387  }
388 
389  lu.col(k).tail(fix<RRows>(rrows)) /= lu.coeff(k,k);
390  }
391  else if(first_zero_pivot==-1)
392  {
393  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
394  // and continue the factorization such we still have A = PLU
395  first_zero_pivot = k;
396  }
397 
398  if(k<rows-1)
399  lu.bottomRightCorner(fix<RRows>(rrows),fix<RCols>(rcols)).noalias() -= lu.col(k).tail(fix<RRows>(rrows)) * lu.row(k).tail(fix<RCols>(rcols));
400  }
401 
402  // special handling of the last entry
403  if(UnBlockedAtCompileTime)
404  {
405  Index k = endk;
406  row_transpositions[k] = PivIndex(k);
407  if (numext::is_exactly_zero(Scoring()(lu(k, k))) && first_zero_pivot == -1)
408  first_zero_pivot = k;
409  }
410 
411  return first_zero_pivot;
412  }
413 
429  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
430  {
431  MatrixTypeRef lu = MatrixType::Map(lu_data,rows, cols, OuterStride<>(luStride));
432 
433  const Index size = (std::min)(rows,cols);
434 
435  // if the matrix is too small, no blocking:
436  if(UnBlockedAtCompileTime || size<=UnBlockedBound)
437  {
438  return unblocked_lu(lu, row_transpositions, nb_transpositions);
439  }
440 
441  // automatically adjust the number of subdivisions to the size
442  // of the matrix so that there is enough sub blocks:
443  Index blockSize;
444  {
445  blockSize = size/8;
446  blockSize = (blockSize/16)*16;
447  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
448  }
449 
450  nb_transpositions = 0;
451  Index first_zero_pivot = -1;
452  for(Index k = 0; k < size; k+=blockSize)
453  {
454  Index bs = (std::min)(size-k,blockSize); // actual size of the block
455  Index trows = rows - k - bs; // trailing rows
456  Index tsize = size - k - bs; // trailing size
457 
458  // partition the matrix:
459  // A00 | A01 | A02
460  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
461  // A20 | A21 | A22
462  BlockType A_0 = lu.block(0,0,rows,k);
463  BlockType A_2 = lu.block(0,k+bs,rows,tsize);
464  BlockType A11 = lu.block(k,k,bs,bs);
465  BlockType A12 = lu.block(k,k+bs,bs,tsize);
466  BlockType A21 = lu.block(k+bs,k,trows,bs);
467  BlockType A22 = lu.block(k+bs,k+bs,trows,tsize);
468 
469  PivIndex nb_transpositions_in_panel;
470  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
471  // with a very small blocking size:
472  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
473  row_transpositions+k, nb_transpositions_in_panel, 16);
474  if(ret>=0 && first_zero_pivot==-1)
475  first_zero_pivot = k+ret;
476 
477  nb_transpositions += nb_transpositions_in_panel;
478  // update permutations and apply them to A_0
479  for(Index i=k; i<k+bs; ++i)
480  {
481  Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
482  A_0.row(i).swap(A_0.row(piv));
483  }
484 
485  if(trows)
486  {
487  // apply permutations to A_2
488  for(Index i=k;i<k+bs; ++i)
489  A_2.row(i).swap(A_2.row(row_transpositions[i]));
490 
491  // A12 = A11^-1 A12
492  A11.template triangularView<UnitLower>().solveInPlace(A12);
493 
494  A22.noalias() -= A21 * A12;
495  }
496  }
497  return first_zero_pivot;
498  }
499 };
500 
503 template<typename MatrixType, typename TranspositionType>
504 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
505 {
506  // Special-case of zero matrix.
507  if (lu.rows() == 0 || lu.cols() == 0) {
508  nb_transpositions = 0;
509  return;
510  }
511  eigen_assert(lu.cols() == row_transpositions.size());
512  eigen_assert(row_transpositions.size() < 2 || (&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
513 
514  partial_lu_impl
515  < typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor,
516  typename TranspositionType::StorageIndex,
517  internal::min_size_prefer_fixed(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)>
518  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
519 }
520 
521 } // end namespace internal
522 
523 template<typename MatrixType>
524 void PartialPivLU<MatrixType>::compute()
525 {
526  // the row permutation is stored as int indices, so just to be sure:
527  eigen_assert(m_lu.rows()<NumTraits<int>::highest());
528 
529  if(m_lu.cols()>0)
530  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
531  else
532  m_l1_norm = RealScalar(0);
533 
534  eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
535  const Index size = m_lu.rows();
536 
537  m_rowsTranspositions.resize(size);
538 
539  typename TranspositionType::StorageIndex nb_transpositions;
540  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
541  m_det_p = (nb_transpositions%2) ? -1 : 1;
542 
543  m_p = m_rowsTranspositions;
544 
545  m_isInitialized = true;
546 }
547 
548 template<typename MatrixType>
549 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
550 {
551  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
552  return Scalar(m_det_p) * m_lu.diagonal().prod();
553 }
554 
558 template<typename MatrixType>
560 {
561  eigen_assert(m_isInitialized && "LU is not initialized.");
562  // LU
563  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
564  * m_lu.template triangularView<Upper>();
565 
566  // P^{-1}(LU)
567  res = m_p.inverse() * res;
568 
569  return res;
570 }
571 
572 /***** Implementation details *****************************************************/
573 
574 namespace internal {
575 
576 /***** Implementation of inverse() *****************************************************/
577 template<typename DstXprType, typename MatrixType>
578 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
579 {
580  typedef PartialPivLU<MatrixType> LuType;
581  typedef Inverse<LuType> SrcXprType;
582  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
583  {
584  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
585  }
586 };
587 } // end namespace internal
588 
589 /******** MatrixBase methods *******/
590 
597 template<typename Derived>
598 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
600 {
601  return PartialPivLU<PlainObject>(eval());
602 }
603 
612 template<typename Derived>
615 {
616  return PartialPivLU<PlainObject>(eval());
617 }
618 
619 } // end namespace Eigen
620 
621 #endif // EIGEN_PARTIALLU_H
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:599
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:614
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
LU decomposition of a matrix with partial pivoting, and related features.
Definition: PartialPivLU.h:80
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
RealScalar rcond() const
Definition: PartialPivLU.h:185
Scalar determinant() const
Definition: PartialPivLU.h:549
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:283
const PermutationType & permutationP() const
Definition: PartialPivLU.h:153
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:198
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:559
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:145
InverseReturnType transpose() const
Definition: PermutationMatrix.h:193
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:642
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:285
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
@ ColMajor
Definition: Constants.h:321
@ RowMajor
Definition: Constants.h:323
const unsigned int RowMajorBit
Definition: Constants.h:68
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
EIGEN_CONSTEXPR Index size() const EIGEN_NOEXCEPT
Definition: EigenBase.h:69