Eigen  3.4.0 (git rev e3e74001f7c4bf95f0dde572e8a08c5b2918a3ab)
HouseholderQR.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6 // Copyright (C) 2010 Vincent Lejeune
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 #ifndef EIGEN_QR_H
13 #define EIGEN_QR_H
14 
15 namespace Eigen {
16 
17 namespace internal {
18 template<typename _MatrixType> struct traits<HouseholderQR<_MatrixType> >
19  : traits<_MatrixType>
20 {
21  typedef MatrixXpr XprKind;
22  typedef SolverStorage StorageKind;
23  typedef int StorageIndex;
24  enum { Flags = 0 };
25 };
26 
27 } // end namespace internal
28 
56 template<typename _MatrixType> class HouseholderQR
57  : public SolverBase<HouseholderQR<_MatrixType> >
58 {
59  public:
60 
61  typedef _MatrixType MatrixType;
63  friend class SolverBase<HouseholderQR>;
64 
65  EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR)
66  enum {
67  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
68  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
69  };
70  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
71  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
72  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
74 
81  HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
82 
90  : m_qr(rows, cols),
91  m_hCoeffs((std::min)(rows,cols)),
92  m_temp(cols),
93  m_isInitialized(false) {}
94 
107  template<typename InputType>
108  explicit HouseholderQR(const EigenBase<InputType>& matrix)
109  : m_qr(matrix.rows(), matrix.cols()),
110  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
111  m_temp(matrix.cols()),
112  m_isInitialized(false)
113  {
114  compute(matrix.derived());
115  }
116 
117 
125  template<typename InputType>
127  : m_qr(matrix.derived()),
128  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
129  m_temp(matrix.cols()),
130  m_isInitialized(false)
131  {
132  computeInPlace();
133  }
134 
135  #ifdef EIGEN_PARSED_BY_DOXYGEN
150  template<typename Rhs>
151  inline const Solve<HouseholderQR, Rhs>
152  solve(const MatrixBase<Rhs>& b) const;
153  #endif
154 
164  {
165  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
166  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
167  }
168 
172  const MatrixType& matrixQR() const
173  {
174  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
175  return m_qr;
176  }
177 
178  template<typename InputType>
179  HouseholderQR& compute(const EigenBase<InputType>& matrix) {
180  m_qr = matrix.derived();
181  computeInPlace();
182  return *this;
183  }
184 
198  typename MatrixType::RealScalar absDeterminant() const;
199 
212  typename MatrixType::RealScalar logAbsDeterminant() const;
213 
214  inline Index rows() const { return m_qr.rows(); }
215  inline Index cols() const { return m_qr.cols(); }
216 
221  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
222 
223  #ifndef EIGEN_PARSED_BY_DOXYGEN
224  template<typename RhsType, typename DstType>
225  void _solve_impl(const RhsType &rhs, DstType &dst) const;
226 
227  template<bool Conjugate, typename RhsType, typename DstType>
228  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
229  #endif
230 
231  protected:
232 
233  static void check_template_parameters()
234  {
235  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
236  }
237 
238  void computeInPlace();
239 
240  MatrixType m_qr;
241  HCoeffsType m_hCoeffs;
242  RowVectorType m_temp;
243  bool m_isInitialized;
244 };
245 
246 template<typename MatrixType>
247 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
248 {
249  using std::abs;
250  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
251  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
252  return abs(m_qr.diagonal().prod());
253 }
254 
255 template<typename MatrixType>
256 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
257 {
258  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
259  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
260  return m_qr.diagonal().cwiseAbs().array().log().sum();
261 }
262 
263 namespace internal {
264 
266 template<typename MatrixQR, typename HCoeffs>
267 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
268 {
269  typedef typename MatrixQR::Scalar Scalar;
270  typedef typename MatrixQR::RealScalar RealScalar;
271  Index rows = mat.rows();
272  Index cols = mat.cols();
273  Index size = (std::min)(rows,cols);
274 
275  eigen_assert(hCoeffs.size() == size);
276 
278  TempType tempVector;
279  if(tempData==0)
280  {
281  tempVector.resize(cols);
282  tempData = tempVector.data();
283  }
284 
285  for(Index k = 0; k < size; ++k)
286  {
287  Index remainingRows = rows - k;
288  Index remainingCols = cols - k - 1;
289 
290  RealScalar beta;
291  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
292  mat.coeffRef(k,k) = beta;
293 
294  // apply H to remaining part of m_qr from the left
295  mat.bottomRightCorner(remainingRows, remainingCols)
296  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
297  }
298 }
299 
301 template<typename MatrixQR, typename HCoeffs,
302  typename MatrixQRScalar = typename MatrixQR::Scalar,
303  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
304 struct householder_qr_inplace_blocked
305 {
306  // This is specialized for LAPACK-supported Scalar types in HouseholderQR_LAPACKE.h
307  static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
308  typename MatrixQR::Scalar* tempData = 0)
309  {
310  typedef typename MatrixQR::Scalar Scalar;
311  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
312 
313  Index rows = mat.rows();
314  Index cols = mat.cols();
315  Index size = (std::min)(rows, cols);
316 
317  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
318  TempType tempVector;
319  if(tempData==0)
320  {
321  tempVector.resize(cols);
322  tempData = tempVector.data();
323  }
324 
325  Index blockSize = (std::min)(maxBlockSize,size);
326 
327  Index k = 0;
328  for (k = 0; k < size; k += blockSize)
329  {
330  Index bs = (std::min)(size-k,blockSize); // actual size of the block
331  Index tcols = cols - k - bs; // trailing columns
332  Index brows = rows-k; // rows of the block
333 
334  // partition the matrix:
335  // A00 | A01 | A02
336  // mat = A10 | A11 | A12
337  // A20 | A21 | A22
338  // and performs the qr dec of [A11^T A12^T]^T
339  // and update [A21^T A22^T]^T using level 3 operations.
340  // Finally, the algorithm continue on A22
341 
342  BlockType A11_21 = mat.block(k,k,brows,bs);
343  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
344 
345  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
346 
347  if(tcols)
348  {
349  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
350  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
351  }
352  }
353  }
354 };
355 
356 } // end namespace internal
357 
358 #ifndef EIGEN_PARSED_BY_DOXYGEN
359 template<typename _MatrixType>
360 template<typename RhsType, typename DstType>
361 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
362 {
363  const Index rank = (std::min)(rows(), cols());
364 
365  typename RhsType::PlainObject c(rhs);
366 
367  c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
368 
369  m_qr.topLeftCorner(rank, rank)
370  .template triangularView<Upper>()
371  .solveInPlace(c.topRows(rank));
372 
373  dst.topRows(rank) = c.topRows(rank);
374  dst.bottomRows(cols()-rank).setZero();
375 }
376 
377 template<typename _MatrixType>
378 template<bool Conjugate, typename RhsType, typename DstType>
379 void HouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
380 {
381  const Index rank = (std::min)(rows(), cols());
382 
383  typename RhsType::PlainObject c(rhs);
384 
385  m_qr.topLeftCorner(rank, rank)
386  .template triangularView<Upper>()
387  .transpose().template conjugateIf<Conjugate>()
388  .solveInPlace(c.topRows(rank));
389 
390  dst.topRows(rank) = c.topRows(rank);
391  dst.bottomRows(rows()-rank).setZero();
392 
393  dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
394 }
395 #endif
396 
403 template<typename MatrixType>
405 {
406  check_template_parameters();
407 
408  Index rows = m_qr.rows();
409  Index cols = m_qr.cols();
410  Index size = (std::min)(rows,cols);
411 
412  m_hCoeffs.resize(size);
413 
414  m_temp.resize(cols);
415 
416  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
417 
418  m_isInitialized = true;
419 }
420 
425 template<typename Derived>
428 {
429  return HouseholderQR<PlainObject>(eval());
430 }
431 
432 } // end namespace Eigen
433 
434 #endif // EIGEN_QR_H
Householder QR decomposition of a matrix.
Definition: HouseholderQR.h:58
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:89
void computeInPlace()
Definition: HouseholderQR.h:404
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:126
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:81
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:247
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:172
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:221
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:256
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:108
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:163
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:427
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Pseudo expression representing a solving operation.
Definition: Solve.h:63
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
HouseholderQR< _MatrixType > & derived()
Definition: EigenBase.h:46
@ ColMajor
Definition: Constants.h:319
@ RowMajor
Definition: Constants.h:321
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39