From dd45c4805ced4ad8ead743875f42259e9b6a2795 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 20 May 2009 15:41:23 +0200 Subject: [PATCH] * add a writable generic coeff wise expression (CwiseUnaryView) * add writable .real() and .imag() functions --- Eigen/Core | 1 + Eigen/src/Core/CwiseUnaryOp.h | 6 +- Eigen/src/Core/CwiseUnaryView.h | 131 ++++++++++++++++++++++ Eigen/src/Core/Functors.h | 2 + Eigen/src/Core/MathFunctions.h | 8 ++ Eigen/src/Core/MatrixBase.h | 21 +++- Eigen/src/Core/util/ForwardDeclarations.h | 1 + test/basicstuff.cpp | 37 ++++++ test/sizeof.cpp | 3 + 9 files changed, 204 insertions(+), 6 deletions(-) create mode 100644 Eigen/src/Core/CwiseUnaryView.h diff --git a/Eigen/Core b/Eigen/Core index 810daf81f..47121fda2 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -157,6 +157,7 @@ namespace Eigen { #include "src/Core/CwiseBinaryOp.h" #include "src/Core/CwiseUnaryOp.h" #include "src/Core/CwiseNullaryOp.h" +#include "src/Core/CwiseUnaryView.h" #include "src/Core/Dot.h" #include "src/Core/DiagonalProduct.h" #include "src/Core/SolveTriangular.h" diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 7289952ed..7507f1ceb 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -165,14 +165,14 @@ MatrixBase::conjugate() const return ConjugateReturnType(derived()); } -/** \returns an expression of the real part of \c *this. +/** \returns a read-only expression of the real part of \c *this. * * \sa imag() */ template -EIGEN_STRONG_INLINE const typename MatrixBase::RealReturnType +EIGEN_STRONG_INLINE typename MatrixBase::RealReturnType MatrixBase::real() const { return derived(); } -/** \returns an expression of the imaginary part of \c *this. +/** \returns an read-only expression of the imaginary part of \c *this. * * \sa real() */ template diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h new file mode 100644 index 000000000..043bc2d67 --- /dev/null +++ b/Eigen/src/Core/CwiseUnaryView.h @@ -0,0 +1,131 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// 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 . + +#ifndef EIGEN_CWISE_UNARY_VIEW_H +#define EIGEN_CWISE_UNARY_VIEW_H + +/** \class CwiseUnaryView + * + * \brief Generic lvalue expression of a coefficient-wise unary operator of a matrix or a vector + * + * \param ViewOp template functor implementing the view + * \param MatrixType the type of the matrix we are applying the unary operator + * + * This class represents a lvalue expression of a generic unary view operator of a matrix or a vector. + * It is the return type of real() and imag(), and most of the time this is the only way it is used. + * + * \sa MatrixBase::unaryViewExpr(const CustomUnaryOp &) const, class CwiseUnaryOp + */ +template +struct ei_traits > + : ei_traits +{ + typedef typename ei_result_of< + ViewOp(typename MatrixType::Scalar) + >::type Scalar; + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename ei_unref::type _MatrixTypeNested; + enum { + Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits::Cost + }; +}; + +template +class CwiseUnaryView : ei_no_assignment_operator, + public MatrixBase > +{ + public: + + EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseUnaryView) + + inline CwiseUnaryView(const MatrixType& mat, const ViewOp& func = ViewOp()) + : m_matrix(mat), m_functor(func) {} + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryView) + + EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_matrix.cols(); } + + EIGEN_STRONG_INLINE const Scalar coeff(int row, int col) const + { + return m_functor(m_matrix.coeff(row, col)); + } + + EIGEN_STRONG_INLINE const Scalar coeff(int index) const + { + return m_functor(m_matrix.coeff(index)); + } + + EIGEN_STRONG_INLINE Scalar& coeffRef(int row, int col) + { + return m_functor(m_matrix.const_cast_derived().coeffRef(row, col)); + } + + EIGEN_STRONG_INLINE Scalar& coeffRef(int index) + { + return m_functor(m_matrix.const_cast_derived().coeffRef(index)); + } + + protected: + const typename MatrixType::Nested m_matrix; + const ViewOp m_functor; +}; + +/** \returns an expression of a custom coefficient-wise unary operator \a func of *this + * + * The template parameter \a CustomUnaryOp is the type of the functor + * of the custom unary operator. + * + * \addexample CustomCwiseUnaryFunctors \label How to use custom coeff wise unary functors + * + * Example: + * \include class_CwiseUnaryOp.cpp + * Output: \verbinclude class_CwiseUnaryOp.out + * + * \sa class CwiseUnaryOp, class CwiseBinarOp, MatrixBase::operator-, Cwise::abs + */ +template +template +EIGEN_STRONG_INLINE const CwiseUnaryView +MatrixBase::unaryViewExpr(const CustomViewOp& func) const +{ + return CwiseUnaryView(derived(), func); +} + +/** \returns a non const expression of the real part of \c *this. + * + * \sa imag() */ +template +EIGEN_STRONG_INLINE typename MatrixBase::NonConstRealReturnType +MatrixBase::real() { return derived(); } + +/** \returns a non const expression of the imaginary part of \c *this. + * + * \sa real() */ +template +EIGEN_STRONG_INLINE typename MatrixBase::NonConstImagReturnType +MatrixBase::imag() { return derived(); } + +#endif // EIGEN_CWISE_UNARY_VIEW_H diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 38371cf2a..631307517 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -280,6 +280,7 @@ template struct ei_scalar_real_op EIGEN_EMPTY_STRUCT { typedef typename NumTraits::Real result_type; EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return ei_real(a); } + EIGEN_STRONG_INLINE result_type& operator() (Scalar& a) const { return ei_real_ref(a); } }; template struct ei_functor_traits > @@ -294,6 +295,7 @@ template struct ei_scalar_imag_op EIGEN_EMPTY_STRUCT { typedef typename NumTraits::Real result_type; EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return ei_imag(a); } + EIGEN_STRONG_INLINE result_type& operator() (Scalar& a) const { return ei_imag_ref(a); } }; template struct ei_functor_traits > diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index e201f98b2..64131ea78 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -53,6 +53,7 @@ template inline typename NumTraits::Real ei_hypot(T x, T y) template<> inline int precision() { return 0; } template<> inline int machine_epsilon() { return 0; } inline int ei_real(int x) { return x; } +inline int& ei_real_ref(int& x) { return x; } inline int ei_imag(int) { return 0; } inline int ei_conj(int x) { return x; } inline int ei_abs(int x) { return abs(x); } @@ -106,6 +107,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision()) template<> inline float precision() { return 1e-5f; } template<> inline float machine_epsilon() { return 1.192e-07f; } inline float ei_real(float x) { return x; } +inline float& ei_real_ref(float& x) { return x; } inline float ei_imag(float) { return 0.f; } inline float ei_conj(float x) { return x; } inline float ei_abs(float x) { return std::abs(x); } @@ -153,6 +155,7 @@ template<> inline double precision() { return 1e-11; } template<> inline double machine_epsilon() { return 2.220e-16; } inline double ei_real(double x) { return x; } +inline double& ei_real_ref(double& x) { return x; } inline double ei_imag(double) { return 0.; } inline double ei_conj(double x) { return x; } inline double ei_abs(double x) { return std::abs(x); } @@ -200,6 +203,8 @@ template<> inline float precision >() { return precision inline float machine_epsilon >() { return machine_epsilon(); } inline float ei_real(const std::complex& x) { return std::real(x); } inline float ei_imag(const std::complex& x) { return std::imag(x); } +inline float& ei_real_ref(std::complex& x) { return reinterpret_cast(&x)[0]; } +inline float& ei_imag_ref(std::complex& x) { return reinterpret_cast(&x)[1]; } inline std::complex ei_conj(const std::complex& x) { return std::conj(x); } inline float ei_abs(const std::complex& x) { return std::abs(x); } inline float ei_abs2(const std::complex& x) { return std::norm(x); } @@ -234,6 +239,8 @@ template<> inline double precision >() { return precision inline double machine_epsilon >() { return machine_epsilon(); } inline double ei_real(const std::complex& x) { return std::real(x); } inline double ei_imag(const std::complex& x) { return std::imag(x); } +inline double& ei_real_ref(std::complex& x) { return reinterpret_cast(&x)[0]; } +inline double& ei_imag_ref(std::complex& x) { return reinterpret_cast(&x)[1]; } inline std::complex ei_conj(const std::complex& x) { return std::conj(x); } inline double ei_abs(const std::complex& x) { return std::abs(x); } inline double ei_abs2(const std::complex& x) { return std::norm(x); } @@ -268,6 +275,7 @@ inline bool ei_isApprox(const std::complex& a, const std::complex inline long double precision() { return precision(); } template<> inline long double machine_epsilon() { return 1.084e-19l; } inline long double ei_real(long double x) { return x; } +inline long double& ei_real_ref(long double& x) { return x; } inline long double ei_imag(long double) { return 0.; } inline long double ei_conj(long double x) { return x; } inline long double ei_abs(long double x) { return std::abs(x); } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 7df9f307b..4e689f614 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -217,10 +217,20 @@ template class MatrixBase const CwiseUnaryOp, Derived>, const Derived& >::ret ConjugateReturnType; + /** \internal the return type of MatrixBase::real() const */ + typedef typename ei_meta_if::IsComplex, + const CwiseUnaryOp, Derived>, + const Derived& + >::ret RealReturnType; /** \internal the return type of MatrixBase::real() */ - typedef CwiseUnaryOp, Derived> RealReturnType; - /** \internal the return type of MatrixBase::imag() */ + typedef typename ei_meta_if::IsComplex, + CwiseUnaryView, Derived>, + Derived& + >::ret NonConstRealReturnType; + /** \internal the return type of MatrixBase::imag() const */ typedef CwiseUnaryOp, Derived> ImagReturnType; + /** \internal the return type of MatrixBase::imag() */ + typedef CwiseUnaryView, Derived> NonConstImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ typedef Eigen::Transpose::type> > AdjointReturnType; @@ -543,11 +553,16 @@ template class MatrixBase ConjugateReturnType conjugate() const; - const RealReturnType real() const; + RealReturnType real() const; + NonConstRealReturnType real(); const ImagReturnType imag() const; + NonConstImagReturnType imag(); template const CwiseUnaryOp unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const; + + template + const CwiseUnaryView unaryViewExpr(const CustomViewOp& func = CustomViewOp()) const; template const CwiseBinaryOp diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 9ef708194..070c6f38e 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -43,6 +43,7 @@ template class Transpose; template class Conjugate; template class CwiseNullaryOp; template class CwiseUnaryOp; +template class CwiseUnaryView; template class CwiseBinaryOp; template class Product; template class DiagonalMatrixBase; diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 9bc1b7459..7b22cfb79 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -107,6 +107,40 @@ template void basicStuff(const MatrixType& m) { VERIFY_IS_NOT_APPROX(m3, m1); } + + m3.real() = m1.real(); + VERIFY_IS_APPROX(static_cast(m3).real(), static_cast(m1).real()); + VERIFY_IS_APPROX(static_cast(m3).real(), m1.real()); +} + +template void basicStuffComplex(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix RealMatrixType; + + int rows = m.rows(); + int cols = m.cols(); + + Scalar s1 = ei_random(), + s2 = ei_random(); + + VERIFY(ei_real(s1)==ei_real_ref(s1)); + VERIFY(ei_imag(s1)==ei_imag_ref(s1)); + ei_real_ref(s1) = ei_real(s2); + ei_imag_ref(s1) = ei_imag(s2); + VERIFY(s1==s2); + + RealMatrixType rm1 = RealMatrixType::Random(rows,cols), + rm2 = RealMatrixType::Random(rows,cols); + MatrixType cm(rows,cols); + cm.real() = rm1; + cm.imag() = rm2; + VERIFY_IS_APPROX(static_cast(cm).real(), rm1); + VERIFY_IS_APPROX(static_cast(cm).imag(), rm2); + cm.real().setZero(); + VERIFY(static_cast(cm).real().isZero()); + VERIFY(!static_cast(cm).imag().isZero()); } void casting() @@ -128,6 +162,9 @@ void test_basicstuff() CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) ); CALL_SUBTEST( basicStuff(Matrix()) ); CALL_SUBTEST( basicStuff(Matrix(10,10)) ); + + CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) ); + CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) ); } CALL_SUBTEST(casting()); diff --git a/test/sizeof.cpp b/test/sizeof.cpp index 6c0485ae1..c713afade 100644 --- a/test/sizeof.cpp +++ b/test/sizeof.cpp @@ -43,4 +43,7 @@ void test_sizeof() CALL_SUBTEST( verifySizeOf(MatrixXi(8, 12)) ); CALL_SUBTEST( verifySizeOf(MatrixXcd(20, 20)) ); CALL_SUBTEST( verifySizeOf(Matrix()) ); + + VERIFY(sizeof(std::complex) == 2*sizeof(float)); + VERIFY(sizeof(std::complex) == 2*sizeof(double)); }