aboutsummaryrefslogtreecommitdiff
path: root/Eigen/src/QR/HouseholderQR.h
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/QR/HouseholderQR.h')
-rw-r--r--Eigen/src/QR/HouseholderQR.h119
1 files changed, 70 insertions, 49 deletions
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 343a66499..3513d995c 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -21,7 +21,7 @@ namespace Eigen {
*
* \brief Householder QR decomposition of a matrix
*
- * \param MatrixType the type of the matrix of which we are computing the QR decomposition
+ * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
*
* This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R
* such that
@@ -37,6 +37,8 @@ namespace Eigen {
* This Householder QR decomposition is faster, but less numerically stable and less feature-full than
* FullPivHouseholderQR or ColPivHouseholderQR.
*
+ * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
+ *
* \sa MatrixBase::householderQr()
*/
template<typename _MatrixType> class HouseholderQR
@@ -47,13 +49,13 @@ template<typename _MatrixType> class HouseholderQR
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
+ // FIXME should be int
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
@@ -91,13 +93,32 @@ template<typename _MatrixType> class HouseholderQR
*
* \sa compute()
*/
- HouseholderQR(const MatrixType& matrix)
+ template<typename InputType>
+ explicit HouseholderQR(const EigenBase<InputType>& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
m_temp(matrix.cols()),
m_isInitialized(false)
{
- compute(matrix);
+ compute(matrix.derived());
+ }
+
+
+ /** \brief Constructs a QR factorization from a given matrix
+ *
+ * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when
+ * \c MatrixType is a Eigen::Ref.
+ *
+ * \sa HouseholderQR(const EigenBase&)
+ */
+ template<typename InputType>
+ explicit HouseholderQR(EigenBase<InputType>& matrix)
+ : m_qr(matrix.derived()),
+ m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
+ m_temp(matrix.cols()),
+ m_isInitialized(false)
+ {
+ computeInPlace();
}
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which
@@ -107,9 +128,6 @@ template<typename _MatrixType> class HouseholderQR
*
* \returns a solution.
*
- * \note The case where b is a matrix is not yet implemented. Also, this
- * code is space inefficient.
- *
* \note_about_checking_solutions
*
* \note_about_arbitrary_choice_of_solution
@@ -118,11 +136,11 @@ template<typename _MatrixType> class HouseholderQR
* Output: \verbinclude HouseholderQR_solve.out
*/
template<typename Rhs>
- inline const internal::solve_retval<HouseholderQR, Rhs>
+ inline const Solve<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
- return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
+ return Solve<HouseholderQR, Rhs>(*this, b.derived());
}
/** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations.
@@ -148,7 +166,12 @@ template<typename _MatrixType> class HouseholderQR
return m_qr;
}
- HouseholderQR& compute(const MatrixType& matrix);
+ template<typename InputType>
+ HouseholderQR& compute(const EigenBase<InputType>& matrix) {
+ m_qr = matrix.derived();
+ computeInPlace();
+ return *this;
+ }
/** \returns the absolute value of the determinant of the matrix of which
* *this is the QR decomposition. It has only linear complexity
@@ -187,6 +210,12 @@ template<typename _MatrixType> class HouseholderQR
* For advanced uses only.
*/
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
+ #ifndef EIGEN_PARSED_BY_DOXYGEN
+ template<typename RhsType, typename DstType>
+ EIGEN_DEVICE_FUNC
+ void _solve_impl(const RhsType &rhs, DstType &dst) const;
+ #endif
protected:
@@ -194,6 +223,8 @@ template<typename _MatrixType> class HouseholderQR
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
+
+ void computeInPlace();
MatrixType m_qr;
HCoeffsType m_hCoeffs;
@@ -224,7 +255,6 @@ namespace internal {
template<typename MatrixQR, typename HCoeffs>
void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
{
- typedef typename MatrixQR::Index Index;
typedef typename MatrixQR::Scalar Scalar;
typedef typename MatrixQR::RealScalar RealScalar;
Index rows = mat.rows();
@@ -263,11 +293,9 @@ template<typename MatrixQR, typename HCoeffs,
struct householder_qr_inplace_blocked
{
// This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
- static void run(MatrixQR& mat, HCoeffs& hCoeffs,
- typename MatrixQR::Index maxBlockSize=32,
+ static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
typename MatrixQR::Scalar* tempData = 0)
{
- typedef typename MatrixQR::Index Index;
typedef typename MatrixQR::Scalar Scalar;
typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
@@ -289,8 +317,8 @@ struct householder_qr_inplace_blocked
for (k = 0; k < size; k += blockSize)
{
Index bs = (std::min)(size-k,blockSize); // actual size of the block
- Index tcols = cols - k - bs; // trailing columns
- Index brows = rows-k; // rows of the block
+ Index tcols = cols - k - bs; // trailing columns
+ Index brows = rows-k; // rows of the block
// partition the matrix:
// A00 | A01 | A02
@@ -308,43 +336,38 @@ struct householder_qr_inplace_blocked
if(tcols)
{
BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
- apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
+ apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
}
}
}
};
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
- : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
-{
- EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- const Index rows = dec().rows(), cols = dec().cols();
- const Index rank = (std::min)(rows, cols);
- eigen_assert(rhs().rows() == rows);
+} // end namespace internal
- typename Rhs::PlainObject c(rhs());
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename _MatrixType>
+template<typename RhsType, typename DstType>
+void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
+{
+ const Index rank = (std::min)(rows(), cols());
+ eigen_assert(rhs.rows() == rows());
- // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
- c.applyOnTheLeft(householderSequence(
- dec().matrixQR().leftCols(rank),
- dec().hCoeffs().head(rank)).transpose()
- );
+ typename RhsType::PlainObject c(rhs);
- dec().matrixQR()
- .topLeftCorner(rank, rank)
- .template triangularView<Upper>()
- .solveInPlace(c.topRows(rank));
+ // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+ c.applyOnTheLeft(householderSequence(
+ m_qr.leftCols(rank),
+ m_hCoeffs.head(rank)).transpose()
+ );
- dst.topRows(rank) = c.topRows(rank);
- dst.bottomRows(cols-rank).setZero();
- }
-};
+ m_qr.topLeftCorner(rank, rank)
+ .template triangularView<Upper>()
+ .solveInPlace(c.topRows(rank));
-} // end namespace internal
+ dst.topRows(rank) = c.topRows(rank);
+ dst.bottomRows(cols()-rank).setZero();
+}
+#endif
/** Performs the QR factorization of the given matrix \a matrix. The result of
* the factorization is stored into \c *this, and a reference to \c *this
@@ -353,15 +376,14 @@ struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
* \sa class HouseholderQR, HouseholderQR(const MatrixType&)
*/
template<typename MatrixType>
-HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
+void HouseholderQR<MatrixType>::computeInPlace()
{
check_template_parameters();
- Index rows = matrix.rows();
- Index cols = matrix.cols();
+ Index rows = m_qr.rows();
+ Index cols = m_qr.cols();
Index size = (std::min)(rows,cols);
- m_qr = matrix;
m_hCoeffs.resize(size);
m_temp.resize(cols);
@@ -369,7 +391,6 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
m_isInitialized = true;
- return *this;
}
/** \return the Householder QR decomposition of \c *this.
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-08 20:07:21 +0000