Eigen  3.3.7
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 
44 template<typename _MatrixType> class HouseholderQR
45 {
46  public:
47 
48  typedef _MatrixType MatrixType;
49  enum {
50  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
51  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
52  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
53  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
54  };
55  typedef typename MatrixType::Scalar Scalar;
56  typedef typename MatrixType::RealScalar RealScalar;
57  // FIXME should be int
58  typedef typename MatrixType::StorageIndex StorageIndex;
59  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
60  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
61  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
62  typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
63 
70  HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
71 
79  : m_qr(rows, cols),
80  m_hCoeffs((std::min)(rows,cols)),
81  m_temp(cols),
82  m_isInitialized(false) {}
83 
96  template<typename InputType>
97  explicit HouseholderQR(const EigenBase<InputType>& matrix)
98  : m_qr(matrix.rows(), matrix.cols()),
99  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
100  m_temp(matrix.cols()),
101  m_isInitialized(false)
102  {
103  compute(matrix.derived());
104  }
105 
106 
114  template<typename InputType>
116  : m_qr(matrix.derived()),
117  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
118  m_temp(matrix.cols()),
119  m_isInitialized(false)
120  {
121  computeInPlace();
122  }
123 
138  template<typename Rhs>
139  inline const Solve<HouseholderQR, Rhs>
140  solve(const MatrixBase<Rhs>& b) const
141  {
142  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
143  return Solve<HouseholderQR, Rhs>(*this, b.derived());
144  }
145 
155  {
156  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
157  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
158  }
159 
163  const MatrixType& matrixQR() const
164  {
165  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
166  return m_qr;
167  }
168 
169  template<typename InputType>
170  HouseholderQR& compute(const EigenBase<InputType>& matrix) {
171  m_qr = matrix.derived();
172  computeInPlace();
173  return *this;
174  }
175 
189  typename MatrixType::RealScalar absDeterminant() const;
190 
203  typename MatrixType::RealScalar logAbsDeterminant() const;
204 
205  inline Index rows() const { return m_qr.rows(); }
206  inline Index cols() const { return m_qr.cols(); }
207 
212  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
213 
214  #ifndef EIGEN_PARSED_BY_DOXYGEN
215  template<typename RhsType, typename DstType>
216  EIGEN_DEVICE_FUNC
217  void _solve_impl(const RhsType &rhs, DstType &dst) const;
218  #endif
219 
220  protected:
221 
222  static void check_template_parameters()
223  {
224  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
225  }
226 
227  void computeInPlace();
228 
229  MatrixType m_qr;
230  HCoeffsType m_hCoeffs;
231  RowVectorType m_temp;
232  bool m_isInitialized;
233 };
234 
235 template<typename MatrixType>
236 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
237 {
238  using std::abs;
239  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
240  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
241  return abs(m_qr.diagonal().prod());
242 }
243 
244 template<typename MatrixType>
245 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
246 {
247  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
248  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
249  return m_qr.diagonal().cwiseAbs().array().log().sum();
250 }
251 
252 namespace internal {
253 
255 template<typename MatrixQR, typename HCoeffs>
256 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
257 {
258  typedef typename MatrixQR::Scalar Scalar;
259  typedef typename MatrixQR::RealScalar RealScalar;
260  Index rows = mat.rows();
261  Index cols = mat.cols();
262  Index size = (std::min)(rows,cols);
263 
264  eigen_assert(hCoeffs.size() == size);
265 
267  TempType tempVector;
268  if(tempData==0)
269  {
270  tempVector.resize(cols);
271  tempData = tempVector.data();
272  }
273 
274  for(Index k = 0; k < size; ++k)
275  {
276  Index remainingRows = rows - k;
277  Index remainingCols = cols - k - 1;
278 
279  RealScalar beta;
280  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
281  mat.coeffRef(k,k) = beta;
282 
283  // apply H to remaining part of m_qr from the left
284  mat.bottomRightCorner(remainingRows, remainingCols)
285  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
286  }
287 }
288 
290 template<typename MatrixQR, typename HCoeffs,
291  typename MatrixQRScalar = typename MatrixQR::Scalar,
292  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
293 struct householder_qr_inplace_blocked
294 {
295  // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
296  static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
297  typename MatrixQR::Scalar* tempData = 0)
298  {
299  typedef typename MatrixQR::Scalar Scalar;
300  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
301 
302  Index rows = mat.rows();
303  Index cols = mat.cols();
304  Index size = (std::min)(rows, cols);
305 
306  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
307  TempType tempVector;
308  if(tempData==0)
309  {
310  tempVector.resize(cols);
311  tempData = tempVector.data();
312  }
313 
314  Index blockSize = (std::min)(maxBlockSize,size);
315 
316  Index k = 0;
317  for (k = 0; k < size; k += blockSize)
318  {
319  Index bs = (std::min)(size-k,blockSize); // actual size of the block
320  Index tcols = cols - k - bs; // trailing columns
321  Index brows = rows-k; // rows of the block
322 
323  // partition the matrix:
324  // A00 | A01 | A02
325  // mat = A10 | A11 | A12
326  // A20 | A21 | A22
327  // and performs the qr dec of [A11^T A12^T]^T
328  // and update [A21^T A22^T]^T using level 3 operations.
329  // Finally, the algorithm continue on A22
330 
331  BlockType A11_21 = mat.block(k,k,brows,bs);
332  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
333 
334  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
335 
336  if(tcols)
337  {
338  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
339  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
340  }
341  }
342  }
343 };
344 
345 } // end namespace internal
346 
347 #ifndef EIGEN_PARSED_BY_DOXYGEN
348 template<typename _MatrixType>
349 template<typename RhsType, typename DstType>
350 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
351 {
352  const Index rank = (std::min)(rows(), cols());
353  eigen_assert(rhs.rows() == rows());
354 
355  typename RhsType::PlainObject c(rhs);
356 
357  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
358  c.applyOnTheLeft(householderSequence(
359  m_qr.leftCols(rank),
360  m_hCoeffs.head(rank)).transpose()
361  );
362 
363  m_qr.topLeftCorner(rank, rank)
364  .template triangularView<Upper>()
365  .solveInPlace(c.topRows(rank));
366 
367  dst.topRows(rank) = c.topRows(rank);
368  dst.bottomRows(cols()-rank).setZero();
369 }
370 #endif
371 
378 template<typename MatrixType>
380 {
381  check_template_parameters();
382 
383  Index rows = m_qr.rows();
384  Index cols = m_qr.cols();
385  Index size = (std::min)(rows,cols);
386 
387  m_hCoeffs.resize(size);
388 
389  m_temp.resize(cols);
390 
391  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
392 
393  m_isInitialized = true;
394 }
395 
400 template<typename Derived>
403 {
404  return HouseholderQR<PlainObject>(eval());
405 }
406 
407 } // end namespace Eigen
408 
409 #endif // EIGEN_QR_H
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:78
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:163
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:236
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Definition: Half.h:150
Derived & derived()
Definition: EigenBase.h:45
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: HouseholderQR.h:140
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:451
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Definition: EigenBase.h:29
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: ForwardDeclarations.h:258
Derived & derived()
Definition: EigenBase.h:45
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:115
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:245
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:402
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:154
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
Householder QR decomposition of a matrix.
Definition: ForwardDeclarations.h:250
void computeInPlace()
Definition: HouseholderQR.h:379
Pseudo expression representing a solving operation.
Definition: Solve.h:62
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:70
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:97
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:212