Eigen  3.3.7
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 namespace Eigen {
15 
16 namespace internal {
17 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
18  : traits<_MatrixType>
19 {
20  typedef MatrixXpr XprKind;
21  typedef SolverStorage StorageKind;
22  typedef traits<_MatrixType> BaseTraits;
23  enum {
24  Flags = BaseTraits::Flags & RowMajorBit,
25  CoeffReadCost = Dynamic
26  };
27 };
28 
29 template<typename T,typename Derived>
30 struct enable_if_ref;
31 // {
32 // typedef Derived type;
33 // };
34 
35 template<typename T,typename Derived>
36 struct enable_if_ref<Ref<T>,Derived> {
37  typedef Derived type;
38 };
39 
40 } // end namespace internal
41 
75 template<typename _MatrixType> class PartialPivLU
76  : public SolverBase<PartialPivLU<_MatrixType> >
77 {
78  public:
79 
80  typedef _MatrixType MatrixType;
81  typedef SolverBase<PartialPivLU> Base;
82  EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
83  // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
84  enum {
85  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
86  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
87  };
88  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
89  typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
90  typedef typename MatrixType::PlainObject PlainObject;
91 
98  PartialPivLU();
99 
106  explicit PartialPivLU(Index size);
107 
115  template<typename InputType>
116  explicit PartialPivLU(const EigenBase<InputType>& matrix);
117 
125  template<typename InputType>
126  explicit PartialPivLU(EigenBase<InputType>& matrix);
127 
128  template<typename InputType>
129  PartialPivLU& compute(const EigenBase<InputType>& matrix) {
130  m_lu = matrix.derived();
131  compute();
132  return *this;
133  }
134 
141  inline const MatrixType& matrixLU() const
142  {
143  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
144  return m_lu;
145  }
146 
149  inline const PermutationType& permutationP() const
150  {
151  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
152  return m_p;
153  }
154 
172  // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
173  template<typename Rhs>
174  inline const Solve<PartialPivLU, Rhs>
175  solve(const MatrixBase<Rhs>& b) const
176  {
177  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
178  return Solve<PartialPivLU, Rhs>(*this, b.derived());
179  }
180 
184  inline RealScalar rcond() const
185  {
186  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
187  return internal::rcond_estimate_helper(m_l1_norm, *this);
188  }
189 
197  inline const Inverse<PartialPivLU> inverse() const
198  {
199  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
200  return Inverse<PartialPivLU>(*this);
201  }
202 
216  Scalar determinant() const;
217 
218  MatrixType reconstructedMatrix() const;
219 
220  inline Index rows() const { return m_lu.rows(); }
221  inline Index cols() const { return m_lu.cols(); }
222 
223  #ifndef EIGEN_PARSED_BY_DOXYGEN
224  template<typename RhsType, typename DstType>
225  EIGEN_DEVICE_FUNC
226  void _solve_impl(const RhsType &rhs, DstType &dst) const {
227  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
228  * So we proceed as follows:
229  * Step 1: compute c = Pb.
230  * Step 2: replace c by the solution x to Lx = c.
231  * Step 3: replace c by the solution x to Ux = c.
232  */
233 
234  eigen_assert(rhs.rows() == m_lu.rows());
235 
236  // Step 1
237  dst = permutationP() * rhs;
238 
239  // Step 2
240  m_lu.template triangularView<UnitLower>().solveInPlace(dst);
241 
242  // Step 3
243  m_lu.template triangularView<Upper>().solveInPlace(dst);
244  }
245 
246  template<bool Conjugate, typename RhsType, typename DstType>
247  EIGEN_DEVICE_FUNC
248  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
249  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
250  * So we proceed as follows:
251  * Step 1: compute c = Pb.
252  * Step 2: replace c by the solution x to Lx = c.
253  * Step 3: replace c by the solution x to Ux = c.
254  */
255 
256  eigen_assert(rhs.rows() == m_lu.cols());
257 
258  if (Conjugate) {
259  // Step 1
260  dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs);
261  // Step 2
262  m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst);
263  } else {
264  // Step 1
265  dst = m_lu.template triangularView<Upper>().transpose().solve(rhs);
266  // Step 2
267  m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst);
268  }
269  // Step 3
270  dst = permutationP().transpose() * dst;
271  }
272  #endif
273 
274  protected:
275 
276  static void check_template_parameters()
277  {
278  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
279  }
280 
281  void compute();
282 
283  MatrixType m_lu;
284  PermutationType m_p;
285  TranspositionType m_rowsTranspositions;
286  RealScalar m_l1_norm;
287  signed char m_det_p;
288  bool m_isInitialized;
289 };
290 
291 template<typename MatrixType>
293  : m_lu(),
294  m_p(),
295  m_rowsTranspositions(),
296  m_l1_norm(0),
297  m_det_p(0),
298  m_isInitialized(false)
299 {
300 }
301 
302 template<typename MatrixType>
304  : m_lu(size, size),
305  m_p(size),
306  m_rowsTranspositions(size),
307  m_l1_norm(0),
308  m_det_p(0),
309  m_isInitialized(false)
310 {
311 }
312 
313 template<typename MatrixType>
314 template<typename InputType>
316  : m_lu(matrix.rows(),matrix.cols()),
317  m_p(matrix.rows()),
318  m_rowsTranspositions(matrix.rows()),
319  m_l1_norm(0),
320  m_det_p(0),
321  m_isInitialized(false)
322 {
323  compute(matrix.derived());
324 }
325 
326 template<typename MatrixType>
327 template<typename InputType>
329  : m_lu(matrix.derived()),
330  m_p(matrix.rows()),
331  m_rowsTranspositions(matrix.rows()),
332  m_l1_norm(0),
333  m_det_p(0),
334  m_isInitialized(false)
335 {
336  compute();
337 }
338 
339 namespace internal {
340 
342 template<typename Scalar, int StorageOrder, typename PivIndex>
343 struct partial_lu_impl
344 {
345  // FIXME add a stride to Map, so that the following mapping becomes easier,
346  // another option would be to create an expression being able to automatically
347  // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
348  // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
349  // and Block.
351  typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
352  typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
353  typedef typename MatrixType::RealScalar RealScalar;
354 
365  static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
366  {
367  typedef scalar_score_coeff_op<Scalar> Scoring;
368  typedef typename Scoring::result_type Score;
369  const Index rows = lu.rows();
370  const Index cols = lu.cols();
371  const Index size = (std::min)(rows,cols);
372  nb_transpositions = 0;
373  Index first_zero_pivot = -1;
374  for(Index k = 0; k < size; ++k)
375  {
376  Index rrows = rows-k-1;
377  Index rcols = cols-k-1;
378 
379  Index row_of_biggest_in_col;
380  Score biggest_in_corner
381  = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
382  row_of_biggest_in_col += k;
383 
384  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
385 
386  if(biggest_in_corner != Score(0))
387  {
388  if(k != row_of_biggest_in_col)
389  {
390  lu.row(k).swap(lu.row(row_of_biggest_in_col));
391  ++nb_transpositions;
392  }
393 
394  // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
395  // overflow but not the actual quotient?
396  lu.col(k).tail(rrows) /= lu.coeff(k,k);
397  }
398  else if(first_zero_pivot==-1)
399  {
400  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
401  // and continue the factorization such we still have A = PLU
402  first_zero_pivot = k;
403  }
404 
405  if(k<rows-1)
406  lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
407  }
408  return first_zero_pivot;
409  }
410 
426  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
427  {
428  MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
429  MatrixType lu(lu1,0,0,rows,cols);
430 
431  const Index size = (std::min)(rows,cols);
432 
433  // if the matrix is too small, no blocking:
434  if(size<=16)
435  {
436  return unblocked_lu(lu, row_transpositions, nb_transpositions);
437  }
438 
439  // automatically adjust the number of subdivisions to the size
440  // of the matrix so that there is enough sub blocks:
441  Index blockSize;
442  {
443  blockSize = size/8;
444  blockSize = (blockSize/16)*16;
445  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
446  }
447 
448  nb_transpositions = 0;
449  Index first_zero_pivot = -1;
450  for(Index k = 0; k < size; k+=blockSize)
451  {
452  Index bs = (std::min)(size-k,blockSize); // actual size of the block
453  Index trows = rows - k - bs; // trailing rows
454  Index tsize = size - k - bs; // trailing size
455 
456  // partition the matrix:
457  // A00 | A01 | A02
458  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
459  // A20 | A21 | A22
460  BlockType A_0(lu,0,0,rows,k);
461  BlockType A_2(lu,0,k+bs,rows,tsize);
462  BlockType A11(lu,k,k,bs,bs);
463  BlockType A12(lu,k,k+bs,bs,tsize);
464  BlockType A21(lu,k+bs,k,trows,bs);
465  BlockType A22(lu,k+bs,k+bs,trows,tsize);
466 
467  PivIndex nb_transpositions_in_panel;
468  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
469  // with a very small blocking size:
470  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
471  row_transpositions+k, nb_transpositions_in_panel, 16);
472  if(ret>=0 && first_zero_pivot==-1)
473  first_zero_pivot = k+ret;
474 
475  nb_transpositions += nb_transpositions_in_panel;
476  // update permutations and apply them to A_0
477  for(Index i=k; i<k+bs; ++i)
478  {
479  Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
480  A_0.row(i).swap(A_0.row(piv));
481  }
482 
483  if(trows)
484  {
485  // apply permutations to A_2
486  for(Index i=k;i<k+bs; ++i)
487  A_2.row(i).swap(A_2.row(row_transpositions[i]));
488 
489  // A12 = A11^-1 A12
490  A11.template triangularView<UnitLower>().solveInPlace(A12);
491 
492  A22.noalias() -= A21 * A12;
493  }
494  }
495  return first_zero_pivot;
496  }
497 };
498 
501 template<typename MatrixType, typename TranspositionType>
502 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
503 {
504  eigen_assert(lu.cols() == row_transpositions.size());
505  eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
506 
507  partial_lu_impl
508  <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndex>
509  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
510 }
511 
512 } // end namespace internal
513 
514 template<typename MatrixType>
515 void PartialPivLU<MatrixType>::compute()
516 {
517  check_template_parameters();
518 
519  // the row permutation is stored as int indices, so just to be sure:
520  eigen_assert(m_lu.rows()<NumTraits<int>::highest());
521 
522  if(m_lu.cols()>0)
523  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
524  else
525  m_l1_norm = RealScalar(0);
526 
527  eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
528  const Index size = m_lu.rows();
529 
530  m_rowsTranspositions.resize(size);
531 
532  typename TranspositionType::StorageIndex nb_transpositions;
533  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
534  m_det_p = (nb_transpositions%2) ? -1 : 1;
535 
536  m_p = m_rowsTranspositions;
537 
538  m_isInitialized = true;
539 }
540 
541 template<typename MatrixType>
542 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
543 {
544  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
545  return Scalar(m_det_p) * m_lu.diagonal().prod();
546 }
547 
551 template<typename MatrixType>
553 {
554  eigen_assert(m_isInitialized && "LU is not initialized.");
555  // LU
556  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
557  * m_lu.template triangularView<Upper>();
558 
559  // P^{-1}(LU)
560  res = m_p.inverse() * res;
561 
562  return res;
563 }
564 
565 /***** Implementation details *****************************************************/
566 
567 namespace internal {
568 
569 /***** Implementation of inverse() *****************************************************/
570 template<typename DstXprType, typename MatrixType>
571 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
572 {
573  typedef PartialPivLU<MatrixType> LuType;
574  typedef Inverse<LuType> SrcXprType;
575  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
576  {
577  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
578  }
579 };
580 } // end namespace internal
581 
582 /******** MatrixBase methods *******/
583 
590 template<typename Derived>
591 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
593 {
594  return PartialPivLU<PlainObject>(eval());
595 }
596 
605 template<typename Derived>
608 {
609  return PartialPivLU<PlainObject>(eval());
610 }
611 
612 } // end namespace Eigen
613 
614 #endif // EIGEN_PARTIALLU_H
Eigen::Inverse
Expression of the inverse of another expression.
Definition: Inverse.h:43
Eigen::MatrixBase::lu
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:607
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Eigen::PartialPivLU
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:246
Eigen::Block
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Eigen::EigenBase::Index
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
Eigen::PermutationBase::transpose
InverseReturnType transpose() const
Definition: PermutationMatrix.h:202
Eigen::EigenBase
Definition: EigenBase.h:29
Eigen::DenseCoeffsBase< Derived, DirectWriteAccessors >::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:61
Eigen::RowMajor
Definition: Constants.h:322
Eigen::PartialPivLU::inverse
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:197
Eigen::EigenBase::size
Index size() const
Definition: EigenBase.h:66
Eigen::PartialPivLU::matrixLU
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:141
Eigen::PartialPivLU::permutationP
const PermutationType & permutationP() const
Definition: PartialPivLU.h:149
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::PartialPivLU::determinant
Scalar determinant() const
Definition: PartialPivLU.h:542
Eigen::PartialPivLU::rcond
RealScalar rcond() const
Definition: PartialPivLU.h:184
Eigen::PartialPivLU::reconstructedMatrix
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:552
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::PartialPivLU::PartialPivLU
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:292
Eigen::Solve
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Eigen::PermutationMatrix< RowsAtCompileTime, MaxRowsAtCompileTime >
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::PartialPivLU::solve
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:175
Eigen::ColMajor
Definition: Constants.h:320
Eigen::MatrixBase::partialPivLu
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:592
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33