HouseholderQR.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) 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;
57  // FIXME should be int
58  typedef typename MatrixType::StorageIndex StorageIndex;
63 
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>
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 
154  HouseholderSequenceType householderQ() const
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>
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 
223  {
225  }
226 
227  void computeInPlace();
228 
229  MatrixType m_qr;
230  HCoeffsType m_hCoeffs;
231  RowVectorType m_temp;
233 };
234 
235 template<typename MatrixType>
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>
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)>
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 
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 {
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 
392 
393  m_isInitialized = true;
394 }
395 
400 template<typename Derived>
403 {
405 }
406 
407 } // end namespace Eigen
408 
409 #endif // EIGEN_QR_H
HouseholderQR & compute(const EigenBase< InputType > &matrix)
SCALAR Scalar
Definition: bench_gemm.cpp:33
HouseholderSequenceType householderQ() const
static void run(MatrixQR &mat, HCoeffs &hCoeffs, Index maxBlockSize=32, typename MatrixQR::Scalar *tempData=0)
MatrixType::StorageIndex StorageIndex
Definition: HouseholderQR.h:58
Scalar * b
Definition: benchVecAdd.cpp:17
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:78
void apply_block_householder_on_the_left(MatrixType &mat, const VectorsType &vectors, const CoeffsType &hCoeffs, bool forward)
#define min(a, b)
Definition: datatypes.h:19
void householder_qr_inplace_unblocked(MatrixQR &mat, HCoeffs &hCoeffs, typename MatrixQR::Scalar *tempData=0)
Matrix< Scalar, RowsAtCompileTime, RowsAtCompileTime,(MatrixType::Flags &RowMajorBit)?RowMajor:ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime > MatrixQType
Definition: HouseholderQR.h:59
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
static void check_template_parameters()
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Definition: Half.h:150
const HCoeffsType & hCoeffs() const
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Sequence of Householder reflections acting on subspaces with decreasing size.
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
MatrixType::RealScalar RealScalar
Definition: HouseholderQR.h:56
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.
#define eigen_assert(x)
Definition: Macros.h:579
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:184
Index cols() const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
internal::plain_row_type< MatrixType >::type RowVectorType
Definition: HouseholderQR.h:61
RowVectorType m_temp
_MatrixType MatrixType
Definition: HouseholderQR.h:48
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Householder QR decomposition of a matrix.
internal::plain_diag_type< MatrixType >::type HCoeffsType
Definition: HouseholderQR.h:60
internal::nested_eval< T, 1 >::type eval(const T &xpr)
EIGEN_DEVICE_FUNC void _solve_impl(const RhsType &rhs, DstType &dst) const
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
MatrixType::RealScalar logAbsDeterminant() const
Pseudo expression representing a solving operation.
Definition: Solve.h:62
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:70
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
The matrix class, also used for vectors and row-vectors.
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Definition: HouseholderQR.h:62
Index rows() const
#define abs(x)
Definition: datatypes.h:17
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
MatrixType::RealScalar absDeterminant() const
virtual EIGEN_DEVICE_FUNC const Scalar * data() const
Definition: TensorRef.h:59
const HouseholderQR< PlainObject > householderQr() const
MatrixType::Scalar Scalar
Definition: HouseholderQR.h:55
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:97
const MatrixType & matrixQR() const


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:12