[eigen3] 01/02: Fix regression in x=y+A*x (aliasing).
Anton Gladky
gladk at moszumanska.debian.org
Sun Jan 10 16:01:11 UTC 2016
This is an automated email from the git hooks/post-receive script.
gladk pushed a commit to branch master
in repository eigen3.
commit 5b7de84bafe0b4ae2fe071a25febedf1ca0c7ae2
Author: Gael Guennebaud <g.gael at free.fr>
Date: Sun Jan 10 16:35:17 2016 +0100
Fix regression in x=y+A*x (aliasing).
---
debian/patches/09_fix_1144.patch | 299 +++++++++++++++++++++++++++++++++++++++
debian/patches/series | 1 +
2 files changed, 300 insertions(+)
diff --git a/debian/patches/09_fix_1144.patch b/debian/patches/09_fix_1144.patch
new file mode 100644
index 0000000..ac3bb67
--- /dev/null
+++ b/debian/patches/09_fix_1144.patch
@@ -0,0 +1,299 @@
+# HG changeset patch
+# User Gael Guennebaud <g.gael at free.fr>
+# Date 1452324638 -3600
+# Node ID 0db197c797c3891768bd8d74fcf85279680942a5
+# Parent dd88c44d411c37a833ecbc0f3e82f58eb801107a
+Bug 1144: fix regression in x=y+A*x (aliasing), and move evaluator_traits::AssumeAliasing to evaluator_assume_aliasing.
+
+Index: eigen3/Eigen/src/Core/AssignEvaluator.h
+===================================================================
+--- eigen3.orig/Eigen/src/Core/AssignEvaluator.h
++++ eigen3/Eigen/src/Core/AssignEvaluator.h
+@@ -682,9 +682,9 @@ template< typename DstXprType, typename
+ 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<typename Dst, typename Src>
+@@ -698,21 +698,21 @@ EIGEN_DEVICE_FUNC void call_assignment(c
+ call_assignment(dst, src, internal::assign_op<typename Dst::Scalar>());
+ }
+
+-// Deal with AssumeAliasing
++// Deal with "assume-aliasing"
+ template<typename Dst, typename Src, typename Func>
+-EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if<evaluator_traits<Src>::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<Src>::value, void*>::type = 0)
+ {
+ typename plain_matrix_type<Src>::type tmp(src);
+ call_assignment_no_alias(dst, tmp, func);
+ }
+
+ template<typename Dst, typename Src, typename Func>
+-EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if<evaluator_traits<Src>::AssumeAliasing==0, void*>::type = 0)
++EIGEN_DEVICE_FUNC void call_assignment(Dst& dst, const Src& src, const Func& func, typename enable_if<!evaluator_assume_aliasing<Src>::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<typename Dst, template <typename> class StorageBase, typename Src, typename Func>
+ EIGEN_DEVICE_FUNC void call_assignment(NoAlias<Dst,StorageBase>& dst, const Src& src, const Func& func)
+Index: eigen3/Eigen/src/Core/CoreEvaluators.h
+===================================================================
+--- eigen3.orig/Eigen/src/Core/CoreEvaluators.h
++++ eigen3/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<typename traits<T>::StorageKind>::Kind Kind;
+ typedef typename storage_kind_to_shape<typename traits<T>::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 evaluat
+ {
+ };
+
++template<typename T, typename Shape = typename evaluator_traits<T>::Shape >
++struct evaluator_assume_aliasing {
++ static const bool value = false;
++};
+
+ // By default, we assume a unary expression:
+ template<typename T>
+Index: eigen3/Eigen/src/Core/ProductEvaluators.h
+===================================================================
+--- eigen3.orig/Eigen/src/Core/ProductEvaluators.h
++++ eigen3/Eigen/src/Core/ProductEvaluators.h
+@@ -38,10 +38,9 @@ struct evaluator<Product<Lhs, Rhs, Optio
+ // Catch scalar * ( A * B ) and transform it to (A*scalar) * B
+ // TODO we should apply that rule only if that's really helpful
+ template<typename Lhs, typename Rhs, typename Scalar>
+-struct evaluator_traits<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const Product<Lhs, Rhs, DefaultProduct> > >
+- : evaluator_traits_base<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const Product<Lhs, Rhs, DefaultProduct> > >
++struct evaluator_assume_aliasing<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const Product<Lhs, Rhs, DefaultProduct> > >
+ {
+- enum { AssumeAliasing = 1 };
++ static const bool value = true;
+ };
+ template<typename Lhs, typename Rhs, typename Scalar>
+ struct evaluator<CwiseUnaryOp<internal::scalar_multiple_op<Scalar>, const Product<Lhs, Rhs, DefaultProduct> > >
+@@ -81,17 +80,8 @@ template< typename Lhs, typename Rhs,
+ struct generic_product_impl;
+
+ template<typename Lhs, typename Rhs>
+-struct evaluator_traits<Product<Lhs, Rhs, DefaultProduct> >
+- : evaluator_traits_base<Product<Lhs, Rhs, DefaultProduct> >
+-{
+- enum { AssumeAliasing = 1 };
+-};
+-
+-template<typename Lhs, typename Rhs>
+-struct evaluator_traits<Product<Lhs, Rhs, AliasFreeProduct> >
+- : evaluator_traits_base<Product<Lhs, Rhs, AliasFreeProduct> >
+-{
+- enum { AssumeAliasing = 0 };
++struct evaluator_assume_aliasing<Product<Lhs, Rhs, DefaultProduct> > {
++ static const bool value = true;
+ };
+
+ // This is the default evaluator implementation for products:
+@@ -189,6 +179,13 @@ struct Assignment<DstXprType, CwiseUnary
+ //----------------------------------------
+ // Catch "Dense ?= xpr + Product<>" 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<typename OtherXpr, typename Lhs, typename Rhs>
++struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar>, const OtherXpr,
++ const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
++ static const bool value = true;
++};
+
+ template<typename DstXprType, typename OtherXpr, typename ProductType, typename Scalar, typename Func1, typename Func2>
+ struct assignment_from_xpr_plus_product
+Index: eigen3/Eigen/src/Core/SelfAdjointView.h
+===================================================================
+--- eigen3.orig/Eigen/src/Core/SelfAdjointView.h
++++ eigen3/Eigen/src/Core/SelfAdjointView.h
+@@ -203,8 +203,6 @@ struct evaluator_traits<SelfAdjointView<
+ {
+ typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
+ typedef SelfAdjointShape Shape;
+-
+- static const int AssumeAliasing = 0;
+ };
+
+ template<int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version>
+Index: eigen3/Eigen/src/Core/TriangularMatrix.h
+===================================================================
+--- eigen3.orig/Eigen/src/Core/TriangularMatrix.h
++++ eigen3/Eigen/src/Core/TriangularMatrix.h
+@@ -711,10 +711,6 @@ struct evaluator_traits<TriangularView<M
+ {
+ typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
+ typedef typename glue_shapes<typename evaluator_traits<MatrixType>::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<typename MatrixType, unsigned int Mode>
+Index: eigen3/Eigen/src/Geometry/Homogeneous.h
+===================================================================
+--- eigen3.orig/Eigen/src/Geometry/Homogeneous.h
++++ eigen3/Eigen/src/Geometry/Homogeneous.h
+@@ -304,7 +304,6 @@ struct evaluator_traits<Homogeneous<ArgT
+ {
+ typedef typename storage_kind_to_evaluator_kind<typename ArgType::StorageKind>::Kind Kind;
+ typedef HomogeneousShape Shape;
+- static const int AssumeAliasing = 0;
+ };
+
+ template<> struct AssignmentKind<DenseShape,HomogeneousShape> { typedef Dense2Dense Kind; };
+Index: eigen3/Eigen/src/SparseCore/SparseSelfAdjointView.h
+===================================================================
+--- eigen3.orig/Eigen/src/SparseCore/SparseSelfAdjointView.h
++++ eigen3/Eigen/src/SparseCore/SparseSelfAdjointView.h
+@@ -211,8 +211,6 @@ struct evaluator_traits<SparseSelfAdjoin
+ {
+ typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
+ typedef SparseSelfAdjointShape Shape;
+-
+- static const int AssumeAliasing = 0;
+ };
+
+ struct SparseSelfAdjoint2Sparse {};
+Index: eigen3/Eigen/src/SparseQR/SparseQR.h
+===================================================================
+--- eigen3.orig/Eigen/src/SparseQR/SparseQR.h
++++ eigen3/Eigen/src/SparseQR/SparseQR.h
+@@ -691,7 +691,6 @@ struct evaluator_traits<SparseQRMatrixQR
+ typedef typename SparseQRType::MatrixType MatrixType;
+ typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
+ typedef SparseShape Shape;
+- static const int AssumeAliasing = 0;
+ };
+
+ template< typename DstXprType, typename SparseQRType>
+Index: eigen3/test/product.h
+===================================================================
+--- eigen3.orig/test/product.h
++++ eigen3/test/product.h
+@@ -145,14 +145,31 @@ template<typename MatrixType> void produ
+ 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);
++ }
+ }
+Index: eigen3/test/product_large.cpp
+===================================================================
+--- eigen3.orig/test/product_large.cpp
++++ eigen3/test/product_large.cpp
+@@ -9,6 +9,27 @@
+
+ #include "product.h"
+
++template<typename T>
++void test_aliasing()
++{
++ int rows = internal::random<int>(1,12);
++ int cols = internal::random<int>(1,12);
++ typedef Matrix<T,Dynamic,Dynamic> MatrixType;
++ typedef Matrix<T,Dynamic,1> 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<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
+ CALL_SUBTEST_4( product(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
+ CALL_SUBTEST_5( product(Matrix<float,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
++
++ CALL_SUBTEST_1( test_aliasing<float>() );
+ }
+
+ #if defined EIGEN_TEST_PART_6
+Index: eigen3/test/product_notemporary.cpp
+===================================================================
+--- eigen3.orig/test/product_notemporary.cpp
++++ eigen3/test/product_notemporary.cpp
+@@ -43,10 +43,16 @@ template<typename MatrixType> void produ
+ r1 = internal::random<Index>(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);
diff --git a/debian/patches/series b/debian/patches/series
index 14ae892..5ba5222 100644
--- a/debian/patches/series
+++ b/debian/patches/series
@@ -4,3 +4,4 @@
06_remove_doc_matrix.patch
07_remove_compressed_doc.patch
08_fix_path_FindEigen3.patch
+09_fix_1144.patch
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/eigen3.git
More information about the debian-science-commits
mailing list