From 8b9dc9f0dfb44c2c4ad6d02fb88ecce0986cd154 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 9 Jan 2016 08:30:38 +0100 Subject: [PATCH] bug #1144: fix regression in x=y+A*x (aliasing), and move evaluator_traits::AssumeAliasing to evaluator_assume_aliasing. --- Eigen/src/Core/AssignEvaluator.h | 14 ++++---- Eigen/src/Core/CoreEvaluators.h | 8 ++--- Eigen/src/Core/ProductEvaluators.h | 25 ++++++-------- Eigen/src/Core/SelfAdjointView.h | 2 -- Eigen/src/Core/TriangularMatrix.h | 4 --- Eigen/src/Geometry/Homogeneous.h | 1 - Eigen/src/SparseCore/SparseSelfAdjointView.h | 2 -- Eigen/src/SparseQR/SparseQR.h | 1 - test/product.h | 35 +++++++++++++++----- test/product_large.cpp | 23 +++++++++++++ test/product_notemporary.cpp | 6 ++++ 11 files changed, 77 insertions(+), 44 deletions(-) diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index 9dfffbcc4..f6632de69 100755 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -682,9 +682,9 @@ template< typename DstXprType, typename SrcXprType, typename Functor, struct Assignment; -// The only purpose of this call_assignment() function is to deal with noalias() / AssumeAliasing and automatic transposition. -// Indeed, I (Gael) think that this concept of AssumeAliasing was a mistake, and it makes thing quite complicated. -// So this intermediate function removes everything related to AssumeAliasing such that Assignment +// The only purpose of this call_assignment() function is to deal with noalias() / "assume-aliasing" and automatic transposition. +// Indeed, I (Gael) think that this concept of "assume-aliasing" was a mistake, and it makes thing quite complicated. +// So this intermediate function removes everything related to "assume-aliasing" such that Assignment // does not has to bother about these annoying details. template @@ -698,21 +698,21 @@ EIGEN_DEVICE_FUNC void call_assignment(const Dst& dst, const Src& src) call_assignment(dst, src, internal::assign_op()); } -// Deal with AssumeAliasing +// Deal with "assume-aliasing" template -EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::AssumeAliasing==1, void*>::type = 0) +EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if< evaluator_assume_aliasing::value, void*>::type = 0) { typename plain_matrix_type::type tmp(src); call_assignment_no_alias(dst, tmp, func); } template -EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::AssumeAliasing==0, void*>::type = 0) +EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if::value, void*>::type = 0) { call_assignment_no_alias(dst, src, func); } -// by-pass AssumeAliasing +// by-pass "assume-aliasing" // When there is no aliasing, we require that 'dst' has been properly resized template class StorageBase, typename Src, typename Func> EIGEN_DEVICE_FUNC void call_assignment(NoAlias& dst, const Src& src, const Func& func) diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index f97dc33de..8bd73b814 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -63,10 +63,6 @@ struct evaluator_traits_base // by default, get evaluator kind and shape from storage typedef typename storage_kind_to_evaluator_kind::StorageKind>::Kind Kind; typedef typename storage_kind_to_shape::StorageKind>::Shape Shape; - - // 1 if assignment A = B assumes aliasing when B is of type T and thus B needs to be evaluated into a - // temporary; 0 if not. - static const int AssumeAliasing = 0; }; // Default evaluator traits @@ -75,6 +71,10 @@ struct evaluator_traits : public evaluator_traits_base { }; +template::Shape > +struct evaluator_assume_aliasing { + static const bool value = false; +}; // By default, we assume a unary expression: template diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index 794038a2a..b2a0a4b4f 100755 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h @@ -38,10 +38,9 @@ struct evaluator > // Catch scalar * ( A * B ) and transform it to (A*scalar) * B // TODO we should apply that rule only if that's really helpful template -struct evaluator_traits, const Product > > - : evaluator_traits_base, const Product > > +struct evaluator_assume_aliasing, const Product > > { - enum { AssumeAliasing = 1 }; + static const bool value = true; }; template struct evaluator, const Product > > @@ -81,17 +80,8 @@ template< typename Lhs, typename Rhs, struct generic_product_impl; template -struct evaluator_traits > - : evaluator_traits_base > -{ - enum { AssumeAliasing = 1 }; -}; - -template -struct evaluator_traits > - : evaluator_traits_base > -{ - enum { AssumeAliasing = 0 }; +struct evaluator_assume_aliasing > { + static const bool value = true; }; // This is the default evaluator implementation for products: @@ -189,6 +179,13 @@ struct Assignment" expression to save one temporary // FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct +// TODO enable it for "Dense ?= xpr - Product<>" as well. + +template +struct evaluator_assume_aliasing, const OtherXpr, + const Product >, DenseShape > { + static const bool value = true; +}; template struct assignment_from_xpr_plus_product diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 87e87ab3a..e709eb213 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -203,8 +203,6 @@ struct evaluator_traits > { typedef typename storage_kind_to_evaluator_kind::Kind Kind; typedef SelfAdjointShape Shape; - - static const int AssumeAliasing = 0; }; template diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 7d6a97848..27845e89c 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -704,10 +704,6 @@ struct evaluator_traits > { typedef typename storage_kind_to_evaluator_kind::Kind Kind; typedef typename glue_shapes::Shape, TriangularShape>::type Shape; - - // 1 if assignment A = B assumes aliasing when B is of type T and thus B needs to be evaluated into a - // temporary; 0 if not. - static const int AssumeAliasing = 0; }; template diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 367fd3930..cd52b5470 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -304,7 +304,6 @@ struct evaluator_traits > { typedef typename storage_kind_to_evaluator_kind::Kind Kind; typedef HomogeneousShape Shape; - static const int AssumeAliasing = 0; }; template<> struct AssignmentKind { typedef Dense2Dense Kind; }; diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 46c6ce1d3..975cefd28 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -211,8 +211,6 @@ struct evaluator_traits > { typedef typename storage_kind_to_evaluator_kind::Kind Kind; typedef SparseSelfAdjointShape Shape; - - static const int AssumeAliasing = 0; }; struct SparseSelfAdjoint2Sparse {}; diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 4f26c19ca..0d448d02e 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -691,7 +691,6 @@ struct evaluator_traits > typedef typename SparseQRType::MatrixType MatrixType; typedef typename storage_kind_to_evaluator_kind::Kind Kind; typedef SparseShape Shape; - static const int AssumeAliasing = 0; }; template< typename DstXprType, typename SparseQRType> diff --git a/test/product.h b/test/product.h index 9dfff9303..bd92309d2 100644 --- a/test/product.h +++ b/test/product.h @@ -145,14 +145,31 @@ template void product(const MatrixType& m) VERIFY_IS_APPROX(res.col(r).noalias() = square * square.col(r), (square * square.col(r)).eval()); // inner product - Scalar x = square2.row(c) * square2.col(c2); - VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum()); - + { + Scalar x = square2.row(c) * square2.col(c2); + VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum()); + } + // outer product - VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); - VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose()); - VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); - VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); - VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols)); - VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols)); + { + VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose()); + VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols)); + VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols)); + } + + // Aliasing + { + ColVectorType x(cols); x.setRandom(); + ColVectorType z(x); + ColVectorType y(cols); y.setZero(); + ColSquareMatrixType A(cols,cols); A.setRandom(); + // CwiseBinaryOp + VERIFY_IS_APPROX(x = y + A*x, A*z); + x = z; + // CwiseUnaryOp + VERIFY_IS_APPROX(x = Scalar(1.)*(A*x), A*z); + } } diff --git a/test/product_large.cpp b/test/product_large.cpp index 7207973c2..98f84c53b 100644 --- a/test/product_large.cpp +++ b/test/product_large.cpp @@ -9,6 +9,27 @@ #include "product.h" +template +void test_aliasing() +{ + int rows = internal::random(1,12); + int cols = internal::random(1,12); + typedef Matrix MatrixType; + typedef Matrix VectorType; + VectorType x(cols); x.setRandom(); + VectorType z(x); + VectorType y(rows); y.setZero(); + MatrixType A(rows,cols); A.setRandom(); + // CwiseBinaryOp + VERIFY_IS_APPROX(x = y + A*x, A*z); // OK because "y + A*x" is marked as "assume-aliasing" + x = z; + // CwiseUnaryOp + VERIFY_IS_APPROX(x = T(1.)*(A*x), A*z); // OK because 1*(A*x) is replaced by (1*A*x) which is a Product<> expression + x = z; + // VERIFY_IS_APPROX(x = y-A*x, -A*z); // Not OK in 3.3 because x is resized before A*x gets evaluated + x = z; +} + void test_product_large() { for(int i = 0; i < g_repeat; i++) { @@ -17,6 +38,8 @@ void test_product_large() CALL_SUBTEST_3( product(MatrixXi(internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_4( product(MatrixXcf(internal::random(1,EIGEN_TEST_MAX_SIZE/2), internal::random(1,EIGEN_TEST_MAX_SIZE/2))) ); CALL_SUBTEST_5( product(Matrix(internal::random(1,EIGEN_TEST_MAX_SIZE), internal::random(1,EIGEN_TEST_MAX_SIZE))) ); + + CALL_SUBTEST_1( test_aliasing() ); } #if defined EIGEN_TEST_PART_6 diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index ff93cb881..5a3f3a01a 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -43,10 +43,16 @@ template void product_notemporary(const MatrixType& m) r1 = internal::random(8,rows-r0); VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1); + VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).transpose(), 1); VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3 = s1 * (m1 * m2.transpose()), 1); +// VERIFY_EVALUATION_COUNT( m3 = m3 + s1 * (m1 * m2.transpose()), 1); VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0); + VERIFY_EVALUATION_COUNT( m3 = m3 + (m1 * m2.adjoint()), 1); + + VERIFY_EVALUATION_COUNT( m3 = m3 + (m1 * m2.adjoint()).transpose(), 1); VERIFY_EVALUATION_COUNT( m3.noalias() = m3 + m1 * m2.transpose(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() += m3 + m1 * m2.transpose(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() -= m3 + m1 * m2.transpose(), 0);