15 #include "./InternalHeaderCheck.h"
20 template<
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;
245 template<
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());
254 template<
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();
265 template<
typename MatrixQR,
typename HCoeffs>
266 void 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);
308 template <
typename MatrixQR,
typename HCoeffs,
typename VectorQR>
309 void householder_qr_inplace_update(MatrixQR& mat, HCoeffs& hCoeffs,
const VectorQR& newColumn,
310 typename MatrixQR::Index k,
typename MatrixQR::Scalar* tempData) {
311 typedef typename MatrixQR::Index
Index;
312 typedef typename MatrixQR::RealScalar RealScalar;
313 Index rows = mat.rows();
315 eigen_assert(k < mat.cols());
316 eigen_assert(k < rows);
317 eigen_assert(hCoeffs.size() == mat.cols());
318 eigen_assert(newColumn.size() == rows);
319 eigen_assert(tempData);
322 mat.col(k) = newColumn;
324 for (
Index i = 0; i < k; ++i) {
325 Index remainingRows = rows - i;
328 .applyHouseholderOnTheLeft(mat.col(i).tail(remainingRows - 1), hCoeffs.coeffRef(i), tempData + i + 1);
332 mat.col(k).tail(rows - k).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
333 mat.coeffRef(k, k) = beta;
337 template<
typename MatrixQR,
typename HCoeffs,
338 typename MatrixQRScalar =
typename MatrixQR::Scalar,
339 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
340 struct householder_qr_inplace_blocked
343 static void run(MatrixQR& mat, HCoeffs& hCoeffs,
Index maxBlockSize=32,
344 typename MatrixQR::Scalar* tempData = 0)
346 typedef typename MatrixQR::Scalar Scalar;
347 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
349 Index rows = mat.rows();
350 Index cols = mat.cols();
351 Index size = (std::min)(rows, cols);
353 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
357 tempVector.resize(cols);
358 tempData = tempVector.data();
361 Index blockSize = (std::min)(maxBlockSize,size);
364 for (k = 0; k < size; k += blockSize)
366 Index bs = (std::min)(size-k,blockSize);
367 Index tcols = cols - k - bs;
368 Index brows = rows-k;
378 BlockType A11_21 = mat.block(k,k,brows,bs);
381 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
385 BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
386 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment,
false);
394 #ifndef EIGEN_PARSED_BY_DOXYGEN
395 template<
typename MatrixType_>
396 template<
typename RhsType,
typename DstType>
397 void HouseholderQR<MatrixType_>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
399 const Index rank = (std::min)(rows(), cols());
401 typename RhsType::PlainObject c(rhs);
403 c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
405 m_qr.topLeftCorner(rank, rank)
406 .template triangularView<Upper>()
407 .solveInPlace(c.topRows(rank));
409 dst.topRows(rank) = c.topRows(rank);
410 dst.bottomRows(cols()-rank).setZero();
413 template<
typename MatrixType_>
414 template<
bool Conjugate,
typename RhsType,
typename DstType>
415 void HouseholderQR<MatrixType_>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
417 const Index rank = (std::min)(rows(), cols());
419 typename RhsType::PlainObject c(rhs);
421 m_qr.topLeftCorner(rank, rank)
422 .template triangularView<Upper>()
423 .transpose().template conjugateIf<Conjugate>()
424 .solveInPlace(c.topRows(rank));
426 dst.topRows(rank) = c.topRows(rank);
427 dst.bottomRows(rows()-rank).setZero();
429 dst.applyOnTheLeft(householderQ().setLength(rank).
template conjugateIf<!Conjugate>() );
439 template<
typename MatrixType>
442 Index rows = m_qr.rows();
443 Index cols = m_qr.cols();
444 Index size = (std::min)(rows,cols);
446 m_hCoeffs.resize(size);
450 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
452 m_isInitialized =
true;
459 template<
typename Derived>
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:107
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 Solve< HouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
void computeInPlace()
Definition: HouseholderQR.h:440
HouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: HouseholderQR.h:110
const HCoeffsType & hCoeffs() const
Definition: HouseholderQR.h:223
HouseholderSequenceType householderQ() const
Definition: HouseholderQR.h:165
HouseholderQR()
Default Constructor.
Definition: HouseholderQR.h:83
MatrixType::RealScalar absDeterminant() const
Definition: HouseholderQR.h:246
const MatrixType & matrixQR() const
Definition: HouseholderQR.h:174
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:461
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: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41