From d9ca0c0d3643f4b777de686a2c0cddde075aa063 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 25 Feb 2010 15:31:15 +0100 Subject: [PATCH] optimize inverse permutations --- Eigen/src/Core/PermutationMatrix.h | 139 +++++++++++++++++++++++++---- test/permutationmatrices.cpp | 4 +- 2 files changed, 122 insertions(+), 21 deletions(-) diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index fcd2e46cc..c42812ec8 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -47,7 +47,7 @@ * \sa class DiagonalMatrix */ template class PermutationMatrix; -template struct ei_permut_matrix_product_retval; +template struct ei_permut_matrix_product_retval; template struct ei_traits > @@ -132,6 +132,9 @@ class PermutationMatrix : public EigenBase void evalTo(MatrixBase& other) const @@ -213,16 +216,29 @@ class PermutationMatrix : public EigenBase inverse() const + { return *this; } + /** \returns the tranpose permutation matrix. + * + * \note \note_try_to_help_rvo + */ + inline Transpose transpose() const + { return *this; } + + /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN - protected: - enum Inverse_t {Inverse}; - PermutationMatrix(Inverse_t, const PermutationMatrix& other) - : m_indices(other.m_indices.size()) + template + PermutationMatrix(const Transpose >& other) + : m_indices(other.nestedPermutation().size()) { - for (int i=0; i& other) const { return PermutationMatrix(Product, *this, other); } + /** \returns the product of a permutation with another inverse permutation. + * + * \note \note_try_to_help_rvo + */ + template + inline PermutationMatrix operator*(const Transpose >& other) const + { return PermutationMatrix(Product, *this, other.eval()); } + + /** \returns the product of an inverse permutation with another permutation. + * + * \note \note_try_to_help_rvo + */ + template friend + inline PermutationMatrix operator*(const Transpose >& other, const PermutationMatrix& perm) + { return PermutationMatrix(Product, other.eval(), perm); } + protected: IndicesType m_indices; @@ -277,15 +304,15 @@ operator*(const PermutationMatrix &perm (permutation, matrix.derived()); } -template -struct ei_traits > +template +struct ei_traits > { typedef typename MatrixType::PlainObject ReturnType; }; -template +template struct ei_permut_matrix_product_retval - : public ReturnByValue > + : public ReturnByValue > { typedef typename ei_cleantype::type MatrixTypeNestedCleaned; @@ -305,7 +332,7 @@ struct ei_permut_matrix_product_retval Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime - >(dst, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i) + >(dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) = @@ -313,7 +340,7 @@ struct ei_permut_matrix_product_retval MatrixTypeNestedCleaned, Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime, Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime - >(m_matrix, Side==OnTheRight ? m_permutation.indices().coeff(i) : i); + >(m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); } } @@ -322,4 +349,78 @@ struct ei_permut_matrix_product_retval const typename MatrixType::Nested m_matrix; }; +/* Template partial specialization for transposed/inverse permutations */ + +template +struct ei_traits > > + : ei_traits > +{}; + +template +class Transpose > + : public EigenBase > > +{ + typedef PermutationMatrix PermutationType; + typedef typename PermutationType::IndicesType IndicesType; + public: + + #ifndef EIGEN_PARSED_BY_DOXYGEN + typedef ei_traits Traits; + typedef Matrix + DenseMatrixType; + enum { + Flags = Traits::Flags, + CoeffReadCost = Traits::CoeffReadCost, + RowsAtCompileTime = Traits::RowsAtCompileTime, + ColsAtCompileTime = Traits::ColsAtCompileTime, + MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = Traits::MaxColsAtCompileTime + }; + typedef typename Traits::Scalar Scalar; + #endif + + Transpose(const PermutationType& p) : m_permutation(p) {} + + inline int rows() const { return m_permutation.rows(); } + inline int cols() const { return m_permutation.cols(); } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template + void evalTo(MatrixBase& other) const + { + other.setZero(); + for (int i=0; i friend + inline const ei_permut_matrix_product_retval + operator*(const MatrixBase& matrix, const Transpose& trPerm) + { + return ei_permut_matrix_product_retval(trPerm.m_permutation, matrix.derived()); + } + + /** \returns the matrix with the inverse permutation applied to the rows. + */ + template + inline const ei_permut_matrix_product_retval + operator*(const MatrixBase& matrix) const + { + return ei_permut_matrix_product_retval(m_permutation, matrix.derived()); + } + + const PermutationType& nestedPermutation() const { return m_permutation; } + + protected: + const PermutationType& m_permutation; +}; + #endif // EIGEN_PERMUTATIONMATRIX_H diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index 0ef0a371a..ae1bd8b85 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp @@ -51,7 +51,7 @@ template void permutationmatrices(const MatrixType& m) typedef Matrix LeftPermutationVectorType; typedef PermutationMatrix RightPermutationType; typedef Matrix RightPermutationVectorType; - + int rows = m.rows(); int cols = m.cols(); @@ -72,7 +72,7 @@ template void permutationmatrices(const MatrixType& m) Matrix rm(rp); VERIFY_IS_APPROX(m_permuted, lm*m_original*rm); - + VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original); VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity());