Eigen  3.2.93
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 
141  template<typename Rhs>
142  inline const Solve<HouseholderQR, Rhs>
143  solve(const MatrixBase<Rhs>& b) const
144  {
145  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
146  return Solve<HouseholderQR, Rhs>(*this, b.derived());
147  }
148 
157  HouseholderSequenceType householderQ() const
158  {
159  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
160  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
161  }
162 
166  const MatrixType& matrixQR() const
167  {
168  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
169  return m_qr;
170  }
171 
172  template<typename InputType>
173  HouseholderQR& compute(const EigenBase<InputType>& matrix) {
174  m_qr = matrix.derived();
175  computeInPlace();
176  return *this;
177  }
178 
192  typename MatrixType::RealScalar absDeterminant() const;
193 
206  typename MatrixType::RealScalar logAbsDeterminant() const;
207 
208  inline Index rows() const { return m_qr.rows(); }
209  inline Index cols() const { return m_qr.cols(); }
210 
215  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
216 
217  #ifndef EIGEN_PARSED_BY_DOXYGEN
218  template<typename RhsType, typename DstType>
219  EIGEN_DEVICE_FUNC
220  void _solve_impl(const RhsType &rhs, DstType &dst) const;
221  #endif
222 
223  protected:
224 
225  static void check_template_parameters()
226  {
227  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
228  }
229 
230  void computeInPlace();
231 
232  MatrixType m_qr;
233  HCoeffsType m_hCoeffs;
234  RowVectorType m_temp;
235  bool m_isInitialized;
236 };
237 
238 template<typename MatrixType>
239 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
240 {
241  using std::abs;
242  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
243  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
244  return abs(m_qr.diagonal().prod());
245 }
246 
247 template<typename MatrixType>
248 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
249 {
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 m_qr.diagonal().cwiseAbs().array().log().sum();
253 }
254 
255 namespace internal {
256 
258 template<typename MatrixQR, typename HCoeffs>
259 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
260 {
261  typedef typename MatrixQR::Scalar Scalar;
262  typedef typename MatrixQR::RealScalar RealScalar;
263  Index rows = mat.rows();
264  Index cols = mat.cols();
265  Index size = (std::min)(rows,cols);
266 
267  eigen_assert(hCoeffs.size() == size);
268 
270  TempType tempVector;
271  if(tempData==0)
272  {
273  tempVector.resize(cols);
274  tempData = tempVector.data();
275  }
276 
277  for(Index k = 0; k < size; ++k)
278  {
279  Index remainingRows = rows - k;
280  Index remainingCols = cols - k - 1;
281 
282  RealScalar beta;
283  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
284  mat.coeffRef(k,k) = beta;
285 
286  // apply H to remaining part of m_qr from the left
287  mat.bottomRightCorner(remainingRows, remainingCols)
288  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
289  }
290 }
291 
293 template<typename MatrixQR, typename HCoeffs,
294  typename MatrixQRScalar = typename MatrixQR::Scalar,
295  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
296 struct householder_qr_inplace_blocked
297 {
298  // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
299  static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
300  typename MatrixQR::Scalar* tempData = 0)
301  {
302  typedef typename MatrixQR::Scalar Scalar;
303  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
304 
305  Index rows = mat.rows();
306  Index cols = mat.cols();
307  Index size = (std::min)(rows, cols);
308 
310  TempType tempVector;
311  if(tempData==0)
312  {
313  tempVector.resize(cols);
314  tempData = tempVector.data();
315  }
316 
317  Index blockSize = (std::min)(maxBlockSize,size);
318 
319  Index k = 0;
320  for (k = 0; k < size; k += blockSize)
321  {
322  Index bs = (std::min)(size-k,blockSize); // actual size of the block
323  Index tcols = cols - k - bs; // trailing columns
324  Index brows = rows-k; // rows of the block
325 
326  // partition the matrix:
327  // A00 | A01 | A02
328  // mat = A10 | A11 | A12
329  // A20 | A21 | A22
330  // and performs the qr dec of [A11^T A12^T]^T
331  // and update [A21^T A22^T]^T using level 3 operations.
332  // Finally, the algorithm continue on A22
333 
334  BlockType A11_21 = mat.block(k,k,brows,bs);
335  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
336 
337  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
338 
339  if(tcols)
340  {
341  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
342  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
343  }
344  }
345  }
346 };
347 
348 } // end namespace internal
349 
350 #ifndef EIGEN_PARSED_BY_DOXYGEN
351 template<typename _MatrixType>
352 template<typename RhsType, typename DstType>
353 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
354 {
355  const Index rank = (std::min)(rows(), cols());
356  eigen_assert(rhs.rows() == rows());
357 
358  typename RhsType::PlainObject c(rhs);
359 
360  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
361  c.applyOnTheLeft(householderSequence(
362  m_qr.leftCols(rank),
363  m_hCoeffs.head(rank)).transpose()
364  );
365 
366  m_qr.topLeftCorner(rank, rank)
367  .template triangularView<Upper>()
368  .solveInPlace(c.topRows(rank));
369 
370  dst.topRows(rank) = c.topRows(rank);
371  dst.bottomRows(cols()-rank).setZero();
372 }
373 #endif
374 
381 template<typename MatrixType>
383 {
384  check_template_parameters();
385 
386  Index rows = m_qr.rows();
387  Index cols = m_qr.cols();
388  Index size = (std::min)(rows,cols);
389 
390  m_hCoeffs.resize(size);
391 
392  m_temp.resize(cols);
393 
394  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
395 
396  m_isInitialized = true;
397 }
398 
399 #ifndef __CUDACC__
400 
404 template<typename Derived>
407 {
408  return HouseholderQR<PlainObject>(eval());
409 }
410 #endif // __CUDACC__
411 
412 } // end namespace Eigen
413 
414 #endif // EIGEN_QR_H
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:157
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:78
Namespace containing all symbols from the Eigen library.
Definition: Core:271
Definition: Half.h:502
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:215
Derived & derived()
Definition: EigenBase.h:44
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:273
Definition: EigenBase.h:28
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:115
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Householder QR decomposition of a matrix.
Definition: ForwardDeclarations.h:254
void computeInPlace()
Definition: HouseholderQR.h:382
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: HouseholderQR.h:143
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:248
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
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:239
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:406
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:97
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:166