diff options
-rw-r--r-- | Eigen/src/Householder/HouseholderSequence.h | 69 | ||||
-rw-r--r-- | Eigen/src/QR/ColPivHouseholderQR.h | 5 | ||||
-rw-r--r-- | Eigen/src/QR/HouseholderQR.h | 6 | ||||
-rw-r--r-- | test/householder.cpp | 50 | ||||
-rw-r--r-- | test/qr.cpp | 7 |
5 files changed, 92 insertions, 45 deletions
diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index 26e0f6a88..7bd403373 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -65,32 +65,44 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> > { typedef typename VectorsType::Scalar Scalar; + typedef Block<VectorsType, Dynamic, 1> EssentialVectorType; + public: typedef HouseholderSequence<VectorsType, typename ei_meta_if<NumTraits<Scalar>::IsComplex, - NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >, + typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type, CoeffsType>::ret> ConjugateReturnType; HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) - : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()) - {} + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()), + m_shift(0) + { + } - HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) - : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors) - {} + HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift) + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors), m_shift(shift) + { + } int rows() const { return m_vectors.rows(); } int cols() const { return m_vectors.rows(); } + const EssentialVectorType essentialVector(int k) const + { + ei_assert(k>= 0 && k < m_actualVectors); + const int start = k+1+m_shift; + return Block<VectorsType,Dynamic,1>(m_vectors, start, k, rows()-start, 1); + } + HouseholderSequence transpose() const - { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); } + { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors, m_shift); } ConjugateReturnType conjugate() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors, m_shift); } ConjugateReturnType adjoint() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors, m_shift); } ConjugateReturnType inverse() const { return adjoint(); } @@ -98,45 +110,41 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence template<typename DestType> void evalTo(DestType& dst) const { int vecs = m_actualVectors; - int length = m_vectors.rows(); - dst.setIdentity(); - Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows()); + dst.setIdentity(rows(), rows()); + Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(rows()); for(int k = vecs-1; k >= 0; --k) { + int cornerSize = rows() - k - m_shift; if(m_trans) - dst.corner(BottomRight, length-k, length-k) - .applyHouseholderOnTheRight(m_vectors.col(k).tail(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + dst.corner(BottomRight, cornerSize, cornerSize) + .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); else - dst.corner(BottomRight, length-k, length-k) - .applyHouseholderOnTheLeft(m_vectors.col(k).tail(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k)); + dst.corner(BottomRight, cornerSize, cornerSize) + .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0)); } } /** \internal */ template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const { - int vecs = m_actualVectors; // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows()); - for(int k = 0; k < vecs; ++k) + for(int k = 0; k < m_actualVectors; ++k) { - int actual_k = m_trans ? vecs-k-1 : k; - dst.corner(BottomRight, dst.rows(), length-actual_k) - .applyHouseholderOnTheRight(m_vectors.col(actual_k).tail(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); + int actual_k = m_trans ? m_actualVectors-k-1 : k; + dst.corner(BottomRight, dst.rows(), rows()-m_shift-actual_k) + .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); } } /** \internal */ template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const { - int vecs = m_actualVectors; // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols()); - for(int k = 0; k < vecs; ++k) + for(int k = 0; k < m_actualVectors; ++k) { - int actual_k = m_trans ? k : vecs-k-1; - dst.corner(BottomRight, length-actual_k, dst.cols()) - .applyHouseholderOnTheLeft(m_vectors.col(actual_k).tail(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); + int actual_k = m_trans ? k : m_actualVectors-k-1; + dst.corner(BottomRight, rows()-m_shift-actual_k, dst.cols()) + .applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); } } @@ -161,6 +169,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence typename CoeffsType::Nested m_coeffs; bool m_trans; int m_actualVectors; + int m_shift; }; template<typename VectorsType, typename CoeffsType> @@ -170,9 +179,9 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp } template<typename VectorsType, typename CoeffsType> -HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) +HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift) { - return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors); + return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors, shift); } #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index f4edc9324..a119d0c2f 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -448,7 +448,8 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> dec().matrixQR(), dec().hCoeffs(), true, - dec().nonzeroPivots() + dec().nonzeroPivots(), + 0 )); dec().matrixQR() @@ -475,7 +476,7 @@ typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHousehol ::householderQ() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots); + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots, 0); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 3afaed9e0..b6b07ea63 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -55,11 +55,11 @@ template<typename _MatrixType> class HouseholderQR RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, - DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime) }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor> MatrixQType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; @@ -191,7 +191,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& RowVectorType temp(cols); - for (int k = 0; k < size; ++k) + for(int k = 0; k < size; ++k) { int remainingRows = rows - k; int remainingCols = cols - k - 1; diff --git a/test/householder.cpp b/test/householder.cpp index 4e4c78863..c4c95e648 100644 --- a/test/householder.cpp +++ b/test/householder.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -23,7 +23,7 @@ // Eigen. If not, see <http://www.gnu.org/licenses/>. #include "main.h" -#include <Eigen/Householder> +#include <Eigen/QR> template<typename MatrixType> void householder(const MatrixType& m) { @@ -40,7 +40,13 @@ template<typename MatrixType> void householder(const MatrixType& m) typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, ei_decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType; + typedef Matrix<Scalar, Dynamic, MatrixType::ColsAtCompileTime> HBlockMatrixType; + typedef Matrix<Scalar, Dynamic, 1> HCoeffsVectorType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RightSquareMatrixType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic> VBlockMatrixType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType; + Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols)); Scalar* tmp = &_tmp.coeffRef(0,0); @@ -85,8 +91,42 @@ template<typename MatrixType> void householder(const MatrixType& m) VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0))); VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha); - // test householder sequence - // TODO test HouseholderSequence + // test householder sequence on the left with a shift + + int shift = ei_random(0, std::max(rows-2,0)); + int brows = rows - shift; + m1.setRandom(rows, cols); + HBlockMatrixType hbm = m1.block(shift,0,brows,cols); + HouseholderQR<HBlockMatrixType> qr(hbm); + m2 = m1; + m2.block(shift,0,brows,cols) = qr.matrixQR(); + HCoeffsVectorType hc = qr.hCoeffs().conjugate(); + HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift); + MatrixType m5 = m2; + m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().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 + +#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<HBlockMatrixType> qr(hbm); + m2 = m1; + m2.block(shift,0,brows,cols) = qr.matrixQR(); + HCoeffsVectorType hc = qr.hCoeffs().conjugate(); + HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift); + MatrixType m5 = m2; + m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().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 } void test_householder() @@ -98,6 +138,6 @@ void test_householder() CALL_SUBTEST_4( householder(Matrix<float,4,4>()) ); CALL_SUBTEST_5( householder(MatrixXd(10,12)) ); CALL_SUBTEST_6( householder(MatrixXcf(16,17)) ); + CALL_SUBTEST_7( householder(MatrixXf(25,7)) ); } - } diff --git a/test/qr.cpp b/test/qr.cpp index 3848ce0a5..f2393c13b 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -36,14 +36,11 @@ template<typename MatrixType> void qr(const MatrixType& m) MatrixType a = MatrixType::Random(rows,cols); HouseholderQR<MatrixType> qrOfA(a); - MatrixType r = qrOfA.matrixQR(); MatrixQType q = qrOfA.householderQ(); VERIFY_IS_UNITARY(q); - // FIXME need better way to construct trapezoid - for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); - + MatrixType r = qrOfA.matrixQR().template triangularView<Upper>(); VERIFY_IS_APPROX(a, qrOfA.householderQ() * r); } @@ -114,7 +111,7 @@ template<typename MatrixType> void qr_verify_assert() void test_qr() { - for(int i = 0; i < 1; i++) { + for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( qr(MatrixXf(47,40)) ); CALL_SUBTEST_2( qr(MatrixXcd(17,7)) ); CALL_SUBTEST_3(( qr_fixedsize<Matrix<float,3,4>, 2 >() )); |