From 3e73f6036c4f28ae8d11ae43641c213e608529e6 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Wed, 2 Dec 2009 11:11:09 -0500 Subject: * HouseholderSequence: * be aware of number of actual householder vectors (optimization in non-full-rank case, no behavior change) * fix applyThisOnTheRight, it was using k instead of actual_k * QR: rename matrixQ() to householderQ() where applicable --- Eigen/src/Householder/HouseholderSequence.h | 37 +++++++++++++++++++---------- Eigen/src/QR/ColPivHouseholderQR.h | 10 ++++---- Eigen/src/QR/HouseholderQR.h | 13 ++-------- 3 files changed, 32 insertions(+), 28 deletions(-) (limited to 'Eigen') diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index afeb758b3..1159f03ed 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -69,31 +69,35 @@ template class HouseholderSequence typedef HouseholderSequence::IsComplex, - typename ei_cleantype::type, + NestByValue::type >, CoeffsType>::ret> ConjugateReturnType; HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) - : m_vectors(v), m_coeffs(h), m_trans(trans) + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()) + {} + + HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors) {} int rows() const { return m_vectors.rows(); } int cols() const { return m_vectors.rows(); } HouseholderSequence transpose() const - { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); } + { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); } ConjugateReturnType conjugate() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); } ConjugateReturnType adjoint() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); } ConjugateReturnType inverse() const { return adjoint(); } /** \internal */ template void evalTo(DestType& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); + int vecs = m_actualVectors; int length = m_vectors.rows(); dst.setIdentity(); Matrix temp(dst.rows()); @@ -111,22 +115,22 @@ template class HouseholderSequence /** \internal */ template inline void applyThisOnTheRight(Dest& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector - Matrix temp(dst.rows()); + int vecs = m_actualVectors; // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix temp(dst.rows()); for(int k = 0; k < vecs; ++k) { int actual_k = m_trans ? vecs-k-1 : k; - dst.corner(BottomRight, dst.rows(), length-k) - .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + dst.corner(BottomRight, dst.rows(), length-actual_k) + .applyHouseholderOnTheRight(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); } } /** \internal */ template inline void applyThisOnTheLeft(Dest& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector + int vecs = m_actualVectors; // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector Matrix temp(dst.cols()); for(int k = 0; k < vecs; ++k) { @@ -156,6 +160,7 @@ template class HouseholderSequence typename VectorsType::Nested m_vectors; typename CoeffsType::Nested m_coeffs; bool m_trans; + int m_actualVectors; private: HouseholderSequence& operator=(const HouseholderSequence&); @@ -167,4 +172,10 @@ HouseholderSequence householderSequence(const VectorsTyp return HouseholderSequence(v, h, trans); } +template +HouseholderSequence householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) +{ + return HouseholderSequence(v, h, trans, actualVectors); +} + #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 9c1e56dc6..1b2472a99 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -105,7 +105,7 @@ template class ColPivHouseholderQR return ei_solve_retval(*this, b.derived()); } - HouseholderSequenceType matrixQ(void) const; + HouseholderSequenceType householderQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -447,7 +447,8 @@ struct ei_solve_retval, Rhs> c.applyOnTheLeft(householderSequence( dec().matrixQR(), dec().hCoeffs(), - true + true, + dec().nonzeroPivots() )); dec().matrixQR() @@ -470,10 +471,11 @@ struct ei_solve_retval, Rhs> /** \returns the matrix Q as a sequence of householder transformations */ template -typename ColPivHouseholderQR::HouseholderSequenceType ColPivHouseholderQR::matrixQ() const +typename ColPivHouseholderQR::HouseholderSequenceType ColPivHouseholderQR + ::householderQ() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false); + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 8d842d129..95496b943 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -104,11 +104,10 @@ template class HouseholderQR ei_assert(m_isInitialized && "HouseholderQR is not initialized."); return ei_solve_retval(*this, b.derived()); } - - MatrixQType matrixQ() const; - HouseholderSequenceType matrixQAsHouseholderSequence() const + HouseholderSequenceType householderQ() const { + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } @@ -240,14 +239,6 @@ struct ei_solve_retval, Rhs> } }; -/** \returns the matrix Q */ -template -typename HouseholderQR::MatrixQType HouseholderQR::matrixQ() const -{ - ei_assert(m_isInitialized && "HouseholderQR is not initialized."); - return matrixQAsHouseholderSequence(); -} - #endif // EIGEN_HIDE_HEAVY_CODE /** \return the Householder QR decomposition of \c *this. -- cgit v1.2.3