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