Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
15namespace Eigen {
16
17namespace internal {
18template<typename _MatrixType> struct traits<HouseholderQR<_MatrixType> >
19 : traits<_MatrixType>
20{
21 typedef MatrixXpr XprKind;
22 typedef SolverStorage StorageKind;
23 typedef int StorageIndex;
24 enum { Flags = 0 };
25};
26
27} // end namespace internal
28
56template<typename _MatrixType> class HouseholderQR
57 : public SolverBase<HouseholderQR<_MatrixType> >
58{
59 public:
60
61 typedef _MatrixType MatrixType;
63 friend class SolverBase<HouseholderQR>;
64
65 EIGEN_GENERIC_PUBLIC_INTERFACE(HouseholderQR)
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;
71 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
72 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
74
81 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
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>
108 explicit HouseholderQR(const EigenBase<InputType>& matrix)
109 : m_qr(matrix.rows(), matrix.cols()),
110 m_hCoeffs((std::min)(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()),
128 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
129 m_temp(matrix.cols()),
130 m_isInitialized(false)
131 {
133 }
134
135 #ifdef EIGEN_PARSED_BY_DOXYGEN
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>
179 HouseholderQR& compute(const EigenBase<InputType>& matrix) {
180 m_qr = matrix.derived();
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
233 static void check_template_parameters()
234 {
235 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
236 }
237
238 void computeInPlace();
239
240 MatrixType m_qr;
241 HCoeffsType m_hCoeffs;
242 RowVectorType m_temp;
243 bool m_isInitialized;
244};
245
246template<typename MatrixType>
247typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
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
255template<typename MatrixType>
256typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
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
263namespace internal {
264
266template<typename MatrixQR, typename HCoeffs>
267void 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
290 RealScalar beta;
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
301template<typename MatrixQR, typename HCoeffs,
302 typename MatrixQRScalar = typename MatrixQR::Scalar,
303 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
304struct householder_qr_inplace_blocked
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
317 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
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
359template<typename _MatrixType>
360template<typename RhsType, typename DstType>
361void 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
377template<typename _MatrixType>
378template<bool Conjugate, typename RhsType, typename DstType>
379void 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
403template<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
416 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
417
418 m_isInitialized = true;
419}
420
425template<typename Derived>
428{
429 return HouseholderQR<PlainObject>(eval());
430}
431
432} // end namespace Eigen
433
434#endif // EIGEN_QR_H
Householder QR decomposition of a matrix.
Definition: HouseholderQR.h:58
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:221
HouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: HouseholderQR.h:89
void computeInPlace()
Definition: HouseholderQR.h:404
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:172
const Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
HouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:126
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:81
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:247
MatrixType::RealScalar logAbsDeterminant() const
Definition: HouseholderQR.h:256
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:108
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:163
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Pseudo expression representing a solving operation.
Definition: Solve.h:63
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
HouseholderQR< _MatrixType > & derived()
Definition: EigenBase.h:46
@ ColMajor
Definition: Constants.h:319
@ RowMajor
Definition: Constants.h:321
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Derived & derived()
Definition: EigenBase.h:46