Eigen  3.2.92
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 } // end namespace internal
30 
62 template<typename _MatrixType> class PartialPivLU
63  : public SolverBase<PartialPivLU<_MatrixType> >
64 {
65  public:
66 
67  typedef _MatrixType MatrixType;
68  typedef SolverBase<PartialPivLU> Base;
69  EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
70  // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
71  enum {
72  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
73  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
74  };
75  typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
76  typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
77  typedef typename MatrixType::PlainObject PlainObject;
78 
79 
86  PartialPivLU();
87 
94  explicit PartialPivLU(Index size);
95 
103  template<typename InputType>
104  explicit PartialPivLU(const EigenBase<InputType>& matrix);
105 
106  template<typename InputType>
107  PartialPivLU& compute(const EigenBase<InputType>& 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  // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion.
147  template<typename Rhs>
148  inline const Solve<PartialPivLU, Rhs>
149  solve(const MatrixBase<Rhs>& b) const
150  {
151  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
152  return Solve<PartialPivLU, Rhs>(*this, b.derived());
153  }
154 
162  inline const Inverse<PartialPivLU> inverse() const
163  {
164  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
165  return Inverse<PartialPivLU>(*this);
166  }
167 
181  typename internal::traits<MatrixType>::Scalar determinant() const;
182 
183  MatrixType reconstructedMatrix() const;
184 
185  inline Index rows() const { return m_lu.rows(); }
186  inline Index cols() const { return m_lu.cols(); }
187 
188  #ifndef EIGEN_PARSED_BY_DOXYGEN
189  template<typename RhsType, typename DstType>
190  EIGEN_DEVICE_FUNC
191  void _solve_impl(const RhsType &rhs, DstType &dst) const {
192  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
193  * So we proceed as follows:
194  * Step 1: compute c = Pb.
195  * Step 2: replace c by the solution x to Lx = c.
196  * Step 3: replace c by the solution x to Ux = c.
197  */
198 
199  eigen_assert(rhs.rows() == m_lu.rows());
200 
201  // Step 1
202  dst = permutationP() * rhs;
203 
204  // Step 2
205  m_lu.template triangularView<UnitLower>().solveInPlace(dst);
206 
207  // Step 3
208  m_lu.template triangularView<Upper>().solveInPlace(dst);
209  }
210 
211  template<bool Conjugate, typename RhsType, typename DstType>
212  EIGEN_DEVICE_FUNC
213  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
214  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
215  * So we proceed as follows:
216  * Step 1: compute c = Pb.
217  * Step 2: replace c by the solution x to Lx = c.
218  * Step 3: replace c by the solution x to Ux = c.
219  */
220 
221  eigen_assert(rhs.rows() == m_lu.cols());
222 
223  if (Conjugate) {
224  // Step 1
225  dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs);
226  // Step 2
227  m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst);
228  } else {
229  // Step 1
230  dst = m_lu.template triangularView<Upper>().transpose().solve(rhs);
231  // Step 2
232  m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst);
233  }
234  // Step 3
235  dst = permutationP().transpose() * dst;
236  }
237  #endif
238 
239  protected:
240 
241  static void check_template_parameters()
242  {
243  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
244  }
245 
246  MatrixType m_lu;
247  PermutationType m_p;
248  TranspositionType m_rowsTranspositions;
249  Index m_det_p;
250  bool m_isInitialized;
251 };
252 
253 template<typename MatrixType>
255  : m_lu(),
256  m_p(),
257  m_rowsTranspositions(),
258  m_det_p(0),
259  m_isInitialized(false)
260 {
261 }
262 
263 template<typename MatrixType>
265  : m_lu(size, size),
266  m_p(size),
267  m_rowsTranspositions(size),
268  m_det_p(0),
269  m_isInitialized(false)
270 {
271 }
272 
273 template<typename MatrixType>
274 template<typename InputType>
276  : m_lu(matrix.rows(), matrix.rows()),
277  m_p(matrix.rows()),
278  m_rowsTranspositions(matrix.rows()),
279  m_det_p(0),
280  m_isInitialized(false)
281 {
282  compute(matrix.derived());
283 }
284 
285 namespace internal {
286 
288 template<typename Scalar, int StorageOrder, typename PivIndex>
289 struct partial_lu_impl
290 {
291  // FIXME add a stride to Map, so that the following mapping becomes easier,
292  // another option would be to create an expression being able to automatically
293  // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
294  // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
295  // and Block.
297  typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
298  typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
299  typedef typename MatrixType::RealScalar RealScalar;
300 
311  static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
312  {
313  typedef scalar_score_coeff_op<Scalar> Scoring;
314  typedef typename Scoring::result_type Score;
315  const Index rows = lu.rows();
316  const Index cols = lu.cols();
317  const Index size = (std::min)(rows,cols);
318  nb_transpositions = 0;
319  Index first_zero_pivot = -1;
320  for(Index k = 0; k < size; ++k)
321  {
322  Index rrows = rows-k-1;
323  Index rcols = cols-k-1;
324 
325  Index row_of_biggest_in_col;
326  Score biggest_in_corner
327  = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
328  row_of_biggest_in_col += k;
329 
330  row_transpositions[k] = PivIndex(row_of_biggest_in_col);
331 
332  if(biggest_in_corner != Score(0))
333  {
334  if(k != row_of_biggest_in_col)
335  {
336  lu.row(k).swap(lu.row(row_of_biggest_in_col));
337  ++nb_transpositions;
338  }
339 
340  // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
341  // overflow but not the actual quotient?
342  lu.col(k).tail(rrows) /= lu.coeff(k,k);
343  }
344  else if(first_zero_pivot==-1)
345  {
346  // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
347  // and continue the factorization such we still have A = PLU
348  first_zero_pivot = k;
349  }
350 
351  if(k<rows-1)
352  lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
353  }
354  return first_zero_pivot;
355  }
356 
372  static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
373  {
374  MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
375  MatrixType lu(lu1,0,0,rows,cols);
376 
377  const Index size = (std::min)(rows,cols);
378 
379  // if the matrix is too small, no blocking:
380  if(size<=16)
381  {
382  return unblocked_lu(lu, row_transpositions, nb_transpositions);
383  }
384 
385  // automatically adjust the number of subdivisions to the size
386  // of the matrix so that there is enough sub blocks:
387  Index blockSize;
388  {
389  blockSize = size/8;
390  blockSize = (blockSize/16)*16;
391  blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
392  }
393 
394  nb_transpositions = 0;
395  Index first_zero_pivot = -1;
396  for(Index k = 0; k < size; k+=blockSize)
397  {
398  Index bs = (std::min)(size-k,blockSize); // actual size of the block
399  Index trows = rows - k - bs; // trailing rows
400  Index tsize = size - k - bs; // trailing size
401 
402  // partition the matrix:
403  // A00 | A01 | A02
404  // lu = A_0 | A_1 | A_2 = A10 | A11 | A12
405  // A20 | A21 | A22
406  BlockType A_0(lu,0,0,rows,k);
407  BlockType A_2(lu,0,k+bs,rows,tsize);
408  BlockType A11(lu,k,k,bs,bs);
409  BlockType A12(lu,k,k+bs,bs,tsize);
410  BlockType A21(lu,k+bs,k,trows,bs);
411  BlockType A22(lu,k+bs,k+bs,trows,tsize);
412 
413  PivIndex nb_transpositions_in_panel;
414  // recursively call the blocked LU algorithm on [A11^T A21^T]^T
415  // with a very small blocking size:
416  Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
417  row_transpositions+k, nb_transpositions_in_panel, 16);
418  if(ret>=0 && first_zero_pivot==-1)
419  first_zero_pivot = k+ret;
420 
421  nb_transpositions += nb_transpositions_in_panel;
422  // update permutations and apply them to A_0
423  for(Index i=k; i<k+bs; ++i)
424  {
425  Index piv = (row_transpositions[i] += k);
426  A_0.row(i).swap(A_0.row(piv));
427  }
428 
429  if(trows)
430  {
431  // apply permutations to A_2
432  for(Index i=k;i<k+bs; ++i)
433  A_2.row(i).swap(A_2.row(row_transpositions[i]));
434 
435  // A12 = A11^-1 A12
436  A11.template triangularView<UnitLower>().solveInPlace(A12);
437 
438  A22.noalias() -= A21 * A12;
439  }
440  }
441  return first_zero_pivot;
442  }
443 };
444 
447 template<typename MatrixType, typename TranspositionType>
448 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
449 {
450  eigen_assert(lu.cols() == row_transpositions.size());
451  eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
452 
453  partial_lu_impl
454  <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndex>
455  ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
456 }
457 
458 } // end namespace internal
459 
460 template<typename MatrixType>
461 template<typename InputType>
462 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const EigenBase<InputType>& matrix)
463 {
464  check_template_parameters();
465 
466  // the row permutation is stored as int indices, so just to be sure:
467  eigen_assert(matrix.rows()<NumTraits<int>::highest());
468 
469  m_lu = matrix.derived();
470 
471  eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
472  const Index size = matrix.rows();
473 
474  m_rowsTranspositions.resize(size);
475 
476  typename TranspositionType::StorageIndex nb_transpositions;
477  internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
478  m_det_p = (nb_transpositions%2) ? -1 : 1;
479 
480  m_p = m_rowsTranspositions;
481 
482  m_isInitialized = true;
483  return *this;
484 }
485 
486 template<typename MatrixType>
487 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
488 {
489  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
490  return Scalar(m_det_p) * m_lu.diagonal().prod();
491 }
492 
496 template<typename MatrixType>
498 {
499  eigen_assert(m_isInitialized && "LU is not initialized.");
500  // LU
501  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
502  * m_lu.template triangularView<Upper>();
503 
504  // P^{-1}(LU)
505  res = m_p.inverse() * res;
506 
507  return res;
508 }
509 
510 /***** Implementation details *****************************************************/
511 
512 namespace internal {
513 
514 /***** Implementation of inverse() *****************************************************/
515 template<typename DstXprType, typename MatrixType, typename Scalar>
516 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
517 {
518  typedef PartialPivLU<MatrixType> LuType;
519  typedef Inverse<LuType> SrcXprType;
520  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
521  {
522  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
523  }
524 };
525 } // end namespace internal
526 
527 /******** MatrixBase methods *******/
528 
535 #ifndef __CUDACC__
536 template<typename Derived>
537 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
539 {
540  return PartialPivLU<PlainObject>(eval());
541 }
542 #endif
543 
552 #ifndef __CUDACC__
553 template<typename Derived>
556 {
557  return PartialPivLU<PlainObject>(eval());
558 }
559 #endif
560 
561 } // end namespace Eigen
562 
563 #endif // EIGEN_PARTIALLU_H
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:162
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:254
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
AdjointReturnType adjoint() const
Definition: SolverBase.h:109
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:248
Definition: Constants.h:320
Definition: LDLT.h:16
Derived & derived()
Definition: EigenBase.h:44
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
const unsigned int RowMajorBit
Definition: Constants.h:61
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:538
Definition: EigenBase.h:28
Expression of the inverse of another expression.
Definition: Inverse.h:43
Definition: Constants.h:322
const PermutationType & permutationP() const
Definition: PartialPivLU.h:123
InverseReturnType transpose() const
Definition: PermutationMatrix.h:203
ConstTransposeReturnType transpose() const
Definition: SolverBase.h:90
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:487
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:497
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:115
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:149
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:555
Index size() const
Definition: EigenBase.h:65