diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h index b291f7b1a..65268b626 100644 --- a/Eigen/src/Core/VectorBlock.h +++ b/Eigen/src/Core/VectorBlock.h @@ -77,11 +77,12 @@ template class VectorBlock typedef Block::RowsAtCompileTime==1 ? 1 : Size, ei_traits::ColsAtCompileTime==1 ? 1 : Size, - PacketAccess> Base; + PacketAccess> _Base; enum { IsColVector = ei_traits::ColsAtCompileTime==1 }; public: + _EIGEN_GENERIC_PUBLIC_INTERFACE(VectorBlock, _Base) using Base::operator=; using Base::operator+=; diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index 86395213b..16e362814 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -80,10 +80,10 @@ template class HouseholderSequence 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); } ConjugateReturnType conjugate() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } ConjugateReturnType adjoint() const { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); } @@ -136,6 +136,22 @@ template class HouseholderSequence } } + template + typename OtherDerived::PlainMatrixType operator*(const MatrixBase& other) const + { + typename OtherDerived::PlainMatrixType res(other); + applyThisOnTheLeft(res); + return res; + } + + template friend + typename OtherDerived::PlainMatrixType operator*(const MatrixBase& other, const HouseholderSequence& h) + { + typename OtherDerived::PlainMatrixType res(other); + h.applyThisOnTheRight(res); + return res; + } + protected: typename VectorsType::Nested m_vectors; @@ -143,4 +159,10 @@ template class HouseholderSequence bool m_trans; }; +template +HouseholderSequence makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) +{ + return HouseholderSequence(v, h, trans); +} + #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h index 883e3449f..c4c7d2d55 100644 --- a/Eigen/src/QR/ColPivotingHouseholderQR.h +++ b/Eigen/src/QR/ColPivotingHouseholderQR.h @@ -45,14 +45,14 @@ template class ColPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix MatrixQType; @@ -62,6 +62,7 @@ template class ColPivotingHouseholderQR typedef Matrix RowVectorType; typedef Matrix ColVectorType; typedef Matrix RealRowVectorType; + typedef typename HouseholderSequence::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -99,7 +100,7 @@ template class ColPivotingHouseholderQR template bool solve(const MatrixBase& b, ResultType *result) const; - MatrixQType matrixQ(void) const; + HouseholderSequenceType matrixQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -110,13 +111,13 @@ template class ColPivotingHouseholderQR } ColPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + /** \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) @@ -145,7 +146,7 @@ template class ColPivotingHouseholderQR * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -268,7 +269,7 @@ ColPivotingHouseholderQR& ColPivotingHouseholderQR::comp int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -279,18 +280,18 @@ ColPivotingHouseholderQR& ColPivotingHouseholderQR::comp IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealRowVectorType colSqNorms(cols); for(int k = 0; k < cols; ++k) colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); - + for (int k = 0; k < size; ++k) { int biggest_col_in_corner; RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); biggest_col_in_corner += k; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) { @@ -302,7 +303,7 @@ ColPivotingHouseholderQR& ColPivotingHouseholderQR::comp } break; } - + cols_transpositions.coeffRef(k) = biggest_col_in_corner; if(k != biggest_col_in_corner) { m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); @@ -316,7 +317,7 @@ ColPivotingHouseholderQR& ColPivotingHouseholderQR::comp m_qr.corner(BottomRight, rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); - + colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); } @@ -326,7 +327,7 @@ ColPivotingHouseholderQR& ColPivotingHouseholderQR::comp m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -352,16 +353,11 @@ bool ColPivotingHouseholderQR::solve( const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - - Matrix temp(cols); - for (int k = 0; k < m_rank; ++k) - { - int remainingSize = rows-k; - c.corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); - } + + // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T + c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose()); if(!isSurjective()) { @@ -381,25 +377,12 @@ bool ColPivotingHouseholderQR::solve( return true; } -/** \returns the matrix Q */ +/** \returns the matrix Q as a sequence of householder transformations */ template -typename ColPivotingHouseholderQR::MatrixQType ColPivotingHouseholderQR::matrixQ() const +typename ColPivotingHouseholderQR::HouseholderSequenceType ColPivotingHouseholderQR::matrixQ() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR 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), ...] - int rows = m_qr.rows(); - int cols = m_qr.cols(); - int size = std::min(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix temp(rows); - for (int k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - } - return res; + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index 0d542cf7a..9fee77803 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -45,14 +45,14 @@ template class FullPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix MatrixQType; @@ -106,13 +106,13 @@ template class FullPivotingHouseholderQR } FullPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + const IntColVectorType& rowsTranspositions() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); @@ -147,7 +147,7 @@ template class FullPivotingHouseholderQR * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -271,7 +271,7 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -283,9 +283,9 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealScalar biggest(0); - + for (int k = 0; k < size; ++k) { int row_of_biggest_in_corner, col_of_biggest_in_corner; @@ -297,7 +297,7 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co row_of_biggest_in_corner += k; col_of_biggest_in_corner += k; if(k==0) biggest = biggest_in_corner; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) { @@ -336,7 +336,7 @@ FullPivotingHouseholderQR& FullPivotingHouseholderQR::co m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -358,13 +358,13 @@ bool FullPivotingHouseholderQR::solve( } else return false; } - + const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - + Matrix temp(cols); for (int k = 0; k < m_rank; ++k) { diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 9c387005a..588a41e56 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -42,18 +42,18 @@ template void qr() VERIFY(!qr.isInjective()); VERIFY(!qr.isInvertible()); VERIFY(!qr.isSurjective()); - + MatrixType r = qr.matrixQR(); // 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 b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 525c669a5..3a37bcb46 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -46,14 +46,14 @@ template void qr() MatrixType r = qr.matrixQR(); // 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 b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); @@ -88,7 +88,7 @@ template void qr_invertible() m3 = MatrixType::Random(size,size); VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); - + // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = ei_random();