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  };
70  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
74 
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()),
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()),
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 
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 
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 
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
Eigen::MatrixXpr
Definition: Constants.h:522
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::Block
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
Eigen::HouseholderQR::HouseholderSequenceType
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Definition: HouseholderQR.h:73
Eigen::HouseholderQR::computeInPlace
void computeInPlace()
Definition: HouseholderQR.h:404
Eigen::HouseholderQR::m_isInitialized
bool m_isInitialized
Definition: HouseholderQR.h:243
Eigen::HouseholderQR::HCoeffsType
internal::plain_diag_type< MatrixType >::type HCoeffsType
Definition: HouseholderQR.h:71
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::EigenBase
Definition: EigenBase.h:29
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:66
Eigen::HouseholderQR::HouseholderQR
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:81
Eigen::HouseholderQR::absDeterminant
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:247
Eigen::HouseholderQR::m_hCoeffs
HCoeffsType m_hCoeffs
Definition: HouseholderQR.h:241
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Eigen::internal::traits< HouseholderQR< _MatrixType > >::StorageKind
SolverStorage StorageKind
Definition: HouseholderQR.h:22
beta
double beta(double a, double b)
Definition: beta.c:61
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Eigen::internal::true_type
Definition: Meta.h:96
Eigen::internal::householder_qr_inplace_unblocked
void householder_qr_inplace_unblocked(MatrixQR &mat, HCoeffs &hCoeffs, typename MatrixQR::Scalar *tempData=0)
Definition: HouseholderQR.h:267
Eigen::internal::traits< HouseholderQR< _MatrixType > >::StorageIndex
int StorageIndex
Definition: HouseholderQR.h:23
Eigen::HouseholderQR::householderQ
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:163
EIGEN_GENERIC_PUBLIC_INTERFACE
#define EIGEN_GENERIC_PUBLIC_INTERFACE(Derived)
Definition: Macros.h:1264
Eigen::HouseholderQR::check_template_parameters
static void check_template_parameters()
Definition: HouseholderQR.h:233
Eigen::SolverStorage
Definition: Constants.h:513
Eigen::HouseholderQR::matrixQR
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:172
Eigen::internal::householder_qr_inplace_blocked
Definition: HouseholderQR.h:304
Eigen::HouseholderQR::RowVectorType
internal::plain_row_type< MatrixType >::type RowVectorType
Definition: HouseholderQR.h:72
gtsam.examples.DogLegOptimizerExample.run
def run(args)
Definition: DogLegOptimizerExample.py:21
Eigen::HouseholderQR::rows
Index rows() const
Definition: HouseholderQR.h:214
Eigen::HouseholderQR::HouseholderQR
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:126
Eigen::HouseholderQR::hCoeffs
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:221
Eigen::HouseholderQR::_solve_impl
void _solve_impl(const RhsType &rhs, DstType &dst) const
Definition: HouseholderQR.h:361
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::internal::apply_block_householder_on_the_left
void apply_block_householder_on_the_left(MatrixType &mat, const VectorsType &vectors, const CoeffsType &hCoeffs, bool forward)
Definition: BlockHouseholder.h:86
Eigen::Solve
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Eigen::HouseholderQR::HouseholderQR
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:108
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::HouseholderQR::HouseholderQR
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:89
Eigen::internal::traits< HouseholderQR< _MatrixType > >::XprKind
MatrixXpr XprKind
Definition: HouseholderQR.h:21
Eigen::internal::householder_qr_inplace_blocked::run
static void run(MatrixQR &mat, HCoeffs &hCoeffs, Index maxBlockSize=32, typename MatrixQR::Scalar *tempData=0)
Definition: HouseholderQR.h:307
Eigen::HouseholderQR::MaxRowsAtCompileTime
@ MaxRowsAtCompileTime
Definition: HouseholderQR.h:67
std
Definition: BFloat16.h:88
Eigen::HouseholderQR::cols
Index cols() const
Definition: HouseholderQR.h:215
Eigen::HouseholderQR::Base
SolverBase< HouseholderQR > Base
Definition: HouseholderQR.h:62
Eigen::HouseholderQR::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: HouseholderQR.h:68
Eigen::HouseholderQR::logAbsDeterminant
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:256
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
abs
#define abs(x)
Definition: datatypes.h:17
Eigen::HouseholderQR::compute
HouseholderQR & compute(const EigenBase< InputType > &matrix)
Definition: HouseholderQR.h:179
internal
Definition: BandTriangularSolver.h:13
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::min
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)
Definition: AutoDiffScalar.h:580
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:319
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::HouseholderQR::m_qr
MatrixType m_qr
Definition: HouseholderQR.h:240
Eigen::HouseholderQR::MatrixType
_MatrixType MatrixType
Definition: HouseholderQR.h:61
Eigen::HouseholderQR::m_temp
RowVectorType m_temp
Definition: HouseholderQR.h:242
Eigen::HouseholderQR::_solve_impl_transposed
void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const
Definition: HouseholderQR.h:379
adjoint
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
Eigen::SolverBase
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:68
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::HouseholderQR
Householder QR decomposition of a matrix.
Definition: ForwardDeclarations.h:273
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Eigen::HouseholderSequence
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: ForwardDeclarations.h:282
EIGEN_STATIC_ASSERT_NON_INTEGER
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:43