diff options
| -rw-r--r-- | Eigen/src/Core/VectorBlock.h | 3 | ||||
| -rw-r--r-- | Eigen/src/Householder/HouseholderSequence.h | 26 | ||||
| -rw-r--r-- | Eigen/src/QR/ColPivotingHouseholderQR.h | 59 | ||||
| -rw-r--r-- | Eigen/src/QR/FullPivotingHouseholderQR.h | 26 | ||||
| -rw-r--r-- | test/qr_colpivoting.cpp | 8 | ||||
| -rw-r--r-- | test/qr_fullpivoting.cpp | 8 |
6 files changed, 68 insertions, 62 deletions
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<typename VectorType, int Size, int PacketAccess> class VectorBlock typedef Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, - PacketAccess> Base; + PacketAccess> _Base; enum { IsColVector = ei_traits<VectorType>::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<typename VectorsType, typename CoeffsType> 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<typename VectorsType, typename CoeffsType> class HouseholderSequence } } + template<typename OtherDerived> + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other) const + { + typename OtherDerived::PlainMatrixType res(other); + applyThisOnTheLeft(res); + return res; + } + + template<typename OtherDerived> friend + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& 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<typename VectorsType, typename CoeffsType> class HouseholderSequence bool m_trans; }; +template<typename VectorsType, typename CoeffsType> +HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) +{ + return HouseholderSequence<VectorsType,CoeffsType>(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<typename MatrixType> 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<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -62,6 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType; + typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -99,7 +100,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR template<typename OtherDerived, typename ResultType> bool solve(const MatrixBase<OtherDerived>& 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<typename MatrixType> 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<typename MatrixType> 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<MatrixType>& ColPivotingHouseholderQR<MatrixType>::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<MatrixType>& ColPivotingHouseholderQR<MatrixType>::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<MatrixType>& ColPivotingHouseholderQR<MatrixType>::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<MatrixType>& ColPivotingHouseholderQR<MatrixType>::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<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -352,16 +353,11 @@ bool ColPivotingHouseholderQR<MatrixType>::solve( const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - - Matrix<Scalar,1,MatrixType::ColsAtCompileTime> 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<MatrixType>::solve( return true; } -/** \returns the matrix Q */ +/** \returns the matrix Q as a sequence of householder transformations */ template<typename MatrixType> -typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const +typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::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<Scalar,1,MatrixType::RowsAtCompileTime> 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<typename MatrixType> 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<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -106,13 +106,13 @@ template<typename MatrixType> 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<typename MatrixType> 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<MatrixType>& FullPivotingHouseholderQR<MatrixType>::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<MatrixType>& FullPivotingHouseholderQR<MatrixType>::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<MatrixType>& FullPivotingHouseholderQR<MatrixType>::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<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -358,13 +358,13 @@ bool FullPivotingHouseholderQR<MatrixType>::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<Scalar,1,MatrixType::ColsAtCompileTime> 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<typename MatrixType> 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<typename MatrixType> 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<typename MatrixType> 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<Scalar>(); |