diff options
author | Adolfo Rodriguez Tsourouksdissian <adolfo.rodriguez@pal-robotics.com> | 2011-03-08 19:04:31 +0100 |
---|---|---|
committer | Adolfo Rodriguez Tsourouksdissian <adolfo.rodriguez@pal-robotics.com> | 2011-03-08 19:04:31 +0100 |
commit | 5e431779f38ffa238e552504732d65628884f483 (patch) | |
tree | 99c546c2863010c048a82a02553a9582c1327ca4 /Eigen/src/QR | |
parent | 7bf0e8cd82271c26f1c37b96132c360b7ce7bfbc (diff) |
bug #206 - part 3: Reimplement FullPivHouseholderQR<T>::matrixQ() using ReturnByValue
Diffstat (limited to 'Eigen/src/QR')
-rw-r--r-- | Eigen/src/QR/FullPivHouseholderQR.h | 97 |
1 files changed, 78 insertions, 19 deletions
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index dde3013be..4992c2c3c 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -26,6 +26,18 @@ #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H +namespace internal { + +template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType; + +template<typename MatrixType> +struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > +{ + typedef typename MatrixType::PlainObject ReturnType; +}; + +} + /** \ingroup QR_Module * * \class FullPivHouseholderQR @@ -62,7 +74,7 @@ template<typename _MatrixType> class FullPivHouseholderQR typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; + typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType; typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; @@ -139,7 +151,9 @@ template<typename _MatrixType> class FullPivHouseholderQR return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived()); } - MatrixQType matrixQ(void) const; + /** \returns Expression object representing the matrix Q + */ + MatrixQReturnType matrixQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -508,28 +522,73 @@ struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> } }; +/** \ingroup QR_Module + * + * \brief Expression type for return value of FullPivHouseholderQR::matrixQ() + * + * \tparam MatrixType type of underlying dense matrix + */ +template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType + : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > +{ +public: + typedef typename MatrixType::Index Index; + typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType; + typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; + typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, + MatrixType::MaxRowsAtCompileTime> WorkVectorType; + + FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr, + const HCoeffsType& hCoeffs, + const IntColVectorType& rowsTranspositions) + : m_qr(qr), + m_hCoeffs(hCoeffs), + m_rowsTranspositions(rowsTranspositions) + {} + + template <typename ResultType> + void evalTo(ResultType& result) const + { + const Index rows = m_qr.rows(); + WorkVectorType workspace(rows); + evalTo(result, workspace); + } + + template <typename ResultType> + void evalTo(ResultType& result, WorkVectorType& workspace) const + { + // compute the product H'_0 H'_1 ... H'_n-1, + // where H_k is the k-th Householder transformation I - h_k v_k v_k' + // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] + const Index rows = m_qr.rows(); + const Index cols = m_qr.cols(); + const Index size = std::min(rows, cols); + workspace.resize(rows); + result.setIdentity(rows, rows); + for (Index k = size-1; k >= 0; k--) + { + result.block(k, k, rows-k, rows-k) + .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k)); + result.row(k).swap(result.row(m_rowsTranspositions.coeff(k))); + } + } + + Index rows() const { return m_qr.rows(); } + Index cols() const { return m_qr.rows(); } + +protected: + const typename MatrixType::Nested m_qr; + const typename HCoeffsType::Nested m_hCoeffs; + const typename IntColVectorType::Nested m_rowsTranspositions; +}; + } // end namespace internal -/** \returns the matrix Q */ template<typename MatrixType> -typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const +inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const { eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - // compute the product H'_0 H'_1 ... H'_n-1, - // where H_k is the k-th Householder transformation I - h_k v_k v_k' - // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] - Index rows = m_qr.rows(); - Index cols = m_qr.cols(); - Index size = (std::min)(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows); - for (Index k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - res.row(k).swap(res.row(m_rows_transpositions.coeff(k))); - } - return res; + return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); } /** \return the full-pivoting Householder QR decomposition of \c *this. |