From 959a1b5d6335833e9ad49a088502705bb6967ff5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 25 Feb 2010 16:30:58 +0100 Subject: [PATCH] detect and implement inplace permutations --- Eigen/src/Core/PermutationMatrix.h | 49 ++++++++++++++++++++++-------- Eigen/src/Core/Transpose.h | 19 ------------ Eigen/src/Core/util/BlasUtil.h | 18 +++++++++++ Eigen/src/LU/FullPivLU.h | 8 ++--- Eigen/src/LU/PartialPivLU.h | 5 ++- test/permutationmatrices.cpp | 19 +++++++++++- 6 files changed, 78 insertions(+), 40 deletions(-) diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index c42812ec8..46884dc3f 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -326,21 +326,46 @@ struct ei_permut_matrix_product_retval template inline void evalTo(Dest& dst) const { const int n = Side==OnTheLeft ? rows() : cols(); - for(int i = 0; i < n; ++i) + + if(ei_is_same_type::ret && ei_extract_data(dst) == ei_extract_data(m_matrix)) { - Block< - Dest, - Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, - Side==OnTheRight ? 1 : Dest::ColsAtCompileTime - >(dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) + // apply the permutation inplace + Matrix mask(m_permutation.size()); + mask.fill(false); + int r = 0; + while(r < m_permutation.size()) + { + // search for the next seed + while(r=m_permutation.size()) + break; + // we got one, let's follow it until we are back to the seed + int k0 = r++; + int kPrev = k0; + mask.coeffRef(k0) = true; + for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k)) + { + Block(dst, k) + .swap(Block + (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev)); - = + mask.coeffRef(k) = true; + kPrev = k; + } + } + } + else + { + for(int i = 0; i < n; ++i) + { + Block + (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i) - Block< - MatrixTypeNestedCleaned, - Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime, - Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime - >(m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); + = + + Block + (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i); + } } } diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index bd06d8464..6c0e50de2 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -295,25 +295,6 @@ struct ei_blas_traits > static inline const XprType extract(const XprType& x) { return x; } }; - -template::ActualAccess> -struct ei_extract_data_selector { - static typename T::Scalar* run(const T& m) - { - return &ei_blas_traits::extract(m).const_cast_derived().coeffRef(0,0); - } -}; - -template -struct ei_extract_data_selector { - static typename T::Scalar* run(const T&) { return 0; } -}; - -template typename T::Scalar* ei_extract_data(const T& m) -{ - return ei_extract_data_selector::run(m); -} - template struct ei_check_transpose_aliasing_selector { diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index 2ca463d5d..4d216d77a 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -236,4 +236,22 @@ struct ei_blas_traits > static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); } }; +template::ActualAccess> +struct ei_extract_data_selector { + static const typename T::Scalar* run(const T& m) + { + return &ei_blas_traits::extract(m).const_cast_derived().coeffRef(0,0); // FIXME this should be .data() + } +}; + +template +struct ei_extract_data_selector { + static typename T::Scalar* run(const T&) { return 0; } +}; + +template const typename T::Scalar* ei_extract_data(const T& m) +{ + return ei_extract_data_selector::run(m); +} + #endif // EIGEN_BLASUTIL_H diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index cd63b9ec7..dea6ec41c 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -119,7 +119,7 @@ template class FullPivLU * diagonal coefficient of U. */ RealScalar maxPivot() const { return m_maxpivot; } - + /** \returns the permutation matrix P * * \sa permutationQ() @@ -506,12 +506,10 @@ MatrixType FullPivLU::reconstructedMatrix() const .template triangularView().toDenseMatrix(); // P^{-1}(LU) - // FIXME implement inplace permutation - res = (m_p.inverse() * res).eval(); + res = m_p.inverse() * res; // (P^{-1}LU)Q^{-1} - // FIXME implement inplace permutation - res = (res * m_q.inverse()).eval(); + res = res * m_q.inverse(); return res; } diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index fcffc2458..ad0d6b810 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -412,10 +412,9 @@ MatrixType PartialPivLU::reconstructedMatrix() const // LU MatrixType res = m_lu.template triangularView().toDenseMatrix() * m_lu.template triangularView(); - + // P^{-1}(LU) - // FIXME implement inplace permutation - res = (m_p.inverse() * res).eval(); + res = m_p.inverse() * res; return res; } diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index ae1bd8b85..89142d910 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp @@ -86,6 +86,23 @@ template void permutationmatrices(const MatrixType& m) identityp.setIdentity(rows); VERIFY_IS_APPROX(m_original, identityp*m_original); + // check inplace permutations + m_permuted = m_original; + m_permuted = lp.inverse() * m_permuted; + VERIFY_IS_APPROX(m_permuted, lp.inverse()*m_original); + + m_permuted = m_original; + m_permuted = m_permuted * rp.inverse(); + VERIFY_IS_APPROX(m_permuted, m_original*rp.inverse()); + + m_permuted = m_original; + m_permuted = lp * m_permuted; + VERIFY_IS_APPROX(m_permuted, lp*m_original); + + m_permuted = m_original; + m_permuted = m_permuted * rp; + VERIFY_IS_APPROX(m_permuted, m_original*rp); + if(rows>1 && cols>1) { lp2 = lp; @@ -114,7 +131,7 @@ void test_permutationmatrices() CALL_SUBTEST_2( permutationmatrices(Matrix3f()) ); CALL_SUBTEST_3( permutationmatrices(Matrix()) ); CALL_SUBTEST_4( permutationmatrices(Matrix4d()) ); - CALL_SUBTEST_5( permutationmatrices(Matrix()) ); + CALL_SUBTEST_5( permutationmatrices(Matrix()) ); CALL_SUBTEST_6( permutationmatrices(Matrix(20, 30)) ); CALL_SUBTEST_7( permutationmatrices(MatrixXcf(15, 10)) ); }