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 
42 template<typename _MatrixType> class HouseholderQR
43 {
44  public:
45 
46  typedef _MatrixType MatrixType;
47  enum {
48  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
49  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
50  Options = MatrixType::Options,
51  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
52  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
53  };
54  typedef typename MatrixType::Scalar Scalar;
55  typedef typename MatrixType::RealScalar RealScalar;
56  typedef typename MatrixType::Index Index;
61 
69 
76  HouseholderQR(Index rows, Index cols)
77  : m_qr(rows, cols),
78  m_hCoeffs((std::min)(rows,cols)),
79  m_temp(cols),
80  m_isInitialized(false) {}
81 
94  HouseholderQR(const MatrixType& matrix)
95  : m_qr(matrix.rows(), matrix.cols()),
96  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
97  m_temp(matrix.cols()),
98  m_isInitialized(false)
99  {
100  compute(matrix);
101  }
102 
120  template<typename Rhs>
122  solve(const MatrixBase<Rhs>& b) const
123  {
124  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
125  return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
126  }
127 
136  HouseholderSequenceType householderQ() const
137  {
138  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
139  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
140  }
141 
145  const MatrixType& matrixQR() const
146  {
147  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
148  return m_qr;
149  }
150 
151  HouseholderQR& compute(const MatrixType& matrix);
152 
166  typename MatrixType::RealScalar absDeterminant() const;
167 
180  typename MatrixType::RealScalar logAbsDeterminant() const;
181 
182  inline Index rows() const { return m_qr.rows(); }
183  inline Index cols() const { return m_qr.cols(); }
184 
189  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
190 
191  protected:
192  MatrixType m_qr;
193  HCoeffsType m_hCoeffs;
194  RowVectorType m_temp;
196 };
197 
198 template<typename MatrixType>
199 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
200 {
201  using std::abs;
202  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
203  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
204  return abs(m_qr.diagonal().prod());
205 }
206 
207 template<typename MatrixType>
208 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
209 {
210  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
211  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
212  return m_qr.diagonal().cwiseAbs().array().log().sum();
213 }
214 
215 namespace internal {
216 
218 template<typename MatrixQR, typename HCoeffs>
219 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
220 {
221  typedef typename MatrixQR::Index Index;
222  typedef typename MatrixQR::Scalar Scalar;
223  typedef typename MatrixQR::RealScalar RealScalar;
224  Index rows = mat.rows();
225  Index cols = mat.cols();
226  Index size = (std::min)(rows,cols);
227 
228  eigen_assert(hCoeffs.size() == size);
229 
231  TempType tempVector;
232  if(tempData==0)
233  {
234  tempVector.resize(cols);
235  tempData = tempVector.data();
236  }
237 
238  for(Index k = 0; k < size; ++k)
239  {
240  Index remainingRows = rows - k;
241  Index remainingCols = cols - k - 1;
242 
243  RealScalar beta;
244  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
245  mat.coeffRef(k,k) = beta;
246 
247  // apply H to remaining part of m_qr from the left
248  mat.bottomRightCorner(remainingRows, remainingCols)
249  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
250  }
251 }
252 
254 template<typename MatrixQR, typename HCoeffs>
255 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
256  typename MatrixQR::Index maxBlockSize=32,
257  typename MatrixQR::Scalar* tempData = 0)
258 {
259  typedef typename MatrixQR::Index Index;
260  typedef typename MatrixQR::Scalar Scalar;
261  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
262 
263  Index rows = mat.rows();
264  Index cols = mat.cols();
265  Index size = (std::min)(rows, cols);
266 
268  TempType tempVector;
269  if(tempData==0)
270  {
271  tempVector.resize(cols);
272  tempData = tempVector.data();
273  }
274 
275  Index blockSize = (std::min)(maxBlockSize,size);
276 
277  Index k = 0;
278  for (k = 0; k < size; k += blockSize)
279  {
280  Index bs = (std::min)(size-k,blockSize); // actual size of the block
281  Index tcols = cols - k - bs; // trailing columns
282  Index brows = rows-k; // rows of the block
283 
284  // partition the matrix:
285  // A00 | A01 | A02
286  // mat = A10 | A11 | A12
287  // A20 | A21 | A22
288  // and performs the qr dec of [A11^T A12^T]^T
289  // and update [A21^T A22^T]^T using level 3 operations.
290  // Finally, the algorithm continue on A22
291 
292  BlockType A11_21 = mat.block(k,k,brows,bs);
293  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
294 
295  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
296 
297  if(tcols)
298  {
299  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
300  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
301  }
302  }
303 }
304 
305 template<typename _MatrixType, typename Rhs>
306 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
307  : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
308 {
310 
311  template<typename Dest> void evalTo(Dest& dst) const
312  {
313  const Index rows = dec().rows(), cols = dec().cols();
314  const Index rank = (std::min)(rows, cols);
315  eigen_assert(rhs().rows() == rows);
316 
317  typename Rhs::PlainObject c(rhs());
318 
319  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
320  c.applyOnTheLeft(householderSequence(
321  dec().matrixQR().leftCols(rank),
322  dec().hCoeffs().head(rank)).transpose()
323  );
324 
325  dec().matrixQR()
326  .topLeftCorner(rank, rank)
327  .template triangularView<Upper>()
328  .solveInPlace(c.topRows(rank));
329 
330  dst.topRows(rank) = c.topRows(rank);
331  dst.bottomRows(cols-rank).setZero();
332  }
333 };
334 
335 } // end namespace internal
336 
343 template<typename MatrixType>
345 {
346  Index rows = matrix.rows();
347  Index cols = matrix.cols();
348  Index size = (std::min)(rows,cols);
349 
350  m_qr = matrix;
351  m_hCoeffs.resize(size);
352 
353  m_temp.resize(cols);
354 
356 
357  m_isInitialized = true;
358  return *this;
359 }
360 
365 template<typename Derived>
368 {
369  return HouseholderQR<PlainObject>(eval());
370 }
371 
372 } // end namespace Eigen
373 
374 #endif // EIGEN_QR_H
HouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:94
HouseholderSequenceType householderQ() const
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:76
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:57
ColsBlockXpr leftCols(Index n)
Definition: BlockMethods.h:514
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
const HCoeffsType & hCoeffs() const
const internal::solve_retval< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Sequence of Householder reflections acting on subspaces with decreasing size.
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
void householder_qr_inplace_blocked(MatrixQR &mat, HCoeffs &hCoeffs, typename MatrixQR::Index maxBlockSize=32, typename MatrixQR::Scalar *tempData=0)
MatrixType::Index Index
Definition: HouseholderQR.h:56
MatrixType::RealScalar RealScalar
Definition: HouseholderQR.h:55
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
Index cols() const
internal::plain_row_type< MatrixType >::type RowVectorType
Definition: HouseholderQR.h:59
void rhs(const real_t *x, real_t *f)
SegmentReturnType head(Index vecSize)
Definition: BlockMethods.h:781
void apply_block_householder_on_the_left(MatrixType &mat, const VectorsType &vectors, const CoeffsType &hCoeffs)
RowVectorType m_temp
_MatrixType MatrixType
Definition: HouseholderQR.h:46
Expression of a fixed-size or dynamic-size block.
Definition: Core/Block.h:102
Householder QR decomposition of a matrix.
internal::plain_diag_type< MatrixType >::type HCoeffsType
Definition: HouseholderQR.h:58
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
HouseholderQR & compute(const MatrixType &matrix)
MatrixType::RealScalar logAbsDeterminant() const
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:68
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
#define eigen_assert(x)
HouseholderSequence< MatrixType, typename internal::remove_all< typename HCoeffsType::ConjugateReturnType >::type > HouseholderSequenceType
Definition: HouseholderQR.h:60
Index rows() const
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
MatrixType::RealScalar absDeterminant() const
const HouseholderQR< PlainObject > householderQr() const
MatrixType::Scalar Scalar
Definition: HouseholderQR.h:54
const MatrixType & matrixQR() const


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:41