* add a writable generic coeff wise expression (CwiseUnaryView)

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

6 files changed, 163 insertions, 6 deletions

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<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>
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 <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
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<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> >
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<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); }
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<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;
-    /** \internal the return type of MatrixBase::imag() */
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+                        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,11 +553,16 @@ 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>
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<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;