From 239ada95b7680c75f793086eaa35fe7ec1047204 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 15 Aug 2009 22:19:29 +0200 Subject: [PATCH] add overloads of lazyAssign to detect common aliasing issue with transpose and adjoint --- Eigen/Core | 2 +- Eigen/src/Core/CwiseBinaryOp.h | 6 +++ Eigen/src/Core/MatrixBase.h | 15 ++++++ Eigen/src/Core/Transpose.h | 88 ++++++++++++++++++++++++++++++++++ test/adjoint.cpp | 14 +++++- 5 files changed, 122 insertions(+), 3 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index ff56c6ca2..854f737d6 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -150,6 +150,7 @@ namespace Eigen { #include "src/Core/Assign.h" #endif +#include "src/Core/util/BlasUtil.h" #include "src/Core/MatrixStorage.h" #include "src/Core/NestByValue.h" #include "src/Core/ReturnByValue.h" @@ -178,7 +179,6 @@ namespace Eigen { #include "src/Core/IO.h" #include "src/Core/Swap.h" #include "src/Core/CommaInitializer.h" -#include "src/Core/util/BlasUtil.h" #include "src/Core/ProductBase.h" #include "src/Core/Product.h" #include "src/Core/TriangularMatrix.h" diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index bcd4ae25f..78de9fbce 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -85,6 +85,8 @@ class CwiseBinaryOp : ei_no_assignment_operator, EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) typedef typename ei_traits::LhsNested LhsNested; typedef typename ei_traits::RhsNested RhsNested; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) : m_lhs(lhs), m_rhs(rhs), m_functor(func) @@ -130,6 +132,10 @@ class CwiseBinaryOp : ei_no_assignment_operator, return m_functor.packetOp(m_lhs.template packet(index), m_rhs.template packet(index)); } + const _LhsNested& lhs() const { return m_lhs; } + const _RhsNested& rhs() const { return m_rhs; } + const BinaryOp& functor() const { return m_functor; } + protected: const LhsNested m_lhs; const RhsNested m_rhs; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index a554de8e3..6ed96e286 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -457,6 +457,21 @@ template class MatrixBase void transposeInPlace(); const AdjointReturnType adjoint() const; void adjointInPlace(); + #ifndef EIGEN_NO_DEBUG + template + Derived& lazyAssign(const Transpose& other); + template + Derived& lazyAssign(const CwiseBinaryOp,Transpose,DerivedB>& other); + template + Derived& lazyAssign(const CwiseBinaryOp,DerivedA,Transpose >& other); + + template + Derived& lazyAssign(const CwiseUnaryOp, NestByValue > >& other); + template + Derived& lazyAssign(const CwiseBinaryOp,CwiseUnaryOp, NestByValue > >,DerivedB>& other); + template + Derived& lazyAssign(const CwiseBinaryOp,DerivedA,CwiseUnaryOp, NestByValue > > >& other); + #endif RowXpr row(int i); const RowXpr row(int i) const; diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 0787daa61..8527edc2a 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -267,4 +267,92 @@ inline void MatrixBase::adjointInPlace() derived() = adjoint().eval(); } +#ifndef EIGEN_NO_DEBUG + +// The following is to detect aliasing problems in the following common cases: +// a = a.transpose() +// a = a.transpose() + X +// a = X + a.transpose() +// a = a.adjoint() +// a = a.adjoint() + X +// a = X + a.adjoint() + +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 +template +Derived& MatrixBase::lazyAssign(const Transpose& other) +{ + ei_assert(ei_extract_data(other) != ei_extract_data(derived()) + && "aliasing detected during tranposition, please use transposeInPlace()"); + return lazyAssign(static_cast >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,Transpose,DerivedB>& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,Transpose,DerivedB> >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,DerivedA,Transpose >& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,DerivedA,Transpose > >& >(other)); +} + +template +template Derived& +MatrixBase:: +lazyAssign(const CwiseUnaryOp, NestByValue > >& other) +{ + ei_assert(ei_extract_data(other) != ei_extract_data(derived()) + && "aliasing detected during tranposition, please use adjointInPlace()"); + return lazyAssign(static_cast, NestByValue > > >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,CwiseUnaryOp, NestByValue > >,DerivedB>& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,CwiseUnaryOp, NestByValue > >,DerivedB> >& >(other)); +} + +template +template +Derived& MatrixBase:: +lazyAssign(const CwiseBinaryOp,DerivedA,CwiseUnaryOp, NestByValue > > >& other) +{ + ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) + && "aliasing detected during tranposition, please evaluate your expression"); + return lazyAssign(static_cast,DerivedA,CwiseUnaryOp, NestByValue > > > >& >(other)); +} +#endif + #endif // EIGEN_TRANSPOSE_H diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 300650338..964658c65 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -21,7 +21,7 @@ // You should have received a copy of the GNU Lesser General Public // License and a copy of the GNU General Public License along with // Eigen. If not, see . -#define EIGEN_NO_ASSERTION_CHECKING + #include "main.h" template void adjoint(const MatrixType& m) @@ -110,7 +110,6 @@ template void adjoint(const MatrixType& m) m3.transposeInPlace(); VERIFY_IS_APPROX(m3,m1.conjugate()); - } void test_adjoint() @@ -125,5 +124,16 @@ void test_adjoint() } // test a large matrix only once CALL_SUBTEST( adjoint(Matrix()) ); + + { + MatrixXcf a(10,10), b(10,10); + VERIFY_RAISES_ASSERT(a = a.transpose()); + VERIFY_RAISES_ASSERT(a = a.transpose() + b); + VERIFY_RAISES_ASSERT(a = b + a.transpose()); + VERIFY_RAISES_ASSERT(a = a.conjugate().transpose()); + VERIFY_RAISES_ASSERT(a = a.adjoint()); + VERIFY_RAISES_ASSERT(a = a.adjoint() + b); + VERIFY_RAISES_ASSERT(a = b + a.adjoint()); + } }