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;
23  enum {
24  Flags = BaseTraits::Flags & RowMajorBit,
25  CoeffReadCost = Dynamic
26  };
27 };
28 
29 template<typename T,typename Derived>
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;
83  // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int
84  enum {
85  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
86  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
87  };
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>
127 
128  template<typename InputType>
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 
277  {
279  }
280 
281  void compute();
282 
283  MatrixType m_lu;
284  PermutationType m_p;
285  TranspositionType m_rowsTranspositions;
287  signed char m_det_p;
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>
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.
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));
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>
503 {
504  eigen_assert(lu.cols() == row_transpositions.size());
505  eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
506 
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>
516 {
518 
519  // the row permutation is stored as int indices, so just to be sure:
521 
522  m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
523 
524  eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
525  const Index size = m_lu.rows();
526 
528 
531  m_det_p = (nb_transpositions%2) ? -1 : 1;
532 
534 
535  m_isInitialized = true;
536 }
537 
538 template<typename MatrixType>
540 {
541  eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
542  return Scalar(m_det_p) * m_lu.diagonal().prod();
543 }
544 
548 template<typename MatrixType>
550 {
551  eigen_assert(m_isInitialized && "LU is not initialized.");
552  // LU
553  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
554  * m_lu.template triangularView<Upper>();
555 
556  // P^{-1}(LU)
557  res = m_p.inverse() * res;
558 
559  return res;
560 }
561 
562 /***** Implementation details *****************************************************/
563 
564 namespace internal {
565 
566 /***** Implementation of inverse() *****************************************************/
567 template<typename DstXprType, typename MatrixType>
568 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
569 {
572  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
573  {
574  dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
575  }
576 };
577 } // end namespace internal
578 
579 /******** MatrixBase methods *******/
580 
587 template<typename Derived>
590 {
592 }
593 
602 template<typename Derived>
605 {
607 }
608 
609 } // end namespace Eigen
610 
611 #endif // EIGEN_PARTIALLU_H
SCALAR Scalar
Definition: bench_gemm.cpp:33
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:852
#define max(a, b)
Definition: datatypes.h:20
RealScalar rcond() const
Definition: PartialPivLU.h:184
MatrixType::PlainObject PlainObject
Definition: PartialPivLU.h:90
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:197
PartialPivLU & compute(const EigenBase< InputType > &matrix)
Definition: PartialPivLU.h:129
EIGEN_DEVICE_FUNC Index rows() const
Definition: Inverse.h:58
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:292
internal::traits< PartialPivLU< Matrix< Scalar, Dynamic, Dynamic > > >::Scalar Scalar
Definition: SolverBase.h:46
PermutationType m_p
Definition: PartialPivLU.h:284
Block< MatrixType, Dynamic, Dynamic > BlockType
Definition: PartialPivLU.h:352
Map< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > > MapLU
Definition: PartialPivLU.h:350
#define min(a, b)
Definition: datatypes.h:19
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
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:572
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:150
EIGEN_DEVICE_FUNC void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
Definition: PartialPivLU.h:248
Transpositions< RowsAtCompileTime, MaxRowsAtCompileTime > TranspositionType
Definition: PartialPivLU.h:89
MatrixType::RealScalar RealScalar
Definition: PartialPivLU.h:353
signed char m_det_p
Definition: PartialPivLU.h:287
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:38
const unsigned int RowMajorBit
Definition: Constants.h:61
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:589
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
Scalar determinant() const
Definition: PartialPivLU.h:539
void partial_lu_inplace(MatrixType &lu, TranspositionType &row_transpositions, typename TranspositionType::StorageIndex &nb_transpositions)
Definition: PartialPivLU.h:502
SolverBase< PartialPivLU > Base
Definition: PartialPivLU.h:81
const PermutationType & permutationP() const
Definition: PartialPivLU.h:149
static void check_template_parameters()
Definition: PartialPivLU.h:276
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
EIGEN_DEVICE_FUNC Index cols() const
Definition: Inverse.h:59
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:184
void resize(Index newSize)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Eigen::FullPivLU< Matrix5x3 > lu(m)
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:192
InverseReturnType inverse() const
PermutationMatrix< RowsAtCompileTime, MaxRowsAtCompileTime > PermutationType
Definition: PartialPivLU.h:88
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
DenseIndex ret
Definition: level1_impl.h:59
Index cols() const
Definition: PartialPivLU.h:221
static Index unblocked_lu(MatrixType &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions)
Definition: PartialPivLU.h:365
Block< MapLU, Dynamic, Dynamic > MatrixType
Definition: PartialPivLU.h:351
TranspositionType m_rowsTranspositions
Definition: PartialPivLU.h:285
EIGEN_DEVICE_FUNC void _solve_impl(const RhsType &rhs, DstType &dst) const
Definition: PartialPivLU.h:226
internal::nested_eval< T, 1 >::type eval(const T &xpr)
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:549
Index rows() const
Definition: PartialPivLU.h:220
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
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)
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:141
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:175
EIGEN_DEVICE_FUNC const XprTypeNestedCleaned & nestedExpression() const
Definition: Inverse.h:61
RealScalar m_l1_norm
Definition: PartialPivLU.h:286
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:41
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:426
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
_MatrixType MatrixType
Definition: PartialPivLU.h:80
EIGEN_DEVICE_FUNC Index size() const
Definition: EigenBase.h:66
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:604


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:23