Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19namespace internal {
20template<typename MatrixType_> struct traits<HouseholderQR<MatrixType_> >
21 : traits<MatrixType_>
22{
23 typedef MatrixXpr XprKind;
24 typedef SolverStorage StorageKind;
25 typedef int StorageIndex;
26 enum { Flags = 0 };
27};
28
29} // end namespace internal
30
58template<typename MatrixType_> class HouseholderQR
59 : public SolverBase<HouseholderQR<MatrixType_> >
60{
61 public:
62
63 typedef MatrixType_ MatrixType;
65 friend class SolverBase<HouseholderQR>;
66
67 EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR)
68 enum {
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71 };
72 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
73 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
76
83 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
84
92 : m_qr(rows, cols),
93 m_hCoeffs((std::min)(rows,cols)),
94 m_temp(cols),
95 m_isInitialized(false) {}
96
109 template<typename InputType>
110 explicit HouseholderQR(const EigenBase<InputType>& matrix)
111 : m_qr(matrix.rows(), matrix.cols()),
112 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
113 m_temp(matrix.cols()),
114 m_isInitialized(false)
115 {
116 compute(matrix.derived());
117 }
118
119
127 template<typename InputType>
129 : m_qr(matrix.derived()),
130 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
131 m_temp(matrix.cols()),
132 m_isInitialized(false)
133 {
135 }
136
137 #ifdef EIGEN_PARSED_BY_DOXYGEN
152 template<typename Rhs>
153 inline const Solve<HouseholderQR, Rhs>
154 solve(const MatrixBase<Rhs>& b) const;
155 #endif
156
166 {
167 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
168 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
169 }
170
174 const MatrixType& matrixQR() const
175 {
176 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
177 return m_qr;
178 }
179
180 template<typename InputType>
181 HouseholderQR& compute(const EigenBase<InputType>& matrix) {
182 m_qr = matrix.derived();
184 return *this;
185 }
186
200 typename MatrixType::RealScalar absDeterminant() const;
201
214 typename MatrixType::RealScalar logAbsDeterminant() const;
215
216 inline Index rows() const { return m_qr.rows(); }
217 inline Index cols() const { return m_qr.cols(); }
218
223 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
224
225 #ifndef EIGEN_PARSED_BY_DOXYGEN
226 template<typename RhsType, typename DstType>
227 void _solve_impl(const RhsType &rhs, DstType &dst) const;
228
229 template<bool Conjugate, typename RhsType, typename DstType>
230 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
231 #endif
232
233 protected:
234
235 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
236
237 void computeInPlace();
238
239 MatrixType m_qr;
240 HCoeffsType m_hCoeffs;
241 RowVectorType m_temp;
242 bool m_isInitialized;
243};
244
245template<typename MatrixType>
246typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
247{
248 using std::abs;
249 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
250 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
251 return abs(m_qr.diagonal().prod());
252}
253
254template<typename MatrixType>
255typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
256{
257 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
258 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
259 return m_qr.diagonal().cwiseAbs().array().log().sum();
260}
261
262namespace internal {
263
265template<typename MatrixQR, typename HCoeffs>
266void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
267{
268 typedef typename MatrixQR::Scalar Scalar;
269 typedef typename MatrixQR::RealScalar RealScalar;
270 Index rows = mat.rows();
271 Index cols = mat.cols();
272 Index size = (std::min)(rows,cols);
273
274 eigen_assert(hCoeffs.size() == size);
275
277 TempType tempVector;
278 if(tempData==0)
279 {
280 tempVector.resize(cols);
281 tempData = tempVector.data();
282 }
283
284 for(Index k = 0; k < size; ++k)
285 {
286 Index remainingRows = rows - k;
287 Index remainingCols = cols - k - 1;
288
289 RealScalar beta;
290 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
291 mat.coeffRef(k,k) = beta;
292
293 // apply H to remaining part of m_qr from the left
294 mat.bottomRightCorner(remainingRows, remainingCols)
295 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
296 }
297}
298
300template<typename MatrixQR, typename HCoeffs,
301 typename MatrixQRScalar = typename MatrixQR::Scalar,
302 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
303struct householder_qr_inplace_blocked
304{
305 // This is specialized for LAPACK-supported Scalar types in HouseholderQR_LAPACKE.h
306 static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
307 typename MatrixQR::Scalar* tempData = 0)
308 {
309 typedef typename MatrixQR::Scalar Scalar;
310 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
311
312 Index rows = mat.rows();
313 Index cols = mat.cols();
314 Index size = (std::min)(rows, cols);
315
316 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
317 TempType tempVector;
318 if(tempData==0)
319 {
320 tempVector.resize(cols);
321 tempData = tempVector.data();
322 }
323
324 Index blockSize = (std::min)(maxBlockSize,size);
325
326 Index k = 0;
327 for (k = 0; k < size; k += blockSize)
328 {
329 Index bs = (std::min)(size-k,blockSize); // actual size of the block
330 Index tcols = cols - k - bs; // trailing columns
331 Index brows = rows-k; // rows of the block
332
333 // partition the matrix:
334 // A00 | A01 | A02
335 // mat = A10 | A11 | A12
336 // A20 | A21 | A22
337 // and performs the qr dec of [A11^T A12^T]^T
338 // and update [A21^T A22^T]^T using level 3 operations.
339 // Finally, the algorithm continue on A22
340
341 BlockType A11_21 = mat.block(k,k,brows,bs);
342 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
343
344 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
345
346 if(tcols)
347 {
348 BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
349 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
350 }
351 }
352 }
353};
354
355} // end namespace internal
356
357#ifndef EIGEN_PARSED_BY_DOXYGEN
358template<typename MatrixType_>
359template<typename RhsType, typename DstType>
360void HouseholderQR<MatrixType_>::_solve_impl(const RhsType &rhs, DstType &dst) const
362 const Index rank = (std::min)(rows(), cols());
363
364 typename RhsType::PlainObject c(rhs);
365
366 c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
367
368 m_qr.topLeftCorner(rank, rank)
369 .template triangularView<Upper>()
370 .solveInPlace(c.topRows(rank));
371
372 dst.topRows(rank) = c.topRows(rank);
373 dst.bottomRows(cols()-rank).setZero();
374}
375
376template<typename MatrixType_>
377template<bool Conjugate, typename RhsType, typename DstType>
378void HouseholderQR<MatrixType_>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
379{
380 const Index rank = (std::min)(rows(), cols());
381
382 typename RhsType::PlainObject c(rhs);
383
384 m_qr.topLeftCorner(rank, rank)
385 .template triangularView<Upper>()
386 .transpose().template conjugateIf<Conjugate>()
387 .solveInPlace(c.topRows(rank));
388
389 dst.topRows(rank) = c.topRows(rank);
390 dst.bottomRows(rows()-rank).setZero();
391
392 dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
393}
394#endif
395
402template<typename MatrixType>
404{
405 Index rows = m_qr.rows();
406 Index cols = m_qr.cols();
407 Index size = (std::min)(rows,cols);
408
409 m_hCoeffs.resize(size);
410
411 m_temp.resize(cols);
412
413 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
414
415 m_isInitialized = true;
416}
417
422template<typename Derived>
425{
426 return HouseholderQR<PlainObject>(eval());
427}
428
429} // end namespace Eigen
430
431#endif // EIGEN_QR_H
Householder QR decomposition of a matrix.
Definition: HouseholderQR.h:60
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:91
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:128
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:223
void computeInPlace()
Definition: HouseholderQR.h:403
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:110
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:165
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:174
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:83
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:246
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:255
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const HouseholderQR< PlainObject > householderQr() const
Definition: HouseholderQR.h:424
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:283
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
HouseholderQR< MatrixType_ > & derived()
Definition: EigenBase.h:48
@ ColMajor
Definition: Constants.h:321
@ RowMajor
Definition: Constants.h:323
const unsigned int RowMajorBit
Definition: Constants.h:68
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Derived & derived()
Definition: EigenBase.h:48