15#include "./InternalHeaderCheck.h"
20template<
typename MatrixType_>
struct traits<HouseholderQR<MatrixType_> >
23 typedef MatrixXpr XprKind;
24 typedef SolverStorage StorageKind;
25 typedef int StorageIndex;
59 :
public SolverBase<HouseholderQR<MatrixType_> >
63 typedef MatrixType_ MatrixType;
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
73 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
83 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
93 m_hCoeffs((std::min)(rows,cols)),
95 m_isInitialized(false) {}
109 template<
typename InputType>
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)
127 template<
typename InputType>
130 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
131 m_temp(matrix.cols()),
132 m_isInitialized(false)
137 #ifdef EIGEN_PARSED_BY_DOXYGEN
152 template<
typename Rhs>
167 eigen_assert(m_isInitialized &&
"HouseholderQR is not initialized.");
176 eigen_assert(m_isInitialized &&
"HouseholderQR is not initialized.");
180 template<
typename InputType>
216 inline Index rows()
const {
return m_qr.rows(); }
217 inline Index cols()
const {
return m_qr.cols(); }
223 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
225 #ifndef EIGEN_PARSED_BY_DOXYGEN
226 template<
typename RhsType,
typename DstType>
227 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
229 template<
bool Conjugate,
typename RhsType,
typename DstType>
230 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
235 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
240 HCoeffsType m_hCoeffs;
241 RowVectorType m_temp;
242 bool m_isInitialized;
245template<
typename MatrixType>
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());
254template<
typename MatrixType>
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();
265template<
typename MatrixQR,
typename HCoeffs>
266void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs,
typename MatrixQR::Scalar* tempData = 0)
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);
274 eigen_assert(hCoeffs.size() == size);
281 tempData = tempVector.data();
284 for(
Index k = 0; k < size; ++k)
286 Index remainingRows = rows - k;
287 Index remainingCols = cols - k - 1;
290 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
291 mat.coeffRef(k,k) = beta;
294 mat.bottomRightCorner(remainingRows, remainingCols)
295 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
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
306 static void run(MatrixQR& mat, HCoeffs& hCoeffs,
Index maxBlockSize=32,
307 typename MatrixQR::Scalar* tempData = 0)
309 typedef typename MatrixQR::Scalar Scalar;
310 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
312 Index rows = mat.rows();
313 Index cols = mat.cols();
314 Index size = (std::min)(rows, cols);
316 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
320 tempVector.resize(cols);
321 tempData = tempVector.data();
324 Index blockSize = (std::min)(maxBlockSize,size);
327 for (k = 0; k < size; k += blockSize)
329 Index bs = (std::min)(size-k,blockSize);
330 Index tcols = cols - k - bs;
331 Index brows = rows-k;
341 BlockType A11_21 = mat.block(k,k,brows,bs);
342 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
344 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
348 BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
349 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment,
false);
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());
364 typename RhsType::PlainObject c(rhs);
366 c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
368 m_qr.topLeftCorner(rank, rank)
369 .template triangularView<Upper>()
370 .solveInPlace(c.topRows(rank));
372 dst.topRows(rank) = c.topRows(rank);
373 dst.bottomRows(cols()-rank).setZero();
376template<
typename MatrixType_>
377template<
bool Conjugate,
typename RhsType,
typename DstType>
380 const Index rank = (std::min)(rows(), cols());
382 typename RhsType::PlainObject c(rhs);
384 m_qr.topLeftCorner(rank, rank)
385 .template triangularView<Upper>()
386 .transpose().template conjugateIf<Conjugate>()
387 .solveInPlace(c.topRows(rank));
389 dst.topRows(rank) = c.topRows(rank);
390 dst.bottomRows(rows()-rank).setZero();
392 dst.applyOnTheLeft(householderQ().setLength(rank).
template conjugateIf<!Conjugate>() );
402template<
typename MatrixType>
405 Index rows = m_qr.rows();
406 Index cols = m_qr.cols();
407 Index size = (std::min)(rows,cols);
409 m_hCoeffs.resize(size);
413 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
415 m_isInitialized =
true;
422template<
typename Derived>
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