* add a writable generic coeff wise expression (CwiseUnaryView)

* add writable .real() and .imag() functions

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"

Eigen/src/Core/CwiseUnaryOp.h

@@ -165,14 +165,14 @@ MatrixBase<Derived>::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<typename Derived>
-EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::RealReturnType
+EIGEN_STRONG_INLINE typename MatrixBase<Derived>::RealReturnType
 MatrixBase<Derived>::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<typename Derived>

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 <g.gael@free.fr>
+//
+// 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 <http://www.gnu.org/licenses/>.
+
+#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<typename ViewOp, typename MatrixType>
+struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> >
+ : ei_traits<MatrixType>
+{
+  typedef typename ei_result_of<
+                     ViewOp(typename MatrixType::Scalar)
+                   >::type Scalar;
+  typedef typename MatrixType::Nested MatrixTypeNested;
+  typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
+  enum {
+    Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)),
+    CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<ViewOp>::Cost
+  };
+};
+
+template<typename ViewOp, typename MatrixType>
+class CwiseUnaryView : ei_no_assignment_operator,
+  public MatrixBase<CwiseUnaryView<ViewOp, MatrixType> >
+{
+  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<typename Derived>
+template<typename CustomViewOp>
+EIGEN_STRONG_INLINE const CwiseUnaryView<CustomViewOp, Derived>
+MatrixBase<Derived>::unaryViewExpr(const CustomViewOp& func) const
+{
+  return CwiseUnaryView<CustomViewOp, Derived>(derived(), func);
+}
+
+/** \returns a non const expression of the real part of \c *this.
+  *
+  * \sa imag() */
+template<typename Derived>
+EIGEN_STRONG_INLINE typename MatrixBase<Derived>::NonConstRealReturnType
+MatrixBase<Derived>::real() { return derived(); }
+
+/** \returns a non const expression of the imaginary part of \c *this.
+  *
+  * \sa real() */
+template<typename Derived>
+EIGEN_STRONG_INLINE typename MatrixBase<Derived>::NonConstImagReturnType
+MatrixBase<Derived>::imag() { return derived(); }
+
+#endif // EIGEN_CWISE_UNARY_VIEW_H

Eigen/src/Core/Functors.h

@@ -280,6 +280,7 @@ template<typename Scalar>
 struct ei_scalar_real_op EIGEN_EMPTY_STRUCT {
   typedef typename NumTraits<Scalar>::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<typename Scalar>
 struct ei_functor_traits<ei_scalar_real_op<Scalar> >
@@ -294,6 +295,7 @@ template<typename Scalar>
 struct ei_scalar_imag_op EIGEN_EMPTY_STRUCT {
   typedef typename NumTraits<Scalar>::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<typename Scalar>
 struct ei_functor_traits<ei_scalar_imag_op<Scalar> >

Eigen/src/Core/MathFunctions.h

@@ -53,6 +53,7 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y)
 template<> inline int precision<int>() { return 0; }
 template<> inline int machine_epsilon<int>() { 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<int>())
 template<> inline float precision<float>() { return 1e-5f; }
 template<> inline float machine_epsilon<float>() { 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<double>() { return 1e-11; }
 template<> inline double machine_epsilon<double>() { 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<std::complex<float> >() { return precision<flo
 template<> inline float machine_epsilon<std::complex<float> >() { return machine_epsilon<float>(); }
 inline float ei_real(const std::complex<float>& x) { return std::real(x); }
 inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
+inline float& ei_real_ref(std::complex<float>& x) { return reinterpret_cast<float*>(&x)[0]; }
+inline float& ei_imag_ref(std::complex<float>& x) { return reinterpret_cast<float*>(&x)[1]; }
 inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); }
 inline float ei_abs(const std::complex<float>& x) { return std::abs(x); }
 inline float ei_abs2(const std::complex<float>& x) { return std::norm(x); }
@@ -234,6 +239,8 @@ template<> inline double precision<std::complex<double> >() { return precision<d
 template<> inline double machine_epsilon<std::complex<double> >() { return machine_epsilon<double>(); }
 inline double ei_real(const std::complex<double>& x) { return std::real(x); }
 inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
+inline double& ei_real_ref(std::complex<double>& x) { return reinterpret_cast<double*>(&x)[0]; }
+inline double& ei_imag_ref(std::complex<double>& x) { return reinterpret_cast<double*>(&x)[1]; }
 inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); }
 inline double ei_abs(const std::complex<double>& x) { return std::abs(x); }
 inline double ei_abs2(const std::complex<double>& x) { return std::norm(x); }
@@ -268,6 +275,7 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double
 template<> inline long double precision<long double>() { return precision<double>(); }
 template<> inline long double machine_epsilon<long double>() { 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); }

Eigen/src/Core/MatrixBase.h

@@ -217,10 +217,20 @@ template<typename Derived> class MatrixBase
                         const CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>,
                         const Derived&
                      >::ret ConjugateReturnType;
+    /** \internal the return type of MatrixBase::real() const */
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+                        const CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived>,
+                        const Derived&
+                     >::ret RealReturnType;
     /** \internal the return type of MatrixBase::real() */
-    typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType;
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
-    /** \internal the return type of MatrixBase::imag() */
+                        CwiseUnaryView<ei_scalar_real_op<Scalar>, Derived>,
+                        Derived&
+                     >::ret NonConstRealReturnType;
+    /** \internal the return type of MatrixBase::imag() const */
     typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
+    /** \internal the return type of MatrixBase::imag() */
+    typedef CwiseUnaryView<ei_scalar_imag_op<Scalar>, Derived> NonConstImagReturnType;
     /** \internal the return type of MatrixBase::adjoint() */
     typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> >
             AdjointReturnType;
@@ -543,12 +553,17 @@ template<typename Derived> class MatrixBase
 
 
     ConjugateReturnType conjugate() const;
-    const RealReturnType real() const;
+    RealReturnType real() const;
+    NonConstRealReturnType real();
     const ImagReturnType imag() const;
+    NonConstImagReturnType imag();
 
     template<typename CustomUnaryOp>
     const CwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const;
     
+    template<typename CustomViewOp>
+    const CwiseUnaryView<CustomViewOp, Derived> unaryViewExpr(const CustomViewOp& func = CustomViewOp()) const;
+
     template<typename CustomBinaryOp, typename OtherDerived>
     const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
     binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const;

Eigen/src/Core/util/ForwardDeclarations.h

@@ -43,6 +43,7 @@ template<typename MatrixType> class Transpose;
 template<typename MatrixType> class Conjugate;
 template<typename NullaryOp, typename MatrixType>         class CwiseNullaryOp;
 template<typename UnaryOp,   typename MatrixType>         class CwiseUnaryOp;
+template<typename ViewOp,    typename MatrixType>         class CwiseUnaryView;
 template<typename BinaryOp,  typename Lhs, typename Rhs>  class CwiseBinaryOp;
 template<typename Lhs, typename Rhs, int ProductMode> class Product;
 template<typename CoeffsVectorType, typename Derived> class DiagonalMatrixBase;

test/basicstuff.cpp

@@ -107,6 +107,40 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   {
     VERIFY_IS_NOT_APPROX(m3, m1);
   }
+  
+  m3.real() = m1.real();
+  VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), static_cast<const MatrixType&>(m1).real());
+  VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), m1.real());
+}
+
+template<typename MatrixType> void basicStuffComplex(const MatrixType& m)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> RealMatrixType;
+
+  int rows = m.rows();
+  int cols = m.cols();
+  
+  Scalar s1 = ei_random<Scalar>(),
+         s2 = ei_random<Scalar>();
+  
+  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<const MatrixType&>(cm).real(), rm1);
+  VERIFY_IS_APPROX(static_cast<const MatrixType&>(cm).imag(), rm2);
+  cm.real().setZero();
+  VERIFY(static_cast<const MatrixType&>(cm).real().isZero());
+  VERIFY(!static_cast<const MatrixType&>(cm).imag().isZero());
 }
 
 void casting()
@@ -128,6 +162,9 @@ void test_basicstuff()
     CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) );
     CALL_SUBTEST( basicStuff(Matrix<float, 100, 100>()) );
     CALL_SUBTEST( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
+    
+    CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) );
+    CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) );
   }
 
   CALL_SUBTEST(casting());

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<float, 100, 100>()) );
+  
+  VERIFY(sizeof(std::complex<float>) == 2*sizeof(float));
+  VERIFY(sizeof(std::complex<double>) == 2*sizeof(double));
 }