diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 6959a1c1b..b6bba04e6 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -134,7 +134,7 @@ template class SVD; template class JacobiSVD; template class LLT; template class LDLT; -template class HouseholderSequence; +template class HouseholderSequence; template class PlanarRotation; // Geometry module: diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 1e549633a..cfd3935fc 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Benoit Jacob +// Copyright (C) 2010 Benoit Jacob // Copyright (C) 2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index 7bd403373..61a8e8c1d 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2010 Benoit Jacob // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -48,31 +49,61 @@ * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ -template -struct ei_traits > +template +struct ei_traits > { typedef typename VectorsType::Scalar Scalar; enum { - RowsAtCompileTime = ei_traits::RowsAtCompileTime, - ColsAtCompileTime = ei_traits::RowsAtCompileTime, - MaxRowsAtCompileTime = ei_traits::MaxRowsAtCompileTime, - MaxColsAtCompileTime = ei_traits::MaxRowsAtCompileTime, + RowsAtCompileTime = Side==OnTheLeft ? ei_traits::RowsAtCompileTime + : ei_traits::ColsAtCompileTime, + ColsAtCompileTime = RowsAtCompileTime, + MaxRowsAtCompileTime = Side==OnTheLeft ? ei_traits::MaxRowsAtCompileTime + : ei_traits::MaxColsAtCompileTime, + MaxColsAtCompileTime = MaxRowsAtCompileTime, Flags = 0 }; }; -template class HouseholderSequence - : public AnyMatrixBase > +template +struct ei_hseq_side_dependent_impl +{ + typedef Block EssentialVectorType; + typedef HouseholderSequence HouseholderSequenceType; + static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, int k) + { + const int start = k+1+h.m_shift; + return Block(h.m_vectors, start, k, h.rows()-start, 1); + } +}; + +template +struct ei_hseq_side_dependent_impl +{ + typedef Transpose > EssentialVectorType; + typedef HouseholderSequence HouseholderSequenceType; + static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, int k) + { + const int start = k+1+h.m_shift; + return Block(h.m_vectors, k, start, 1, h.rows()-start).transpose(); + } +}; + +template class HouseholderSequence + : public AnyMatrixBase > { typedef typename VectorsType::Scalar Scalar; - typedef Block EssentialVectorType; + typedef typename ei_hseq_side_dependent_impl::EssentialVectorType + EssentialVectorType; public: - typedef HouseholderSequence::IsComplex, typename ei_cleantype::type, - CoeffsType>::ret> ConjugateReturnType; + CoeffsType>::ret, + Side + > ConjugateReturnType; HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()), @@ -85,14 +116,13 @@ template class HouseholderSequence { } - int rows() const { return m_vectors.rows(); } - int cols() const { return m_vectors.rows(); } + int rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); } + int cols() const { return rows(); } const EssentialVectorType essentialVector(int k) const { - ei_assert(k>= 0 && k < m_actualVectors); - const int start = k+1+m_shift; - return Block(m_vectors, start, k, rows()-start, 1); + ei_assert(k >= 0 && k < m_actualVectors); + return ei_hseq_side_dependent_impl::essentialVector(*this, k); } HouseholderSequence transpose() const @@ -164,6 +194,8 @@ template class HouseholderSequence return res; } + template friend struct ei_hseq_side_dependent_impl; + protected: typename VectorsType::Nested m_vectors; typename CoeffsType::Nested m_coeffs; @@ -175,13 +207,25 @@ template class HouseholderSequence template HouseholderSequence householderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) { - return HouseholderSequence(v, h, trans); + return HouseholderSequence(v, h, trans); } template HouseholderSequence householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift) { - return HouseholderSequence(v, h, trans, actualVectors, shift); + return HouseholderSequence(v, h, trans, actualVectors, shift); +} + +template +HouseholderSequence rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) +{ + return HouseholderSequence(v, h, trans); +} + +template +HouseholderSequence rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift) +{ + return HouseholderSequence(v, h, trans, actualVectors, shift); } #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/test/householder.cpp b/test/householder.cpp index c4c95e648..d98780872 100644 --- a/test/householder.cpp +++ b/test/householder.cpp @@ -106,27 +106,15 @@ template void householder(const MatrixType& m) m5.block(shift,0,brows,cols).template triangularView().setZero(); VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly m3 = hseq; - VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying + VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying -#if 0 // test householder sequence on the right with a shift - TMatrixType tm1 = m1.transpose(); TMatrixType tm2 = m2.transpose(); - - int bcols = cols - shift; - VBlockMatrixType vbm = - HouseholderQR qr(hbm); - m2 = m1; - m2.block(shift,0,brows,cols) = qr.matrixQR(); - HCoeffsVectorType hc = qr.hCoeffs().conjugate(); - HouseholderSequence hseq(m2, hc, false, hc.size(), shift); - MatrixType m5 = m2; - m5.block(shift,0,brows,cols).template triangularView().setZero(); - VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly - m3 = hseq; - VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying -#endif + HouseholderSequence rhseq(tm2, hc, false, hc.size(), shift); + VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly + m3 = rhseq; + VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying } void test_householder()