Eigen  3.2.91
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 
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  // FIXME should be int
57  typedef typename MatrixType::StorageIndex StorageIndex;
58  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
59  typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
60  typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
61  typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
62 
69  HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
70 
77  HouseholderQR(Index rows, Index cols)
78  : m_qr(rows, cols),
79  m_hCoeffs((std::min)(rows,cols)),
80  m_temp(cols),
81  m_isInitialized(false) {}
82 
95  explicit HouseholderQR(const MatrixType& matrix)
96  : m_qr(matrix.rows(), matrix.cols()),
97  m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
98  m_temp(matrix.cols()),
99  m_isInitialized(false)
100  {
101  compute(matrix);
102  }
103 
121  template<typename Rhs>
122  inline const Solve<HouseholderQR, Rhs>
123  solve(const MatrixBase<Rhs>& b) const
124  {
125  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
126  return Solve<HouseholderQR, Rhs>(*this, b.derived());
127  }
128 
137  HouseholderSequenceType householderQ() const
138  {
139  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
140  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
141  }
142 
146  const MatrixType& matrixQR() const
147  {
148  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
149  return m_qr;
150  }
151 
152  HouseholderQR& compute(const MatrixType& matrix);
153 
167  typename MatrixType::RealScalar absDeterminant() const;
168 
181  typename MatrixType::RealScalar logAbsDeterminant() const;
182 
183  inline Index rows() const { return m_qr.rows(); }
184  inline Index cols() const { return m_qr.cols(); }
185 
190  const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
191 
192  #ifndef EIGEN_PARSED_BY_DOXYGEN
193  template<typename RhsType, typename DstType>
194  EIGEN_DEVICE_FUNC
195  void _solve_impl(const RhsType &rhs, DstType &dst) const;
196  #endif
197 
198  protected:
199 
200  static void check_template_parameters()
201  {
202  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
203  }
204 
205  MatrixType m_qr;
206  HCoeffsType m_hCoeffs;
207  RowVectorType m_temp;
208  bool m_isInitialized;
209 };
210 
211 template<typename MatrixType>
212 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
213 {
214  using std::abs;
215  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
216  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
217  return abs(m_qr.diagonal().prod());
218 }
219 
220 template<typename MatrixType>
221 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
222 {
223  eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
224  eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
225  return m_qr.diagonal().cwiseAbs().array().log().sum();
226 }
227 
228 namespace internal {
229 
231 template<typename MatrixQR, typename HCoeffs>
232 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
233 {
234  typedef typename MatrixQR::Scalar Scalar;
235  typedef typename MatrixQR::RealScalar RealScalar;
236  Index rows = mat.rows();
237  Index cols = mat.cols();
238  Index size = (std::min)(rows,cols);
239 
240  eigen_assert(hCoeffs.size() == size);
241 
243  TempType tempVector;
244  if(tempData==0)
245  {
246  tempVector.resize(cols);
247  tempData = tempVector.data();
248  }
249 
250  for(Index k = 0; k < size; ++k)
251  {
252  Index remainingRows = rows - k;
253  Index remainingCols = cols - k - 1;
254 
255  RealScalar beta;
256  mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
257  mat.coeffRef(k,k) = beta;
258 
259  // apply H to remaining part of m_qr from the left
260  mat.bottomRightCorner(remainingRows, remainingCols)
261  .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
262  }
263 }
264 
266 template<typename MatrixQR, typename HCoeffs,
267  typename MatrixQRScalar = typename MatrixQR::Scalar,
268  bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
269 struct householder_qr_inplace_blocked
270 {
271  // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
272  static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
273  typename MatrixQR::Scalar* tempData = 0)
274  {
275  typedef typename MatrixQR::Scalar Scalar;
276  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
277 
278  Index rows = mat.rows();
279  Index cols = mat.cols();
280  Index size = (std::min)(rows, cols);
281 
282  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
283  TempType tempVector;
284  if(tempData==0)
285  {
286  tempVector.resize(cols);
287  tempData = tempVector.data();
288  }
289 
290  Index blockSize = (std::min)(maxBlockSize,size);
291 
292  Index k = 0;
293  for (k = 0; k < size; k += blockSize)
294  {
295  Index bs = (std::min)(size-k,blockSize); // actual size of the block
296  Index tcols = cols - k - bs; // trailing columns
297  Index brows = rows-k; // rows of the block
298 
299  // partition the matrix:
300  // A00 | A01 | A02
301  // mat = A10 | A11 | A12
302  // A20 | A21 | A22
303  // and performs the qr dec of [A11^T A12^T]^T
304  // and update [A21^T A22^T]^T using level 3 operations.
305  // Finally, the algorithm continue on A22
306 
307  BlockType A11_21 = mat.block(k,k,brows,bs);
308  Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
309 
310  householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
311 
312  if(tcols)
313  {
314  BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
315  apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
316  }
317  }
318  }
319 };
320 
321 } // end namespace internal
322 
323 #ifndef EIGEN_PARSED_BY_DOXYGEN
324 template<typename _MatrixType>
325 template<typename RhsType, typename DstType>
326 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
327 {
328  const Index rank = (std::min)(rows(), cols());
329  eigen_assert(rhs.rows() == rows());
330 
331  typename RhsType::PlainObject c(rhs);
332 
333  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
334  c.applyOnTheLeft(householderSequence(
335  m_qr.leftCols(rank),
336  m_hCoeffs.head(rank)).transpose()
337  );
338 
339  m_qr.topLeftCorner(rank, rank)
340  .template triangularView<Upper>()
341  .solveInPlace(c.topRows(rank));
342 
343  dst.topRows(rank) = c.topRows(rank);
344  dst.bottomRows(cols()-rank).setZero();
345 }
346 #endif
347 
354 template<typename MatrixType>
356 {
357  check_template_parameters();
358 
359  Index rows = matrix.rows();
360  Index cols = matrix.cols();
361  Index size = (std::min)(rows,cols);
362 
363  m_qr = matrix;
364  m_hCoeffs.resize(size);
365 
366  m_temp.resize(cols);
367 
368  internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
369 
370  m_isInitialized = true;
371  return *this;
372 }
373 
374 #ifndef __CUDACC__
375 
379 template<typename Derived>
382 {
383  return HouseholderQR<PlainObject>(eval());
384 }
385 #endif // __CUDACC__
386 
387 } // end namespace Eigen
388 
389 #endif // EIGEN_QR_H
HouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:95
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:137
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:77
Definition: LDLT.h:16
Definition: StdDeque.h:58
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:190
HouseholderSequence< VectorsType, CoeffsType > householderSequence(const VectorsType &v, const CoeffsType &h)
Convenience function for constructing a Householder sequence.
Definition: HouseholderSequence.h:452
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Definition: Eigen_Colamd.h:54
Householder QR decomposition of a matrix.
Definition: ForwardDeclarations.h:251
HouseholderQR & compute(const MatrixType &matrix)
Definition: HouseholderQR.h:355
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: HouseholderQR.h:123
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:221
Pseudo expression representing a solving operation.
Definition: Solve.h:63
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:69
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:212
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:381
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:146