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 
17 namespace internal {
18 template<typename _MatrixType> struct traits<HouseholderQR<_MatrixType> >
19  : traits<_MatrixType>
20 {
21  typedef MatrixXpr XprKind;
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 
66  enum {
67  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
68  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
69  };
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>
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
136 
150  template<typename Rhs>
151  inline const Solve<HouseholderQR, Rhs>
152  solve(const MatrixBase<Rhs>& b) const;
153  #endif
154 
163  HouseholderSequenceType householderQ() const
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>
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 
234  {
236  }
237 
238  void computeInPlace();
239 
240  MatrixType m_qr;
241  HCoeffsType m_hCoeffs;
242  RowVectorType m_temp;
244 };
245 
246 template<typename MatrixType>
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>
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)>
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 
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 
417 
418  m_isInitialized = true;
419 }
420 
425 template<typename Derived>
428 {
430 }
431 
432 } // end namespace Eigen
433 
434 #endif // EIGEN_QR_H
HouseholderQR & compute(const EigenBase< InputType > &matrix)
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:1264
static void run(MatrixQR &mat, HCoeffs &hCoeffs, Index maxBlockSize=32, typename MatrixQR::Scalar *tempData=0)
SolverBase< HouseholderQR > Base
Definition: HouseholderQR.h:62
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:89
void apply_block_householder_on_the_left(MatrixType &mat, const VectorsType &vectors, const CoeffsType &hCoeffs, bool forward)
internal::traits< HouseholderQR< TransposeTypeWithSameStorageOrder > >::Scalar Scalar
Definition: SolverBase.h:73
void _solve_impl(const RhsType &rhs, DstType &dst) const
#define min(a, b)
Definition: datatypes.h:19
const MatrixType & matrixQR() const
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:70
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
static void check_template_parameters()
MatrixType::RealScalar absDeterminant() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Definition: BFloat16.h:88
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Sequence of Householder reflections acting on subspaces with decreasing size.
void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Index rows() const
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
#define eigen_assert(x)
Definition: Macros.h:1037
MatrixType::RealScalar logAbsDeterminant() const
const HouseholderQR< PlainObject > householderQr() const
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
HouseholderSequenceType householderQ() const
Index cols() const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
internal::plain_row_type< MatrixType >::type RowVectorType
Definition: HouseholderQR.h:72
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
RowVectorType m_temp
_MatrixType MatrixType
Definition: HouseholderQR.h:61
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:71
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:81
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)
The matrix class, also used for vectors and row-vectors.
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Definition: HouseholderQR.h:73
#define abs(x)
Definition: datatypes.h:17
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:68
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
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.
const HCoeffsType & hCoeffs() const


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:20