Eigen  3.2.91
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 traits<_MatrixType> BaseTraits;
21  enum {
22  Flags = BaseTraits::Flags & RowMajorBit,
23  CoeffReadCost = Dynamic
24  };
25 };
26 
27 } // end namespace internal
28 
60 template<typename _MatrixType> class PartialPivLU
61 {
62  public:
63 
64  typedef _MatrixType MatrixType;
65  enum {
66  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68  Options = MatrixType::Options,
69  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71  };
72  typedef typename MatrixType::Scalar Scalar;
73  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
74  typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
75  // FIXME should be int
76  typedef typename MatrixType::StorageIndex StorageIndex;
77  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
78  typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
79  typedef typename MatrixType::PlainObject PlainObject;
80 
81 
88  PartialPivLU();
89 
96  explicit PartialPivLU(Index size);
97 
105  explicit PartialPivLU(const MatrixType& matrix);
106 
107  PartialPivLU& compute(const MatrixType& matrix);
108 
115  inline const MatrixType& matrixLU() const
116  {
117  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
118  return m_lu;
119  }
120 
123  inline const PermutationType& permutationP() const
124  {
125  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
126  return m_p;
127  }
128 
146  template<typename Rhs>
147  inline const Solve<PartialPivLU, Rhs>
148  solve(const MatrixBase<Rhs>& b) const
149  {
150  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
151  return Solve<PartialPivLU, Rhs>(*this, b.derived());
152  }
153 
161  inline const Inverse<PartialPivLU> inverse() const
162  {
163  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
164  return Inverse<PartialPivLU>(*this);
165  }
166 
180  typename internal::traits<MatrixType>::Scalar determinant() const;
181 
182  MatrixType reconstructedMatrix() const;
183 
184  inline Index rows() const { return m_lu.rows(); }
185  inline Index cols() const { return m_lu.cols(); }
186 
187  #ifndef EIGEN_PARSED_BY_DOXYGEN
188  template<typename RhsType, typename DstType>
189  EIGEN_DEVICE_FUNC
190  void _solve_impl(const RhsType &rhs, DstType &dst) const {
191  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
192  * So we proceed as follows:
193  * Step 1: compute c = Pb.
194  * Step 2: replace c by the solution x to Lx = c.
195  * Step 3: replace c by the solution x to Ux = c.
196  */
197 
198  eigen_assert(rhs.rows() == m_lu.rows());
199 
200  // Step 1
201  dst = permutationP() * rhs;
202 
203  // Step 2
204  m_lu.template triangularView<UnitLower>().solveInPlace(dst);
205 
206  // Step 3
207  m_lu.template triangularView<Upper>().solveInPlace(dst);
208  }
209  #endif
210 
211  protected:
212 
213  static void check_template_parameters()
214  {
215  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
216  }
217 
218  MatrixType m_lu;
219  PermutationType m_p;
220  TranspositionType m_rowsTranspositions;
221  Index m_det_p;
222  bool m_isInitialized;
223 };
224 
225 template<typename MatrixType>
227  : m_lu(),
228  m_p(),
229  m_rowsTranspositions(),
230  m_det_p(0),
231  m_isInitialized(false)
232 {
233 }
234 
235 template<typename MatrixType>
237  : m_lu(size, size),
238  m_p(size),
239  m_rowsTranspositions(size),
240  m_det_p(0),
241  m_isInitialized(false)
242 {
243 }
244 
245 template<typename MatrixType>
246 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
247  : m_lu(matrix.rows(), matrix.rows()),
248  m_p(matrix.rows()),
249  m_rowsTranspositions(matrix.rows()),
250  m_det_p(0),
251  m_isInitialized(false)
252 {
253  compute(matrix);
254 }
255 
256 namespace internal {
257 
259 template<typename Scalar, int StorageOrder, typename PivIndex>
260 struct partial_lu_impl
261 {
262  // FIXME add a stride to Map, so that the following mapping becomes easier,
263  // another option would be to create an expression being able to automatically
264  // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
265  // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
266  // and Block.
268  typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
269  typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
270  typedef typename MatrixType::RealScalar RealScalar;
271 
282  static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
283  {
284  typedef scalar_score_coeff_op<Scalar> Scoring;
285  typedef typename Scoring::result_type Score;
286  const Index rows = lu.rows();
287  const Index cols = lu.cols();
288  const Index size = (std::min)(rows,cols);
289  nb_transpositions = 0;
290  Index first_zero_pivot = -1;
291  for(Index k = 0; k < size; ++k)
292  {
293  Index rrows = rows-k-1;
294  Index rcols = cols-k-1;
295 
296  Index row_of_biggest_in_col;
297  Score biggest_in_corner
298  = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
299  row_of_biggest_in_col += k;
300 
301  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
302 
303  if(biggest_in_corner != Score(0))
304  {
305  if(k != row_of_biggest_in_col)
306  {
307  lu.row(k).swap(lu.row(row_of_biggest_in_col));
308  ++nb_transpositions;
309  }
310 
311  // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
312  // overflow but not the actual quotient?
313  lu.col(k).tail(rrows) /= lu.coeff(k,k);
314  }
315  else if(first_zero_pivot==-1)
316  {
317  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
318  // and continue the factorization such we still have A = PLU
319  first_zero_pivot = k;
320  }
321 
322  if(k<rows-1)
323  lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
324  }
325  return first_zero_pivot;
326  }
327 
343  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
344  {
345  MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
346  MatrixType lu(lu1,0,0,rows,cols);
347 
348  const Index size = (std::min)(rows,cols);
349 
350  // if the matrix is too small, no blocking:
351  if(size<=16)
352  {
353  return unblocked_lu(lu, row_transpositions, nb_transpositions);
354  }
355 
356  // automatically adjust the number of subdivisions to the size
357  // of the matrix so that there is enough sub blocks:
358  Index blockSize;
359  {
360  blockSize = size/8;
361  blockSize = (blockSize/16)*16;
362  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
363  }
364 
365  nb_transpositions = 0;
366  Index first_zero_pivot = -1;
367  for(Index k = 0; k < size; k+=blockSize)
368  {
369  Index bs = (std::min)(size-k,blockSize); // actual size of the block
370  Index trows = rows - k - bs; // trailing rows
371  Index tsize = size - k - bs; // trailing size
372 
373  // partition the matrix:
374  // A00 | A01 | A02
375  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
376  // A20 | A21 | A22
377  BlockType A_0(lu,0,0,rows,k);
378  BlockType A_2(lu,0,k+bs,rows,tsize);
379  BlockType A11(lu,k,k,bs,bs);
380  BlockType A12(lu,k,k+bs,bs,tsize);
381  BlockType A21(lu,k+bs,k,trows,bs);
382  BlockType A22(lu,k+bs,k+bs,trows,tsize);
383 
384  PivIndex nb_transpositions_in_panel;
385  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
386  // with a very small blocking size:
387  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
388  row_transpositions+k, nb_transpositions_in_panel, 16);
389  if(ret>=0 && first_zero_pivot==-1)
390  first_zero_pivot = k+ret;
391 
392  nb_transpositions += nb_transpositions_in_panel;
393  // update permutations and apply them to A_0
394  for(Index i=k; i<k+bs; ++i)
395  {
396  Index piv = (row_transpositions[i] += k);
397  A_0.row(i).swap(A_0.row(piv));
398  }
399 
400  if(trows)
401  {
402  // apply permutations to A_2
403  for(Index i=k;i<k+bs; ++i)
404  A_2.row(i).swap(A_2.row(row_transpositions[i]));
405 
406  // A12 = A11^-1 A12
407  A11.template triangularView<UnitLower>().solveInPlace(A12);
408 
409  A22.noalias() -= A21 * A12;
410  }
411  }
412  return first_zero_pivot;
413  }
414 };
415 
418 template<typename MatrixType, typename TranspositionType>
419 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
420 {
421  eigen_assert(lu.cols() == row_transpositions.size());
422  eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
423 
424  partial_lu_impl
425  <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndex>
426  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
427 }
428 
429 } // end namespace internal
430 
431 template<typename MatrixType>
432 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
433 {
434  check_template_parameters();
435 
436  // the row permutation is stored as int indices, so just to be sure:
437  eigen_assert(matrix.rows()<NumTraits<int>::highest());
438 
439  m_lu = matrix;
440 
441  eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
442  const Index size = matrix.rows();
443 
444  m_rowsTranspositions.resize(size);
445 
446  typename TranspositionType::StorageIndex nb_transpositions;
447  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
448  m_det_p = (nb_transpositions%2) ? -1 : 1;
449 
450  m_p = m_rowsTranspositions;
451 
452  m_isInitialized = true;
453  return *this;
454 }
455 
456 template<typename MatrixType>
457 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
458 {
459  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
460  return Scalar(m_det_p) * m_lu.diagonal().prod();
461 }
462 
466 template<typename MatrixType>
468 {
469  eigen_assert(m_isInitialized && "LU is not initialized.");
470  // LU
471  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
472  * m_lu.template triangularView<Upper>();
473 
474  // P^{-1}(LU)
475  res = m_p.inverse() * res;
476 
477  return res;
478 }
479 
480 /***** Implementation of solve() *****************************************************/
481 
482 namespace internal {
483 
484 /***** Implementation of inverse() *****************************************************/
485 template<typename DstXprType, typename MatrixType, typename Scalar>
486 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
487 {
488  typedef PartialPivLU<MatrixType> LuType;
489  typedef Inverse<LuType> SrcXprType;
490  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
491  {
492  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
493  }
494 };
495 } // end namespace internal
496 
497 /******** MatrixBase methods *******/
498 
505 #ifndef __CUDACC__
506 template<typename Derived>
507 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
509 {
510  return PartialPivLU<PlainObject>(eval());
511 }
512 #endif
513 
522 #ifndef __CUDACC__
523 template<typename Derived>
526 {
527  return PartialPivLU<PlainObject>(eval());
528 }
529 #endif
530 
531 } // end namespace Eigen
532 
533 #endif // EIGEN_PARTIALLU_H
Definition: Constants.h:314
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:161
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:226
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:247
Definition: LDLT.h:16
const unsigned int RowMajorBit
Definition: Constants.h:53
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:508
Expression of the inverse of another expression.
Definition: Inverse.h:45
const PermutationType & permutationP() const
Definition: PartialPivLU.h:123
Definition: Eigen_Colamd.h:54
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:104
internal::traits< MatrixType >::Scalar determinant() const
Definition: PartialPivLU.h:457
Definition: Constants.h:312
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:467
Pseudo expression representing a solving operation.
Definition: Solve.h:63
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:115
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:148
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:525