the big Array/Cwise rework as discussed on the mailing list. The new API

can be seen in Eigen/src/Core/Cwise.h.

26 files changed, 465 insertions, 389 deletions

diff --git a/Eigen/Array b/Eigen/Array
index 33ae6d404..51c3abe31 100644
--- a/Eigen/Array
+++ b/Eigen/Array
@@ -10,7 +10,6 @@ namespace Eigen {
 #include "src/Array/AllAndAny.h"
 #include "src/Array/PartialRedux.h"
 #include "src/Array/Random.h"
-#include "src/Array/Array.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/Core b/Eigen/Core
index b9c7e6c0f..ecfd141b2 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -36,6 +36,7 @@ namespace Eigen {
 #include "src/Core/NestByValue.h"
 #include "src/Core/Flagged.h"
 #include "src/Core/Matrix.h"
+#include "src/Core/Cwise.h"
 #include "src/Core/CwiseBinaryOp.h"
 #include "src/Core/CwiseUnaryOp.h"
 #include "src/Core/CwiseNullaryOp.h"
diff --git a/Eigen/src/Array/Array.h b/Eigen/src/Array/Array.h
deleted file mode 100644
index 5fd074e76..000000000
--- a/Eigen/src/Array/Array.h
+++ /dev/null
@@ -1,152 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
-//
-// Copyright (C) 2008 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_ARRAY_H
-#define EIGEN_ARRAY_H
-
-/** \class Array
-  *
-  * \brief Pseudo expression offering additional features to an expression
-  *
-  * \param ExpressionType the type of the object of which we want array related features
-  *
-  * This class represents an expression with additional, array related, features.
-  * It is the return type of MatrixBase::array()
-  * and most of the time this is the only way it is used.
-  *
-  * \sa MatrixBase::array()
-  */
-template<typename ExpressionType> class Array
-{
-  public:
-
-    typedef typename ei_traits<ExpressionType>::Scalar Scalar;
-    typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
-        ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
-//     typedef NestByValue<typename ExpressionType::ConstantReturnType> ConstantReturnType;
-    typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType;
-
-    inline Array(const ExpressionType& matrix) : m_matrix(matrix) {}
-
-    /** \internal */
-    inline const ExpressionType& _expression() const { return m_matrix; }
-
-    const ScalarAddReturnType
-    operator+(const Scalar& scalar) const;
-
-    /** \relates Array */
-    friend const ScalarAddReturnType
-    operator+(const Scalar& scalar, const Array& mat)
-    { return mat + scalar; }
-
-    ExpressionType& operator+=(const Scalar& scalar);
-
-    const ScalarAddReturnType
-    operator-(const Scalar& scalar) const;
-
-    ExpressionType& operator-=(const Scalar& scalar);
-
-    /** \returns true if each coeff of \c *this is less than its respective coeff of \a other */
-    template<typename OtherDerived> bool operator<(const Array<OtherDerived>& other) const
-    { return m_matrix.cwiseLessThan(other._expression()).all(); }
-
-    /** \returns true if each coeff of \c *this is less or equal to its respective coeff of \a other */
-    template<typename OtherDerived> bool operator<=(const Array<OtherDerived>& other) const
-    { return m_matrix.cwiseLessEqual(other._expression()).all(); }
-
-    /** \returns true if each coeff of \c *this is greater to its respective coeff of \a other */
-    template<typename OtherDerived>
-    bool operator>(const Array<OtherDerived>& other) const
-    { return m_matrix.cwiseGreaterThan(other._expression()).all(); }
-
-    /** \returns true if each coeff of \c *this is greater or equal to its respective coeff of \a other */
-    template<typename OtherDerived>
-    bool operator>=(const Array<OtherDerived>& other) const
-    { return m_matrix.cwiseGreaterEqual(other._expression()).all(); }
-
-  protected:
-    ExpressionTypeNested m_matrix;
-};
-
-/** \returns an expression of \c *this with each coeff incremented by the constant \a scalar */
-template<typename ExpressionType>
-const typename Array<ExpressionType>::ScalarAddReturnType
-Array<ExpressionType>::operator+(const Scalar& scalar) const
-{
-  return CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType>(m_matrix, ei_scalar_add_op<Scalar>(scalar));
-}
-
-/** \see operator+ */
-template<typename ExpressionType>
-ExpressionType& Array<ExpressionType>::operator+=(const Scalar& scalar)
-{
-  m_matrix.const_cast_derived() = *this + scalar;
-  return m_matrix.const_cast_derived();
-}
-
-/** \returns an expression of \c *this with each coeff decremented by the constant \a scalar */
-template<typename ExpressionType>
-const typename Array<ExpressionType>::ScalarAddReturnType
-Array<ExpressionType>::operator-(const Scalar& scalar) const
-{
-  return *this + (-scalar);
-}
-
-/** \see operator- */
-template<typename ExpressionType>
-ExpressionType& Array<ExpressionType>::operator-=(const Scalar& scalar)
-{
-  m_matrix.const_cast_derived() = *this - scalar;
-  return m_matrix.const_cast_derived();
-}
-
-/** \array_module
-  *
-  * \returns an Array expression of *this providing additional,
-  * array related, features.
-  *
-  * \sa class Array
-  */
-template<typename Derived>
-inline const Array<Derived>
-MatrixBase<Derived>::array() const
-{
-  return derived();
-}
-
-/** \array_module
-  *
-  * \returns an Array expression of *this providing additional,
-  * array related, features.
-  *
-  * \sa class Array
-  */
-template<typename Derived>
-inline Array<Derived>
-MatrixBase<Derived>::array()
-{
-  return derived();
-}
-
-#endif // EIGEN_FLAGGED_H
diff --git a/Eigen/src/Array/CwiseOperators.h b/Eigen/src/Array/CwiseOperators.h
index d1f34b0e3..7481b4c0b 100644
--- a/Eigen/src/Array/CwiseOperators.h
+++ b/Eigen/src/Array/CwiseOperators.h
@@ -30,74 +30,97 @@
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise square root of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_sqrt_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseSqrt() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_sqrt_op>::ReturnType
+Cwise<ExpressionType>::sqrt() const
 {
-  return derived();
+  return _expression();
 }
 
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise exponential of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_exp_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseExp() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_exp_op>::ReturnType
+Cwise<ExpressionType>::exp() const
 {
-  return derived();
+  return _expression();
 }
 
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise logarithm of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_log_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseLog() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_log_op>::ReturnType
+Cwise<ExpressionType>::log() const
 {
-  return derived();
+  return _expression();
 }
 
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise cosine of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_cos_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseCos() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_cos_op>::ReturnType
+Cwise<ExpressionType>::cos() const
 {
-  return derived();
+  return _expression();
 }
 
+
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise sine of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_sin_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseSin() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_sin_op>::ReturnType
+Cwise<ExpressionType>::sin() const
 {
-  return derived();
+  return _expression();
 }
 
+
 /** \array_module
   * 
   * \returns an expression of the coefficient-wise power of *this to the given exponent. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_pow_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwisePow(const Scalar& exponent) const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_pow_op>::ReturnType
+Cwise<ExpressionType>::pow(const Scalar& exponent) const
 {
-  return CwiseUnaryOp<ei_scalar_pow_op<Scalar>, Derived>
-    (derived(), ei_scalar_pow_op<Scalar>(exponent));
+  return typename UnOp<ei_scalar_pow_op>::ReturnType(_expression(), ei_scalar_pow_op<Scalar>(exponent));
 }
 
+
 /** \array_module
   * 
-  * \returns an expression of the coefficient-wise reciprocal of *this. */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_inverse_op<typename ei_traits<Derived>::Scalar>, Derived>
-MatrixBase<Derived>::cwiseInverse() const
+  * \returns an expression of the coefficient-wise inverse of *this. */
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_inverse_op>::ReturnType
+Cwise<ExpressionType>::inverse() const
+{
+  return _expression();
+}
+
+/** \array_module
+  *
+  * \returns an expression of the coefficient-wise square of *this. */
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_square_op>::ReturnType
+Cwise<ExpressionType>::square() const
 {
-  return derived();
+  return _expression();
 }
 
+/** \array_module
+  *
+  * \returns an expression of the coefficient-wise cube of *this. */
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_cube_op>::ReturnType
+Cwise<ExpressionType>::cube() const
+{
+  return _expression();
+}
+
+
 // -- binary operators --
 
 /** \array_module
@@ -106,12 +129,12 @@ MatrixBase<Derived>::cwiseInverse() const
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::less<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseLessThan(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::less, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator<(const MatrixBase<OtherDerived> &other) const
 {
-  return cwise(other, std::less<Scalar>());
+  return typename BinOp<std::less, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \array_module
@@ -120,12 +143,12 @@ MatrixBase<Derived>::cwiseLessThan(const MatrixBase<OtherDerived> &other) const
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::less_equal<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseLessEqual(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::less_equal, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator<=(const MatrixBase<OtherDerived> &other) const
 {
-  return cwise(other, std::less_equal<Scalar>());
+  return typename BinOp<std::less_equal, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \array_module
@@ -134,12 +157,12 @@ MatrixBase<Derived>::cwiseLessEqual(const MatrixBase<OtherDerived> &other) const
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::greater<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseGreaterThan(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::greater, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator>(const MatrixBase<OtherDerived> &other) const
 {
-  return cwise(other, std::greater<Scalar>());
+  return typename BinOp<std::greater, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \array_module
@@ -148,12 +171,12 @@ MatrixBase<Derived>::cwiseGreaterThan(const MatrixBase<OtherDerived> &other) con
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::greater_equal<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseGreaterEqual(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::greater_equal, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator>=(const MatrixBase<OtherDerived> &other) const
 {
-  return cwise(other, std::greater_equal<Scalar>());
+  return typename BinOp<std::greater_equal, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \array_module
@@ -162,12 +185,12 @@ MatrixBase<Derived>::cwiseGreaterEqual(const MatrixBase<OtherDerived> &other) co
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::equal_to<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseEqualTo(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::equal_to, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator==(const MatrixBase<OtherDerived> &other) const
 {
-  return cwise(other, std::equal_to<Scalar>());
+  return typename BinOp<std::equal_to, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \array_module
@@ -176,12 +199,43 @@ MatrixBase<Derived>::cwiseEqualTo(const MatrixBase<OtherDerived> &other) const
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<std::not_equal_to<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseNotEqualTo(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<std::not_equal_to, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator!=(const MatrixBase<OtherDerived> &other) const
+{
+  return typename BinOp<std::not_equal_to, OtherDerived>::ReturnType(_expression(), other.derived());
+}
+
+
+/** \returns an expression of \c *this with each coeff incremented by the constant \a scalar */
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::ScalarAddReturnType
+Cwise<ExpressionType>::operator+(const Scalar& scalar) const
+{
+  return typename Cwise<ExpressionType>::ScalarAddReturnType(m_matrix, ei_scalar_add_op<Scalar>(scalar));
+}
+
+/** \see operator+ */
+template<typename ExpressionType>
+inline ExpressionType& Cwise<ExpressionType>::operator+=(const Scalar& scalar)
+{
+  return m_matrix.const_cast_derived() = *this + scalar;
+}
+
+/** \returns an expression of \c *this with each coeff decremented by the constant \a scalar */
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::ScalarAddReturnType
+Cwise<ExpressionType>::operator-(const Scalar& scalar) const
+{
+  return *this + (-scalar);
+}
+
+/** \see operator- */
+template<typename ExpressionType>
+inline ExpressionType& Cwise<ExpressionType>::operator-=(const Scalar& scalar)
 {
-  return cwise(other, std::not_equal_to<Scalar>());
+  return m_matrix.const_cast_derived() = *this - scalar;
 }
 
 #endif // EIGEN_ARRAY_CWISE_OPERATORS_H
diff --git a/Eigen/src/Array/Functors.h b/Eigen/src/Array/Functors.h
index 03ae98090..2da93d8e9 100644
--- a/Eigen/src/Array/Functors.h
+++ b/Eigen/src/Array/Functors.h
@@ -25,15 +25,6 @@
 #ifndef EIGEN_ARRAY_FUNCTORS_H
 #define EIGEN_ARRAY_FUNCTORS_H
 
-/** \internal
-  * \array_module
-  *
-  * \brief Template functor to add a scalar to a fixed other one
-  *
-  * \sa class CwiseUnaryOp, Array::operator+
-  */
-template<typename Scalar, bool PacketAccess = (int(ei_packet_traits<Scalar>::size)>1?true:false) > struct ei_scalar_add_op;
-
 template<typename Scalar>
 struct ei_scalar_add_op<Scalar,true> {
   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
@@ -150,17 +141,59 @@ struct ei_functor_traits<ei_scalar_pow_op<Scalar> >
   *
   * \array_module
   *
-  * \brief Template functor to compute the reciprocal of a scalar
+  * \brief Template functor to compute the inverse of a scalar
   *
-  * \sa class CwiseUnaryOp, MatrixBase::cwiseInverse
+  * \sa class CwiseUnaryOp, Cwise::inverse()
   */
 template<typename Scalar>
 struct ei_scalar_inverse_op {
   inline Scalar operator() (const Scalar& a) const { return Scalar(1)/a; }
+  template<typename PacketScalar>
+  inline const PacketScalar packetOp(const PacketScalar& a) const
+  { return ei_div(ei_pset1(Scalar(1)),a); }
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_inverse_op<Scalar> >
-{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = false }; };
+{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = int(ei_packet_traits<Scalar>::size)>1 }; };
+
+/** \internal
+  *
+  * \array_module
+  *
+  * \brief Template functor to compute the square of a scalar
+  *
+  * \sa class CwiseUnaryOp, Cwise::square()
+  */
+template<typename Scalar>
+struct ei_scalar_square_op {
+  inline Scalar operator() (const Scalar& a) const { return a*a; }
+  template<typename PacketScalar>
+  inline const PacketScalar packetOp(const PacketScalar& a) const
+  { return ei_pmul(a,a); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_square_op<Scalar> >
+{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = int(ei_packet_traits<Scalar>::size)>1 }; };
+
+/** \internal
+  *
+  * \array_module
+  *
+  * \brief Template functor to compute the cube of a scalar
+  *
+  * \sa class CwiseUnaryOp, Cwise::cube()
+  */
+template<typename Scalar>
+struct ei_scalar_cube_op {
+  inline Scalar operator() (const Scalar& a) const { return a*a*a; }
+  template<typename PacketScalar>
+  inline const PacketScalar packetOp(const PacketScalar& a) const
+  { return ei_pmul(a,ei_pmul(a,a)); }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_cube_op<Scalar> >
+{ enum { Cost = 2*NumTraits<Scalar>::MulCost, PacketAccess = int(ei_packet_traits<Scalar>::size)>1 }; };
+
 
 // default ei_functor_traits for STL functors:
 
diff --git a/Eigen/src/Array/Random.h b/Eigen/src/Array/Random.h
index e4b61af33..5cafb878a 100644
--- a/Eigen/src/Array/Random.h
+++ b/Eigen/src/Array/Random.h
@@ -53,7 +53,7 @@ template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
 MatrixBase<Derived>::random(int rows, int cols)
 {
-  return create(rows, cols, ei_scalar_random_op<Scalar>());
+  return NullaryExpr(rows, cols, ei_scalar_random_op<Scalar>());
 }
 
 /** \array_module
@@ -78,7 +78,7 @@ template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
 MatrixBase<Derived>::random(int size)
 {
-  return create(size, ei_scalar_random_op<Scalar>());
+  return NullaryExpr(size, ei_scalar_random_op<Scalar>());
 }
 
 /** \array_module
@@ -98,7 +98,7 @@ template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_random_op<typename ei_traits<Derived>::Scalar>, Derived>
 MatrixBase<Derived>::random()
 {
-  return create(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_random_op<Scalar>());
+  return NullaryExpr(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_random_op<Scalar>());
 }
 
 /** \array_module
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index 652cf670a..8905385cc 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -136,7 +136,7 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBase<D
     .inverseProduct(
       (matrixL()
         .inverseProduct(vecB))
-        .cwiseQuotient(m_matrix.diagonal())
+        .cwise()/m_matrix.diagonal()
       );
 }
 
diff --git a/Eigen/src/Core/Cwise.h b/Eigen/src/Core/Cwise.h
new file mode 100644
index 000000000..63b9834c5
--- /dev/null
+++ b/Eigen/src/Core/Cwise.h
@@ -0,0 +1,173 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.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_H
+#define EIGEN_CWISE_H
+
+/** \internal
+  * \array_module
+  *
+  * \brief Template functor to add a scalar to a fixed other one
+  *
+  * \sa class CwiseUnaryOp, Array::operator+
+  */
+template<typename Scalar, bool PacketAccess = (int(ei_packet_traits<Scalar>::size)>1?true:false) > struct ei_scalar_add_op;
+
+/** \class Cwise
+  *
+  * \brief Pseudo expression providing additional coefficient-wise operations
+  *
+  * \param ExpressionType the type of the object on which to do coefficient-wise operations
+  *
+  * This class represents an expression with additional coefficient-wise features.
+  * It is the return type of MatrixBase::cwise()
+  * and most of the time this is the only way it is used.
+  *
+  * \sa MatrixBase::cwise()
+  */
+template<typename ExpressionType> class Cwise
+{
+  public:
+
+    typedef typename ei_traits<ExpressionType>::Scalar Scalar;
+    typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
+        ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
+//     typedef NestByValue<typename ExpressionType::ConstantReturnType> ConstantReturnType;
+    typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType;
+
+    template<template<typename _Scalar> class Functor, typename OtherDerived> struct BinOp
+    {
+      typedef CwiseBinaryOp<Functor<typename ei_traits<ExpressionType>::Scalar>,
+                            ExpressionType,
+                            OtherDerived
+                           > ReturnType;
+    };
+
+    template<template<typename _Scalar> class Functor> struct UnOp
+    {
+      typedef CwiseUnaryOp<Functor<typename ei_traits<ExpressionType>::Scalar>,
+                           ExpressionType
+                          > ReturnType;
+    };
+
+    inline Cwise(const ExpressionType& matrix) : m_matrix(matrix) {}
+
+    /** \internal */
+    inline const ExpressionType& _expression() const { return m_matrix; }
+
+    template<typename OtherDerived>
+    const typename BinOp<ei_scalar_product_op, OtherDerived>::ReturnType
+    operator*(const MatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const typename BinOp<ei_scalar_quotient_op, OtherDerived>::ReturnType
+    operator/(const MatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const typename BinOp<ei_scalar_min_op, OtherDerived>::ReturnType
+    min(const MatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const typename BinOp<ei_scalar_max_op, OtherDerived>::ReturnType
+    max(const MatrixBase<OtherDerived> &other) const;
+
+    const typename UnOp<ei_scalar_abs_op>::ReturnType abs() const;
+    const typename UnOp<ei_scalar_abs2_op>::ReturnType abs2() const;
+    const typename UnOp<ei_scalar_square_op>::ReturnType square() const;
+    const typename UnOp<ei_scalar_cube_op>::ReturnType cube() const;
+    const typename UnOp<ei_scalar_inverse_op>::ReturnType inverse() const;
+    const typename UnOp<ei_scalar_sqrt_op>::ReturnType sqrt() const;
+    const typename UnOp<ei_scalar_exp_op>::ReturnType exp() const;
+    const typename UnOp<ei_scalar_log_op>::ReturnType log() const;
+    const typename UnOp<ei_scalar_cos_op>::ReturnType cos() const;
+    const typename UnOp<ei_scalar_sin_op>::ReturnType sin() const;
+    const typename UnOp<ei_scalar_pow_op>::ReturnType pow(const Scalar& exponent) const;
+
+
+    const ScalarAddReturnType
+    operator+(const Scalar& scalar) const;
+
+    /** \relates Cwise */
+    friend const ScalarAddReturnType
+    operator+(const Scalar& scalar, const Cwise& mat)
+    { return mat + scalar; }
+
+    ExpressionType& operator+=(const Scalar& scalar);
+
+    const ScalarAddReturnType
+    operator-(const Scalar& scalar) const;
+
+    ExpressionType& operator-=(const Scalar& scalar);
+
+    template<typename OtherDerived> const typename BinOp<std::less, OtherDerived>::ReturnType
+    operator<(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const typename BinOp<std::less_equal, OtherDerived>::ReturnType
+    operator<=(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const typename BinOp<std::greater, OtherDerived>::ReturnType
+    operator>(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const typename BinOp<std::greater_equal, OtherDerived>::ReturnType
+    operator>=(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const typename BinOp<std::equal_to, OtherDerived>::ReturnType
+    operator==(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const typename BinOp<std::not_equal_to, OtherDerived>::ReturnType
+    operator!=(const MatrixBase<OtherDerived>& other) const;
+
+
+  protected:
+    ExpressionTypeNested m_matrix;
+};
+
+/** \array_module
+  *
+  * \returns a Cwise expression of *this providing additional coefficient-wise operations
+  *
+  * \sa class Cwise
+  */
+template<typename Derived>
+inline const Cwise<Derived>
+MatrixBase<Derived>::cwise() const
+{
+  return derived();
+}
+
+/** \array_module
+  *
+  * \returns a Cwise expression of *this providing additional coefficient-wise operations
+  *
+  * \sa class Cwise
+  */
+template<typename Derived>
+inline Cwise<Derived>
+MatrixBase<Derived>::cwise()
+{
+  return derived();
+}
+
+#endif // EIGEN_CWISE_H
diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h
index 2441b30f6..063ea7853 100644
--- a/Eigen/src/Core/CwiseBinaryOp.h
+++ b/Eigen/src/Core/CwiseBinaryOp.h
@@ -41,7 +41,7 @@
   * However, if you want to write a function returning such an expression, you
   * will need to use this class.
   *
-  * \sa MatrixBase::cwise(const MatrixBase<OtherDerived> &,const CustomBinaryOp &) const, class CwiseUnaryOp, class CwiseNullaryOp
+  * \sa MatrixBase::binaryExpr(const MatrixBase<OtherDerived> &,const CustomBinaryOp &) const, class CwiseUnaryOp, class CwiseNullaryOp
   */
 template<typename BinaryOp, typename Lhs, typename Rhs>
 struct ei_traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
@@ -179,48 +179,48 @@ MatrixBase<Derived>::operator+=(const MatrixBase<OtherDerived>& other)
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<ei_scalar_product_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseProduct(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<ei_scalar_product_op, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator*(const MatrixBase<OtherDerived> &other) const
 {
-  return CwiseBinaryOp<ei_scalar_product_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
+  return typename BinOp<ei_scalar_product_op, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \returns an expression of the coefficient-wise quotient of *this and \a other
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<ei_scalar_quotient_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseQuotient(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<ei_scalar_quotient_op, OtherDerived>::ReturnType
+Cwise<ExpressionType>::operator/(const MatrixBase<OtherDerived> &other) const
 {
-  return CwiseBinaryOp<ei_scalar_quotient_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
+  return typename BinOp<ei_scalar_quotient_op, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \returns an expression of the coefficient-wise min of *this and \a other
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<ei_scalar_min_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseMin(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<ei_scalar_min_op, OtherDerived>::ReturnType
+Cwise<ExpressionType>::min(const MatrixBase<OtherDerived> &other) const
 {
-  return CwiseBinaryOp<ei_scalar_min_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
+  return typename BinOp<ei_scalar_min_op, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \returns an expression of the coefficient-wise max of *this and \a other
   *
   * \sa class CwiseBinaryOp
   */
-template<typename Derived>
+template<typename ExpressionType>
 template<typename OtherDerived>
-inline const CwiseBinaryOp<ei_scalar_max_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-MatrixBase<Derived>::cwiseMax(const MatrixBase<OtherDerived> &other) const
+inline const typename Cwise<ExpressionType>::template BinOp<ei_scalar_max_op, OtherDerived>::ReturnType
+Cwise<ExpressionType>::max(const MatrixBase<OtherDerived> &other) const
 {
-  return CwiseBinaryOp<ei_scalar_max_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
+  return typename BinOp<ei_scalar_max_op, OtherDerived>::ReturnType(_expression(), other.derived());
 }
 
 /** \returns an expression of a custom coefficient-wise operator \a func of *this and \a other
@@ -237,7 +237,7 @@ MatrixBase<Derived>::cwiseMax(const MatrixBase<OtherDerived> &other) const
 template<typename Derived>
 template<typename CustomBinaryOp, typename OtherDerived>
 inline const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
-MatrixBase<Derived>::cwise(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
+MatrixBase<Derived>::binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
 {
   return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func);
 }
diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h
index 78ecb878b..9b15d4240 100644
--- a/Eigen/src/Core/CwiseNullaryOp.h
+++ b/Eigen/src/Core/CwiseNullaryOp.h
@@ -38,7 +38,7 @@
   * However, if you want to write a function returning such an expression, you
   * will need to use this class.
   *
-  * \sa class CwiseUnaryOp, class CwiseBinaryOp, MatrixBase::create()
+  * \sa class CwiseUnaryOp, class CwiseBinaryOp, MatrixBase::NullaryExpr()
   */
 template<typename NullaryOp, typename MatrixType>
 struct ei_traits<CwiseNullaryOp<NullaryOp, MatrixType> >
@@ -123,7 +123,7 @@ class CwiseNullaryOp : ei_no_assignment_operator,
 template<typename Derived>
 template<typename CustomNullaryOp>
 const CwiseNullaryOp<CustomNullaryOp, Derived>
-MatrixBase<Derived>::create(int rows, int cols, const CustomNullaryOp& func)
+MatrixBase<Derived>::NullaryExpr(int rows, int cols, const CustomNullaryOp& func)
 {
   return CwiseNullaryOp<CustomNullaryOp, Derived>(rows, cols, func);
 }
@@ -146,7 +146,7 @@ MatrixBase<Derived>::create(int rows, int cols, const CustomNullaryOp& func)
 template<typename Derived>
 template<typename CustomNullaryOp>
 const CwiseNullaryOp<CustomNullaryOp, Derived>
-MatrixBase<Derived>::create(int size, const CustomNullaryOp& func)
+MatrixBase<Derived>::NullaryExpr(int size, const CustomNullaryOp& func)
 {
   ei_assert(IsVectorAtCompileTime);
   if(RowsAtCompileTime == 1) return CwiseNullaryOp<CustomNullaryOp, Derived>(1, size, func);
@@ -165,7 +165,7 @@ MatrixBase<Derived>::create(int size, const CustomNullaryOp& func)
 template<typename Derived>
 template<typename CustomNullaryOp>
 const CwiseNullaryOp<CustomNullaryOp, Derived>
-MatrixBase<Derived>::create(const CustomNullaryOp& func)
+MatrixBase<Derived>::NullaryExpr(const CustomNullaryOp& func)
 {
   return CwiseNullaryOp<CustomNullaryOp, Derived>(RowsAtCompileTime, ColsAtCompileTime, func);
 }
@@ -187,7 +187,7 @@ template<typename Derived>
 const typename MatrixBase<Derived>::ConstantReturnType
 MatrixBase<Derived>::constant(int rows, int cols, const Scalar& value)
 {
-  return create(rows, cols, ei_scalar_constant_op<Scalar>(value));
+  return NullaryExpr(rows, cols, ei_scalar_constant_op<Scalar>(value));
 }
 
 /** \returns an expression of a constant matrix of value \a value
@@ -209,7 +209,7 @@ template<typename Derived>
 const typename MatrixBase<Derived>::ConstantReturnType
 MatrixBase<Derived>::constant(int size, const Scalar& value)
 {
-  return create(size, ei_scalar_constant_op<Scalar>(value));
+  return NullaryExpr(size, ei_scalar_constant_op<Scalar>(value));
 }
 
 /** \returns an expression of a constant matrix of value \a value
@@ -225,7 +225,7 @@ template<typename Derived>
 const typename MatrixBase<Derived>::ConstantReturnType
 MatrixBase<Derived>::constant(const Scalar& value)
 {
-  return create(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_constant_op<Scalar>(value));
+  return NullaryExpr(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_constant_op<Scalar>(value));
 }
 
 template<typename Derived>
@@ -455,7 +455,7 @@ template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_identity_op<typename ei_traits<Derived>::Scalar>, Derived>
 MatrixBase<Derived>::identity(int rows, int cols)
 {
-  return create(rows, cols, ei_scalar_identity_op<Scalar>());
+  return NullaryExpr(rows, cols, ei_scalar_identity_op<Scalar>());
 }
 
 /** \returns an expression of the identity matrix (not necessarily square).
@@ -472,7 +472,7 @@ template<typename Derived>
 inline const CwiseNullaryOp<ei_scalar_identity_op<typename ei_traits<Derived>::Scalar>, Derived>
 MatrixBase<Derived>::identity()
 {
-  return create(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_identity_op<Scalar>());
+  return NullaryExpr(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_identity_op<Scalar>());
 }
 
 /** \returns true if *this is approximately equal to the identity matrix
diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h
index c80124b47..68ef9b4d4 100644
--- a/Eigen/src/Core/CwiseUnaryOp.h
+++ b/Eigen/src/Core/CwiseUnaryOp.h
@@ -37,7 +37,7 @@
   * It is the return type of the unary operator-, of a matrix or a vector, and most
   * of the time this is the only way it is used.
   *
-  * \sa MatrixBase::cwise(const CustomUnaryOp &) const, class CwiseBinaryOp, class CwiseNullaryOp
+  * \sa MatrixBase::unaryExpr(const CustomUnaryOp &) const, class CwiseBinaryOp, class CwiseNullaryOp
   */
 template<typename UnaryOp, typename MatrixType>
 struct ei_traits<CwiseUnaryOp<UnaryOp, MatrixType> >
@@ -118,7 +118,7 @@ class CwiseUnaryOp : ei_no_assignment_operator,
 template<typename Derived>
 template<typename CustomUnaryOp>
 inline const CwiseUnaryOp<CustomUnaryOp, Derived>
-MatrixBase<Derived>::cwise(const CustomUnaryOp& func) const
+MatrixBase<Derived>::unaryExpr(const CustomUnaryOp& func) const
 {
   return CwiseUnaryOp<CustomUnaryOp, Derived>(derived(), func);
 }
@@ -134,20 +134,20 @@ MatrixBase<Derived>::operator-() const
 
 /** \returns an expression of the coefficient-wise absolute value of \c *this
   */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_abs_op<typename ei_traits<Derived>::Scalar>,Derived>
-MatrixBase<Derived>::cwiseAbs() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_abs_op>::ReturnType
+Cwise<ExpressionType>::abs() const
 {
-  return derived();
+  return _expression();
 }
 
 /** \returns an expression of the coefficient-wise squared absolute value of \c *this
   */
-template<typename Derived>
-inline const CwiseUnaryOp<ei_scalar_abs2_op<typename ei_traits<Derived>::Scalar>,Derived>
-MatrixBase<Derived>::cwiseAbs2() const
+template<typename ExpressionType>
+inline const typename Cwise<ExpressionType>::template UnOp<ei_scalar_abs2_op>::ReturnType
+Cwise<ExpressionType>::abs2() const
 {
-  return derived();
+  return _expression();
 }
 
 /** \returns an expression of the complex conjugate of \c *this.
diff --git a/Eigen/src/Core/Fuzzy.h b/Eigen/src/Core/Fuzzy.h
index 0b2305c57..b48c198b3 100644
--- a/Eigen/src/Core/Fuzzy.h
+++ b/Eigen/src/Core/Fuzzy.h
@@ -54,7 +54,7 @@ bool MatrixBase<Derived>::isApprox(
 {
   const typename ei_nested<Derived,2>::type nested(derived());
   const typename ei_nested<OtherDerived,2>::type otherNested(other.derived());
-  return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * std::min(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
+  return (nested - otherNested).cwise().abs2().sum() <= prec * prec * std::min(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum());
 }
 
 /** \returns \c true if the norm of \c *this is much smaller than \a other,
@@ -76,7 +76,7 @@ bool MatrixBase<Derived>::isMuchSmallerThan(
   typename NumTraits<Scalar>::Real prec
 ) const
 {
-  return cwiseAbs2().sum() <= prec * prec * other * other;
+  return cwise().abs2().sum() <= prec * prec * other * other;
 }
 
 /** \returns \c true if the norm of \c *this is much smaller than the norm of \a other,
@@ -96,7 +96,7 @@ bool MatrixBase<Derived>::isMuchSmallerThan(
   typename NumTraits<Scalar>::Real prec
 ) const
 {
-  return this->cwiseAbs2().sum() <= prec * prec * other.cwiseAbs2().sum();
+  return this->cwise().abs2().sum() <= prec * prec * other.cwise().abs2().sum();
 }
 
 #else
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 63817102f..95c188621 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -371,13 +371,13 @@ template<typename Derived> class MatrixBase
 
     template<typename CustomNullaryOp>
     static const CwiseNullaryOp<CustomNullaryOp, Derived>
-    create(int rows, int cols, const CustomNullaryOp& func);
+    NullaryExpr(int rows, int cols, const CustomNullaryOp& func);
     template<typename CustomNullaryOp>
     static const CwiseNullaryOp<CustomNullaryOp, Derived>
-    create(int size, const CustomNullaryOp& func);
+    NullaryExpr(int size, const CustomNullaryOp& func);
     template<typename CustomNullaryOp>
     static const CwiseNullaryOp<CustomNullaryOp, Derived>
-    create(const CustomNullaryOp& func);
+    NullaryExpr(const CustomNullaryOp& func);
 
     static const ConstantReturnType zero(int rows, int cols);
     static const ConstantReturnType zero(int size);
@@ -457,31 +457,12 @@ template<typename Derived> class MatrixBase
     const ConjugateReturnType conjugate() const;
     const RealReturnType real() const;
 
-    template<typename OtherDerived>
-    const CwiseBinaryOp<ei_scalar_product_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseProduct(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<ei_scalar_quotient_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseQuotient(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<ei_scalar_min_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseMin(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<ei_scalar_max_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseMax(const MatrixBase<OtherDerived> &other) const;
-
-    const CwiseUnaryOp<ei_scalar_abs_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseAbs() const;
-    const CwiseUnaryOp<ei_scalar_abs2_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseAbs2() const;
-
     template<typename CustomUnaryOp>
-    const CwiseUnaryOp<CustomUnaryOp, Derived> cwise(const CustomUnaryOp& func = CustomUnaryOp()) const;
+    const CwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const;
 
     template<typename CustomBinaryOp, typename OtherDerived>
     const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
-    cwise(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const;
+    binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const;
 
 
     Scalar sum() const;
@@ -506,45 +487,11 @@ template<typename Derived> class MatrixBase
     inline Derived& const_cast_derived() const
     { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
 
+    const Cwise<Derived> cwise() const;
+    Cwise<Derived> cwise();
 
 /////////// Array module ///////////
 
-    const Array<Derived> array() const;
-    Array<Derived> array();
-
-    const CwiseUnaryOp<ei_scalar_sqrt_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseSqrt() const;
-    const CwiseUnaryOp<ei_scalar_exp_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseExp() const;
-    const CwiseUnaryOp<ei_scalar_log_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseLog() const;
-    const CwiseUnaryOp<ei_scalar_cos_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseCos() const;
-    const CwiseUnaryOp<ei_scalar_sin_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseSin() const;
-    const CwiseUnaryOp<ei_scalar_pow_op<typename ei_traits<Derived>::Scalar>, Derived>
-    cwisePow(const Scalar& exponent) const;
-    const CwiseUnaryOp<ei_scalar_inverse_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseInverse() const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::less<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseLessThan(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::less_equal<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseLessEqual(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::greater<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseGreaterThan(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::greater_equal<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseGreaterEqual(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::equal_to<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseEqualTo(const MatrixBase<OtherDerived> &other) const;
-
-    template<typename OtherDerived>
-    const CwiseBinaryOp<std::not_equal_to<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
-    cwiseNotEqualTo(const MatrixBase<OtherDerived> &other) const;
-
     bool all(void) const;
     bool any(void) const;
 
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index bad5cacdc..922a3716f 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -55,7 +55,7 @@ template<typename MatrixType, int Alignment = Unaligned> class Map;
 template<int Direction, typename UnaryOp, typename MatrixType> class PartialRedux;
 template<typename MatrixType, unsigned int Mode> class Part;
 template<typename MatrixType, unsigned int Mode> class Extract;
-template<typename MatrixType> class Array;
+template<typename ExpressionType> class Cwise;
 
 template<typename Lhs, typename Rhs> struct ei_product_mode;
 template<typename Lhs, typename Rhs, int ProductMode = ei_product_mode<Lhs,Rhs>::value> struct ProductReturnType;
@@ -76,6 +76,8 @@ template<typename Scalar> struct ei_scalar_cos_op;
 template<typename Scalar> struct ei_scalar_sin_op;
 template<typename Scalar> struct ei_scalar_pow_op;
 template<typename Scalar> struct ei_scalar_inverse_op;
+template<typename Scalar> struct ei_scalar_square_op;
+template<typename Scalar> struct ei_scalar_cube_op;
 template<typename Scalar, typename NewType> struct ei_scalar_cast_op;
 template<typename Scalar, bool PacketAccess> struct ei_scalar_multiple_op;
 template<typename Scalar> struct ei_scalar_quotient1_op;
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index b7796bb29..e7bfc195d 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -61,10 +61,10 @@ template<> class ei_int_if_dynamic<Dynamic>
 };
 
 
-template <bool Condition, class Then, class Else>
+template<bool Condition, typename Then, typename Else>
 struct ei_meta_if { typedef Then ret; };
 
-template <class Then, class Else>
+template<typename Then, typename Else>
 struct ei_meta_if <false, Then, Else> { typedef Else ret; };
 
 template<typename T, typename U> struct ei_is_same_type { enum { ret = 0 }; };
diff --git a/Eigen/src/Geometry/AngleAxis.h b/Eigen/src/Geometry/AngleAxis.h
index 2b0ee6196..c77264e1c 100644
--- a/Eigen/src/Geometry/AngleAxis.h
+++ b/Eigen/src/Geometry/AngleAxis.h
@@ -139,7 +139,7 @@ AngleAxis<Scalar>::toRotationMatrix(void) const
   res.coeffRef(1,2) = tmp - sin_axis.x();
   res.coeffRef(2,1) = tmp + sin_axis.x();
 
-  res.diagonal() = Vector3::constant(c) + cos1_axis.cwiseProduct(m_axis);
+  res.diagonal() = (cos1_axis.cwise() * m_axis).cwise() + c;
 
   return res;
 }
diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h
index f08f82151..d8ccf4917 100644
--- a/Eigen/src/Geometry/Quaternion.h
+++ b/Eigen/src/Geometry/Quaternion.h
@@ -212,8 +212,8 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(EulerAnglesType ea)
 {
   ea.coeffs() *= 0.5;
 
-  Vector3 cosines = ea.coeffs().cwiseCos();
-  Vector3 sines   = ea.coeffs().cwiseSin();
+  Vector3 cosines = ea.coeffs().cwise().cos();
+  Vector3 sines   = ea.coeffs().cwise().sin();
 
   Scalar cYcZ = cosines.y() * cosines.z();
   Scalar sYsZ = sines.y() * sines.z();
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index a4f703383..df8a22ebe 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -107,14 +107,14 @@ void Inverse<MatrixType, CheckExistence>
 ::_compute_in_general_case(const MatrixType& _matrix)
 {
   MatrixType matrix(_matrix);
-  const RealScalar max = CheckExistence ? matrix.cwiseAbs().maxCoeff()
+  const RealScalar max = CheckExistence ? matrix.cwise().abs().maxCoeff()
                                         : static_cast<RealScalar>(0);
   const int size = matrix.rows();
   for(int k = 0; k < size-1; k++)
   {
     int rowOfBiggest;
     const RealScalar max_in_this_col
-      = matrix.col(k).end(size-k).cwiseAbs().maxCoeff(&rowOfBiggest);
+      = matrix.col(k).end(size-k).cwise().abs().maxCoeff(&rowOfBiggest);
     if(CheckExistence && ei_isMuchSmallerThan(max_in_this_col, max))
     { m_exists = false; return; }
 
@@ -150,7 +150,7 @@ bool ei_compute_size2_inverse(const ExpressionType& xpr, typename ExpressionType
   typedef typename ExpressionType::Scalar Scalar;
   const typename ei_nested<ExpressionType, 1+CheckExistence>::type matrix(xpr);
   const Scalar det = matrix.determinant();
-  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff()))
+  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff()))
     return false;
   const Scalar invdet = static_cast<Scalar>(1) / det;
   result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
@@ -169,7 +169,7 @@ void Inverse<MatrixType, CheckExistence>::_compute_in_size3_case(const MatrixTyp
   const Scalar det = det_minor00 * matrix.coeff(0,0)
                    - det_minor10 * matrix.coeff(1,0)
                    + det_minor20 * matrix.coeff(2,0);
-  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwiseAbs().maxCoeff()))
+  if(CheckExistence && ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff()))
     m_exists = false;
   else
   {
diff --git a/bench/vdw_new.cpp b/bench/vdw_new.cpp
index b1237c4c7..a1bd422ae 100644
--- a/bench/vdw_new.cpp
+++ b/bench/vdw_new.cpp
@@ -1,35 +1,5 @@
-#include <Eigen/Core>
+#include <Eigen/Array>
 
-namespace Eigen {
-
-template<typename Scalar> struct pow12_op EIGEN_EMPTY_STRUCT {
-  inline const Scalar operator() (const Scalar& a) const
-  {
-      Scalar b = a*a*a;
-      Scalar c = b*b;
-      return c*c;
-  }
-  template<typename PacketScalar>
-  inline const PacketScalar packetOp(const PacketScalar& a) const
-  {
-      PacketScalar b = ei_pmul(a, ei_pmul(a, a));
-      PacketScalar c = ei_pmul(b, b);
-      return ei_pmul(c, c);
-  }
-};
-
-template<typename Scalar>
-struct ei_functor_traits<pow12_op<Scalar> >
-{
-  enum {
-    Cost = 4 * NumTraits<Scalar>::MulCost,
-    PacketAccess = int(ei_packet_traits<Scalar>::size) > 1
-  };
-};
-
-} // namespace Eigen
-
-using Eigen::pow12_op;
 USING_PART_OF_NAMESPACE_EIGEN
 
 #ifndef SCALAR
@@ -50,9 +20,10 @@ using namespace std;
 
 SCALAR E_VDW(const Vec &interactions1, const Vec &interactions2)
 {
-  return interactions2
-         .cwiseQuotient(interactions1)
-         .cwise(pow12_op<SCALAR>())
+  return (interactions2.cwise()/interactions1)
+         .cwise().cube()
+         .cwise().square()
+         .cwise().square()
          .sum();
 }
 
diff --git a/demos/mandelbrot/mandelbrot.cpp b/demos/mandelbrot/mandelbrot.cpp
index 42f2a484e..2aaeab1e5 100644
--- a/demos/mandelbrot/mandelbrot.cpp
+++ b/demos/mandelbrot/mandelbrot.cpp
@@ -1,3 +1,27 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.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/>.
+
 #include "mandelbrot.h"
 #include<QtGui/QPainter>
 #include<QtGui/QImage>
diff --git a/demos/mandelbrot/mandelbrot.h b/demos/mandelbrot/mandelbrot.h
index c8fad5505..63eb81a68 100644
--- a/demos/mandelbrot/mandelbrot.h
+++ b/demos/mandelbrot/mandelbrot.h
@@ -1,3 +1,27 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.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 MANDELBROT_H
 #define MANDELBROT_H
 
diff --git a/test/array.cpp b/test/array.cpp
index d05c4f72e..6233aaec8 100644
--- a/test/array.cpp
+++ b/test/array.cpp
@@ -44,15 +44,15 @@ template<typename MatrixType> void scalarAdd(const MatrixType& m)
   Scalar  s1 = ei_random<Scalar>(),
           s2 = ei_random<Scalar>();
 
-  VERIFY_IS_APPROX(m1.array() + s1, s1 + m1.array());
-  VERIFY_IS_APPROX(m1.array() + s1, MatrixType::constant(rows,cols,s1) + m1);
-  VERIFY_IS_APPROX((m1*Scalar(2)).array() - s2, (m1+m1) - MatrixType::constant(rows,cols,s2) );
+  VERIFY_IS_APPROX(m1.cwise() + s1, s1 + m1.cwise());
+  VERIFY_IS_APPROX(m1.cwise() + s1, MatrixType::constant(rows,cols,s1) + m1);
+  VERIFY_IS_APPROX((m1*Scalar(2)).cwise() - s2, (m1+m1) - MatrixType::constant(rows,cols,s2) );
   m3 = m1;
-  m3.array() += s2;
-  VERIFY_IS_APPROX(m3, m1.array() + s2);
+  m3.cwise() += s2;
+  VERIFY_IS_APPROX(m3, m1.cwise() + s2);
   m3 = m1;
-  m3.array() -= s1;
-  VERIFY_IS_APPROX(m3, m1.array() - s1);
+  m3.cwise() -= s1;
+  VERIFY_IS_APPROX(m3, m1.cwise() - s1);
 }
 
 template<typename MatrixType> void comparisons(const MatrixType& m)
@@ -70,14 +70,14 @@ template<typename MatrixType> void comparisons(const MatrixType& m)
              m2 = MatrixType::random(rows, cols),
              m3(rows, cols);
 
-  VERIFY((m1.array() + Scalar(1)).array() > m1.array());
-  VERIFY((m1.array() - Scalar(1)).array() < m1.array());
+  VERIFY(((m1.cwise() + Scalar(1)).cwise() > m1).all());
+  VERIFY(((m1.cwise() - Scalar(1)).cwise() < m1).all());
   if (rows*cols>1)
   {
     m3 = m1;
     m3(r,c) += 1;
-    VERIFY(! (m1.array() < m3.array()) );
-    VERIFY(! (m1.array() > m3.array()) );
+    VERIFY(! (m1.cwise() < m3).all() );
+    VERIFY(! (m1.cwise() > m3).all() );
   }
 }
 
diff --git a/test/cwiseop.cpp b/test/cwiseop.cpp
index 8f7ce9953..9278dd08d 100644
--- a/test/cwiseop.cpp
+++ b/test/cwiseop.cpp
@@ -55,7 +55,7 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
              v2 = VectorType::random(rows),
              vzero = VectorType::zero(rows);
 
-  m2 = m2.template cwise<AddIfNull<Scalar> >(mones);
+  m2 = m2.template binaryExpr<AddIfNull<Scalar> >(mones);
 
   VERIFY_IS_APPROX(               mzero,    m1-m1);
   VERIFY_IS_APPROX(               m2,       m1+m2-m1);
@@ -63,13 +63,13 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
   if(NumTraits<Scalar>::HasFloatingPoint)
 #endif
   {
-    VERIFY_IS_APPROX(               mones,    m2.cwiseQuotient(m2));
+    VERIFY_IS_APPROX(               mones,    m2.cwise()/m2);
   }
-  VERIFY_IS_APPROX(               m1.cwiseProduct(m2),    m2.cwiseProduct(m1));
+  VERIFY_IS_APPROX(               m1.cwise() * m2,    m2.cwise() * m1);
 
-  VERIFY( m1.cwiseLessThan(m1.cwise(bind2nd(plus<Scalar>(), Scalar(1)))).all() );
-  VERIFY( !m1.cwiseLessThan(m1.cwise(bind2nd(minus<Scalar>(), Scalar(1)))).all() );
-  VERIFY( !m1.cwiseGreaterThan(m1.cwise(bind2nd(plus<Scalar>(), Scalar(1)))).any() );
+  VERIFY( (m1.cwise()<m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).all() );
+  VERIFY( !(m1.cwise()<m1.unaryExpr(bind2nd(minus<Scalar>(), Scalar(1)))).all() );
+  VERIFY( !(m1.cwise()>m1.unaryExpr(bind2nd(plus<Scalar>(), Scalar(1)))).any() );
   //VERIFY_IS_APPROX(               m1,       m2.cwiseProduct(m1).cwiseQuotient(m2));
 
 //   VERIFY_IS_APPROX( cwiseMin(m1,m2), cwiseMin(m2,m1) );
diff --git a/test/geometry.cpp b/test/geometry.cpp
index fd20c137f..64b096aa8 100644
--- a/test/geometry.cpp
+++ b/test/geometry.cpp
@@ -65,7 +65,7 @@ template<typename Scalar> void geometry(void)
   VERIFY_IS_APPROX(Quaternion(EulerAngles(q1)) * v1, q1 * v1);
   EulerAngles ea = q2;
   VERIFY_IS_APPROX(EulerAngles(Quaternion(ea)).coeffs(), ea.coeffs());
-  VERIFY_IS_NOT_APPROX(EulerAngles(Quaternion(EulerAngles(v2.cwiseProduct(Vector3(0.2,-0.2,1))))).coeffs(), v2);
+  VERIFY_IS_NOT_APPROX(EulerAngles(Quaternion(EulerAngles(v2.cwise() * Vector3(0.2,-0.2,1)))).coeffs(), v2);
 
   // angle-axis conversion
   AngleAxis aa = q1;
@@ -128,7 +128,7 @@ template<typename Scalar> void geometry(void)
   t0.pretranslate(v0);
   t0.scale(v1);
   t1.affine() = q1.conjugate().toRotationMatrix();
-  t1.prescale(v1.cwiseInverse());
+  t1.prescale(v1.cwise().inverse());
   t1.translate(-v0);
 
   VERIFY((t0.matrix() * t1.matrix()).isIdentity());
@@ -147,7 +147,7 @@ template<typename Scalar> void geometry(void)
 
   t21.setIdentity();
   t21.affine() = Rotation2D<Scalar>(-a).toRotationMatrix();
-  VERIFY( (t20.fromPositionOrientationScale(v20,a,v21) * (t21.prescale(v21.cwiseInverse()).translate(-v20))).isIdentity() );
+  VERIFY( (t20.fromPositionOrientationScale(v20,a,v21) * (t21.prescale(v21.cwise().inverse()).translate(-v20))).isIdentity() );
 }
 
 void test_geometry()
diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp
index 4c2bdf620..8096d2f3e 100644
--- a/test/linearstructure.cpp
+++ b/test/linearstructure.cpp
@@ -85,7 +85,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
 
   // use .block to disable vectorization and compare to the vectorized version
   VERIFY_IS_APPROX(m1+m1.block(0,0,rows,cols), m1+m1);
-  VERIFY_IS_APPROX(m1.cwiseProduct(m1.block(0,0,rows,cols)), m1.cwiseProduct(m1));
+  VERIFY_IS_APPROX(m1.cwise() * m1.block(0,0,rows,cols), m1.cwise() * m1);
   VERIFY_IS_APPROX(m1 - m1.block(0,0,rows,cols), m1 - m1);
   VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1);
 }
diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp
index d6c566bef..d9653588d 100644
--- a/test/nomalloc.cpp
+++ b/test/nomalloc.cpp
@@ -75,7 +75,7 @@ template<typename MatrixType> void nomalloc(const MatrixType& m)
 
   VERIFY_IS_APPROX((m1+m2)*s1,              s1*m1+s1*m2);
   VERIFY_IS_APPROX((m1+m2)(r,c), (m1(r,c))+(m2(r,c)));
-  VERIFY_IS_APPROX(m1.cwiseProduct(m1.block(0,0,rows,cols)), m1.cwiseProduct(m1));
+  VERIFY_IS_APPROX(m1.cwise() * m1.block(0,0,rows,cols), m1.cwise() * m1);
   VERIFY_IS_APPROX((m1*m1.transpose())*m2,  m1*(m1.transpose()*m2));
 }