GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-10 00:03:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make permutation products aliasing by default.

Browse Source
This commit is contained in:
Antonio Sanchez 2025-08-25 05:28:25 -07:00 committed by Rasmus Munk Larsen
parent 4ae5647355
commit d2a70fe4e2
2 changed files with 33 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
Eigen/src/Core/PermutationMatrix.h
View File

@ -468,17 +468,17 @@ class PermutationWrapper : public PermutationBase<PermutationWrapper<IndicesType
/** \returns the matrix with the permutation applied to the columns.
*/
template <typename MatrixDerived, typename PermutationDerived>
EIGEN_DEVICE_FUNC const Product<MatrixDerived, PermutationDerived, AliasFreeProduct> operator*(
EIGEN_DEVICE_FUNC const Product<MatrixDerived, PermutationDerived, DefaultProduct> operator*(
const MatrixBase<MatrixDerived>& matrix, const PermutationBase<PermutationDerived>& permutation) {
return Product<MatrixDerived, PermutationDerived, AliasFreeProduct>(matrix.derived(), permutation.derived());
return Product<MatrixDerived, PermutationDerived, DefaultProduct>(matrix.derived(), permutation.derived());
}
/** \returns the matrix with the permutation applied to the rows.
*/
template <typename PermutationDerived, typename MatrixDerived>
EIGEN_DEVICE_FUNC const Product<PermutationDerived, MatrixDerived, AliasFreeProduct> operator*(
EIGEN_DEVICE_FUNC const Product<PermutationDerived, MatrixDerived, DefaultProduct> operator*(
const PermutationBase<PermutationDerived>& permutation, const MatrixBase<MatrixDerived>& matrix) {
return Product<PermutationDerived, MatrixDerived, AliasFreeProduct>(permutation.derived(), matrix.derived());
return Product<PermutationDerived, MatrixDerived, DefaultProduct>(permutation.derived(), matrix.derived());
}
template <typename PermutationType>
@ -520,16 +520,16 @@ class InverseImpl<PermutationType, PermutationStorage> : public EigenBase<Invers
/** \returns the matrix with the inverse permutation applied to the columns.
*/
template <typename OtherDerived>
friend const Product<OtherDerived, InverseType, AliasFreeProduct> operator*(const MatrixBase<OtherDerived>& matrix,
const InverseType& trPerm) {
return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
friend const Product<OtherDerived, InverseType, DefaultProduct> operator*(const MatrixBase<OtherDerived>& matrix,
const InverseType& trPerm) {
return Product<OtherDerived, InverseType, DefaultProduct>(matrix.derived(), trPerm.derived());
}
/** \returns the matrix with the inverse permutation applied to the rows.
*/
template <typename OtherDerived>
const Product<InverseType, OtherDerived, AliasFreeProduct> operator*(const MatrixBase<OtherDerived>& matrix) const {
return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
const Product<InverseType, OtherDerived, DefaultProduct> operator*(const MatrixBase<OtherDerived>& matrix) const {
return Product<InverseType, OtherDerived, DefaultProduct>(derived(), matrix.derived());
}
};

26
test/permutationmatrices.cpp
View File

@ -39,7 +39,8 @@ void permutationmatrices(const MatrixType& m) {
RightTranspositionsType rt(rv);
MatrixType m_permuted = MatrixType::Random(rows, cols);
VERIFY_EVALUATION_COUNT(m_permuted = lp * m_original * rp, 1); // 1 temp for sub expression "lp * m_original"
VERIFY_EVALUATION_COUNT(m_permuted.noalias() = lp * m_original * rp,
1); // 1 temp for sub expression "lp * m_original"
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++) VERIFY_IS_APPROX(m_permuted(lv(i), j), m_original(i, rv(j)));
@ -50,7 +51,7 @@ void permutationmatrices(const MatrixType& m) {
VERIFY_IS_APPROX(m_permuted, lm * m_original * rm);
m_permuted = m_original;
VERIFY_EVALUATION_COUNT(m_permuted = lp * m_permuted * rp, 1);
VERIFY_EVALUATION_COUNT(m_permuted.noalias() = lp * m_permuted * rp, 1);
VERIFY_IS_APPROX(m_permuted, lm * m_original * rm);
LeftPermutationType lpi;
@ -169,6 +170,26 @@ void bug890() {
VERIFY_IS_APPROX(v1, (P.inverse() * rhs).eval());
}
void test_aliasing() {
// Bug #2869.
auto P = Eigen::PermutationMatrix<4>{Eigen::Vector4i{0, 2, 3, 1}};
{
Eigen::Vector<float, 5> z = {0.0f, 1.1f, 2.2f, 3.3f, 4.4f};
Eigen::Vector<float, 5> expected = z;
z.tail<4>() = P * z.head<4>();
expected.tail<4>() = (P * expected.head<4>()).eval();
VERIFY_IS_APPROX(z, expected);
}
{
// Obfuscate the aliasing in the RHS expression.
Eigen::Vector4f a = {1.1f, 2.2f, 3.3f, 4.4f};
Eigen::Vector4f expected = P * (a + Eigen::Vector4f::Zero()).cwiseSqrt();
a = P * (a + Eigen::Vector4f::Zero()).cwiseSqrt();
VERIFY_IS_APPROX(a, expected);
}
}
EIGEN_DECLARE_TEST(permutationmatrices) {
for (int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1(permutationmatrices(Matrix<float, 1, 1>()));
@ -182,4 +203,5 @@ EIGEN_DECLARE_TEST(permutationmatrices) {
MatrixXcf(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
}
CALL_SUBTEST_5(bug890<double>());
CALL_SUBTEST_4(test_aliasing());
}