17 files changed, 257 insertions, 338 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 7fc048729..eda3fd82a 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -16,11 +16,9 @@ namespace Eigen {
 #include "src/Core/MatrixRef.h"
 #include "src/Core/MatrixStorage.h"
 #include "src/Core/Matrix.h"
-#include "src/Core/Cast.h"
 #include "src/Core/Eval.h"
 #include "src/Core/CwiseBinaryOp.h"
 #include "src/Core/CwiseUnaryOp.h"
-#include "src/Core/ScalarMultiple.h"
 #include "src/Core/Product.h"
 #include "src/Core/Row.h"
 #include "src/Core/Column.h"
diff --git a/Eigen/src/Core/Cast.h b/Eigen/src/Core/Cast.h
deleted file mode 100644
index 29d24967f..000000000
--- a/Eigen/src/Core/Cast.h
+++ /dev/null
@@ -1,98 +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) 2006-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_CAST_H
-#define EIGEN_CAST_H
-
-/** \class Cast
-  *
-  * \brief Expression with casted scalar type
-  *
-  * \param NewScalar the new scalar type
-  * \param MatrixType the type of the object in which we are casting the scalar type
-  *
-  * This class represents an expression where we are casting the scalar type to a new
-  * type. It is the return type of MatrixBase::cast() and most of the time this is the
-  * only way it is used.
-  *
-  * However, if you want to write a function returning such an expression, you
-  * will need to use this class.
-  *
-  * Here is an example illustrating this:
-  * \include class_Cast.cpp
-  * Output: \verbinclude class_Cast.out
-  *
-  * \sa MatrixBase::cast()
-  */
-template<typename NewScalar, typename MatrixType> class Cast : NoOperatorEquals,
-  public MatrixBase<NewScalar, Cast<NewScalar, MatrixType> >
-{
-  public:
-    typedef NewScalar Scalar;
-    typedef typename MatrixType::AsArg MatRef;
-    friend class MatrixBase<Scalar, Cast>;
-    friend class MatrixBase<Scalar, Cast>::Traits;
-    typedef MatrixBase<Scalar, Cast> Base;
-
-    Cast(const MatRef& matrix) : m_matrix(matrix) {}
-
-  private:
-    enum {
-      RowsAtCompileTime = MatrixType::Traits::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::Traits::ColsAtCompileTime,
-      MaxRowsAtCompileTime = MatrixType::Traits::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::Traits::MaxColsAtCompileTime
-    };
-    const Cast& _asArg() const { return *this; }
-    int _rows() const { return m_matrix.rows(); }
-    int _cols() const { return m_matrix.cols(); }
-
-    NewScalar _coeff(int row, int col) const
-    {
-      return static_cast<NewScalar>(m_matrix.coeff(row, col));
-    }
-
-  protected:
-    const MatRef m_matrix;
-};
-
-/** \returns an expression of *this with the \a Scalar type casted to
-  * \a NewScalar.
-  *
-  * The template parameter \a NewScalar is the type we are casting the scalars to.
-  *
-  * Example: \include MatrixBase_cast.cpp
-  * Output: \verbinclude MatrixBase_cast.out
-  *
-  * \sa class Cast
-  */
-template<typename Scalar, typename Derived>
-template<typename NewScalar>
-const Cast<NewScalar, Derived>
-MatrixBase<Scalar, Derived>::cast() const
-{
-  return Cast<NewScalar, Derived>(asArg());
-}
-
-#endif // EIGEN_CAST_H
diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h
index a0e9b3cd3..cdcc329b5 100644
--- a/Eigen/src/Core/CwiseBinaryOp.h
+++ b/Eigen/src/Core/CwiseBinaryOp.h
@@ -48,18 +48,20 @@
   */
 template<typename BinaryOp, typename Lhs, typename Rhs>
 class CwiseBinaryOp : NoOperatorEquals,
-  public MatrixBase<typename Lhs::Scalar, CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
+  public MatrixBase<
+            typename ei_result_of<BinaryOp(typename Lhs::Scalar,typename Rhs::Scalar)>::type,
+            CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
 {
   public:
-    typedef typename Lhs::Scalar Scalar;
+    typedef typename ei_result_of<BinaryOp(typename Lhs::Scalar,typename Rhs::Scalar)>::type Scalar;
     typedef typename Lhs::AsArg LhsRef;
     typedef typename Rhs::AsArg RhsRef;
     friend class MatrixBase<Scalar, CwiseBinaryOp>;
     friend class MatrixBase<Scalar, CwiseBinaryOp>::Traits;
     typedef MatrixBase<Scalar, CwiseBinaryOp> Base;
 
-    CwiseBinaryOp(const LhsRef& lhs, const RhsRef& rhs)
-      : m_lhs(lhs), m_rhs(rhs)
+    CwiseBinaryOp(const LhsRef& lhs, const RhsRef& rhs, const BinaryOp& func = BinaryOp())
+      : m_lhs(lhs), m_rhs(rhs), m_functor(func)
     {
       assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
     }
@@ -78,44 +80,45 @@ class CwiseBinaryOp : NoOperatorEquals,
 
     Scalar _coeff(int row, int col) const
     {
-      return BinaryOp::template op<Scalar>(m_lhs.coeff(row, col), m_rhs.coeff(row, col));
+      return m_functor(m_lhs.coeff(row, col), m_rhs.coeff(row, col));
     }
 
   protected:
     const LhsRef m_lhs;
     const RhsRef m_rhs;
+    const BinaryOp m_functor;
 };
 
 /** \brief Template functor to compute the sum of two scalars
   *
   * \sa class CwiseBinaryOp, MatrixBase::operator+
   */
-struct CwiseSumOp {
-    template<typename Scalar> static Scalar op(const Scalar& a, const Scalar& b) { return a + b; }
+struct CwiseSumOp EIGEN_EMPTY_STRUCT {
+    template<typename Scalar> Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; }
 };
 
 /** \brief Template functor to compute the difference of two scalars
   *
   * \sa class CwiseBinaryOp, MatrixBase::operator-
   */
-struct CwiseDifferenceOp {
-    template<typename Scalar> static Scalar op(const Scalar& a, const Scalar& b) { return a - b; }
+struct CwiseDifferenceOp EIGEN_EMPTY_STRUCT {
+    template<typename Scalar> Scalar operator() (const Scalar& a, const Scalar& b) const { return a - b; }
 };
 
 /** \brief Template functor to compute the product of two scalars
   *
   * \sa class CwiseBinaryOp, MatrixBase::cwiseProduct()
   */
-struct ScalarProductOp {
-    template<typename Scalar> static Scalar op(const Scalar& a, const Scalar& b) { return a * b; }
+struct ScalarProductOp EIGEN_EMPTY_STRUCT {
+    template<typename Scalar> Scalar operator() (const Scalar& a, const Scalar& b) const { return a * b; }
 };
 
 /** \brief Template functor to compute the quotient of two scalars
   *
   * \sa class CwiseBinaryOp, MatrixBase::cwiseQuotient()
   */
-struct ScalarQuotientOp {
-    template<typename Scalar> static Scalar op(const Scalar& a, const Scalar& b) { return a / b; }
+struct ScalarQuotientOp EIGEN_EMPTY_STRUCT {
+    template<typename Scalar> Scalar operator() (const Scalar& a, const Scalar& b) const { return a / b; }
 };
 
 /** \relates MatrixBase
@@ -196,43 +199,19 @@ MatrixBase<Scalar, Derived>::cwiseQuotient(const MatrixBase<Scalar, OtherDerived
 }
 
 
-/** \relates MatrixBase
-  *
-  * \returns an expression of a custom coefficient-wise operator of \a mat1 and \a mat2
-  *
-  * The template parameter \a CustomBinaryOp is the template functor
-  * of the custom operator (see class CwiseBinaryOp for an example)
-  *
-  * \sa class CwiseBinaryOp, MatrixBase::operator+, MatrixBase::operator-, MatrixBase::cwiseProduct, MatrixBase::cwiseQuotient
-  */
-template<typename CustomBinaryOp, typename Scalar, typename Derived1, typename Derived2>
-const CwiseBinaryOp<CustomBinaryOp, Derived1, Derived2>
-cwise(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
-{
-  return CwiseBinaryOp<CustomBinaryOp, Derived1, Derived2>(mat1.asArg(), mat2.asArg());
-}
-
-/** \returns an expression of a custom coefficient-wise operator of *this and \a other
+/** \returns an expression of a custom coefficient-wise operator \a func of *this and \a other
   *
-  * The template parameter \a CustomBinaryOp is the template functor
+  * The template parameter \a CustomBinaryOp is the type of the functor
   * of the custom operator (see class CwiseBinaryOp for an example)
-  *
-  * \note since cwise is a templated member with a mandatory template parameter,
-  * the keyword template as to be used if the matrix type is also a template parameter:
-  * \code
-  * template <typename MatrixType> void foo(const MatrixType& m1, const MatrixType& m2) {
-  *   m1.template cwise<ScalarProductOp>(m2);
-  * }
-  * \endcode
-  *
+  * 
   * \sa class CwiseBinaryOp, MatrixBase::operator+, MatrixBase::operator-, MatrixBase::cwiseProduct, MatrixBase::cwiseQuotient
   */
 template<typename Scalar, typename Derived>
 template<typename CustomBinaryOp, typename OtherDerived>
 const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
-MatrixBase<Scalar, Derived>::cwise(const MatrixBase<Scalar, OtherDerived> &other) const
+MatrixBase<Scalar, Derived>::cwise(const MatrixBase<Scalar, OtherDerived> &other, const CustomBinaryOp& func) const
 {
-  return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(asArg(), other.asArg());
+  return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(asArg(), other.asArg(), func);
 }
 
 
diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h
index 7dc912416..c51ad26bc 100644
--- a/Eigen/src/Core/CwiseUnaryOp.h
+++ b/Eigen/src/Core/CwiseUnaryOp.h
@@ -41,16 +41,18 @@
   */
 template<typename UnaryOp, typename MatrixType>
 class CwiseUnaryOp : NoOperatorEquals,
-  public MatrixBase<typename MatrixType::Scalar, CwiseUnaryOp<UnaryOp, MatrixType> >
+  public MatrixBase<
+    typename ei_result_of<UnaryOp(typename MatrixType::Scalar)>::type,
+    CwiseUnaryOp<UnaryOp, MatrixType> >
 {
   public:
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename ei_result_of<UnaryOp(typename MatrixType::Scalar)>::type Scalar;
     typedef typename MatrixType::AsArg MatRef;
     friend class MatrixBase<Scalar, CwiseUnaryOp>;
     friend class MatrixBase<Scalar, CwiseUnaryOp>::Traits;
     typedef MatrixBase<Scalar, CwiseUnaryOp> Base;
 
-    CwiseUnaryOp(const MatRef& mat) : m_matrix(mat) {}
+    CwiseUnaryOp(const MatRef& mat, const UnaryOp& func = UnaryOp()) : m_matrix(mat), m_functor(func) {}
 
   private:
     enum {
@@ -66,27 +68,28 @@ class CwiseUnaryOp : NoOperatorEquals,
 
     Scalar _coeff(int row, int col) const
     {
-      return UnaryOp::template op<Scalar>(m_matrix.coeff(row, col));
+      return m_functor(m_matrix.coeff(row, col));
     }
 
   protected:
     const MatRef m_matrix;
+    const UnaryOp m_functor;
 };
 
 /** \brief Template functor to compute the opposite of a scalar
   *
   * \sa class CwiseUnaryOp, MatrixBase::operator-
   */
-struct ScalarOppositeOp {
-  template<typename Scalar> static Scalar op(const Scalar& a) { return -a; }
+struct ScalarOppositeOp EIGEN_EMPTY_STRUCT {
+  template<typename Scalar> Scalar operator() (const Scalar& a) const { return -a; }
 };
 
 /** \brief Template functor to compute the absolute value of a scalar
   *
   * \sa class CwiseUnaryOp, MatrixBase::cwiseAbs
   */
-struct ScalarAbsOp {
-  template<typename Scalar> static Scalar op(const Scalar& a) { return ei_abs(a); }
+struct ScalarAbsOp EIGEN_EMPTY_STRUCT {
+  template<typename Scalar> Scalar operator() (const Scalar& a) const { return ei_abs(a); }
 };
 
 
@@ -109,43 +112,22 @@ MatrixBase<Scalar, Derived>::cwiseAbs() const
 }
 
 
-/** \relates MatrixBase
+/** \returns an expression of a custom coefficient-wise unary operator \a func of *this
   *
-  * \returns an expression of a custom coefficient-wise unary operator of \a mat
-  *
-  * The template parameter \a CustomUnaryOp is the template functor
+  * The template parameter \a CustomUnaryOp is the type of the functor
   * of the custom unary operator.
-  *
-  * \sa class CwiseUnaryOp, class CwiseBinarOp, MatrixBase::cwise, MatrixBase::operator-, MatrixBase::cwiseAbs
-  */
-template<typename CustomUnaryOp, typename Scalar, typename Derived>
-const CwiseUnaryOp<CustomUnaryOp, Derived>
-cwise(const MatrixBase<Scalar, Derived> &mat)
-{
-  return CwiseUnaryOp<CustomUnaryOp, Derived>(mat.asArg());
-}
-
-/** \returns an expression of a custom coefficient-wise unary operator of *this
-  *
-  * The template parameter \a CustomUnaryOp is the template functor
-  * of the custom unary operator.
-  *
-  * \note since cwise is a templated member with a mandatory template parameter,
-  * the keyword template as to be used if the matrix type is also a template parameter:
-  * \code
-  * template <typename MatrixType> void foo(const MatrixType& m) {
-  *   m.template cwise<ScalarAbsOp>();
-  * }
-  * \endcode
+  * 
+  * Here is an example:
+  * \include class_CwiseUnaryOp.cpp
   *
   * \sa class CwiseUnaryOp, class CwiseBinarOp, MatrixBase::operator-, MatrixBase::cwiseAbs
   */
 template<typename Scalar, typename Derived>
 template<typename CustomUnaryOp>
 const CwiseUnaryOp<CustomUnaryOp, Derived>
-MatrixBase<Scalar, Derived>::cwise() const
+MatrixBase<Scalar, Derived>::cwise(const CustomUnaryOp& func) const
 {
-  return CwiseUnaryOp<CustomUnaryOp, Derived>(asArg());
+  return CwiseUnaryOp<CustomUnaryOp, Derived>(asArg(), func);
 }
 
 
@@ -153,13 +135,13 @@ MatrixBase<Scalar, Derived>::cwise() const
   *
   * \sa class CwiseUnaryOp, MatrixBase::conjugate()
   */
-struct ScalarConjugateOp {
-    template<typename Scalar> static Scalar op(const Scalar& a) { return ei_conj(a); }
+struct ScalarConjugateOp EIGEN_EMPTY_STRUCT {
+    template<typename Scalar> Scalar operator() (const Scalar& a) const { return ei_conj(a); }
 };
 
 /** \returns an expression of the complex conjugate of *this.
   *
-  * \sa adjoint(), class Conjugate */
+  * \sa adjoint() */
 template<typename Scalar, typename Derived>
 const CwiseUnaryOp<ScalarConjugateOp, Derived>
 MatrixBase<Scalar, Derived>::conjugate() const
@@ -167,6 +149,77 @@ MatrixBase<Scalar, Derived>::conjugate() const
   return CwiseUnaryOp<ScalarConjugateOp, Derived>(asArg());
 }
 
+/** \brief Template functor to cast a scalar to another
+  *
+  * \sa class CwiseUnaryOp, MatrixBase::cast()
+  */
+template<typename NewType>
+struct ScalarCastOp EIGEN_EMPTY_STRUCT {
+    typedef NewType result_type;
+    template<typename Scalar> NewType operator() (const Scalar& a) const { return static_cast<NewType>(a); }
+};
+
+/** \returns an expression of *this with the \a Scalar type casted to
+  * \a NewScalar.
+  *
+  * The template parameter \a NewScalar is the type we are casting the scalars to.
+  *
+  * Example: \include MatrixBase_cast.cpp
+  * Output: \verbinclude MatrixBase_cast.out
+  *
+  * \sa class CwiseUnaryOp, class ScalarCastOp
+  */
+template<typename Scalar, typename Derived>
+template<typename NewType>
+const CwiseUnaryOp<ScalarCastOp<NewType>, Derived>
+MatrixBase<Scalar, Derived>::cast() const
+{
+  return CwiseUnaryOp<ScalarCastOp<NewType>, Derived>(asArg());
+}
 
 
+/** \brief Template functor to multiply a scalar by a fixed another one
+  *
+  * \sa class CwiseUnaryOp, MatrixBase::operator*, MatrixBase::operator/
+  */
+template<typename Scalar>
+struct ScalarMultipleOp {
+    ScalarMultipleOp(const Scalar& other) : m_other(other) {}
+    Scalar operator() (const Scalar& a) const { return a * m_other; }
+    const Scalar m_other;
+};
+
+/** \relates MatrixBase \sa class ScalarMultipleOp */
+template<typename Scalar, typename Derived>
+const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>
+MatrixBase<Scalar, Derived>::operator*(const Scalar& scalar) const
+{
+  return CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>(asArg(), ScalarMultipleOp<Scalar>(scalar));
+}
+
+/** \relates MatrixBase \sa class ScalarMultipleOp */
+template<typename Scalar, typename Derived>
+const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>
+MatrixBase<Scalar, Derived>::operator/(const Scalar& scalar) const
+{
+  assert(NumTraits<Scalar>::HasFloatingPoint);
+  return CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>(asArg(), ScalarMultipleOp<Scalar>(static_cast<Scalar>(1) / scalar));
+}
+
+/** \sa ScalarMultipleOp */
+template<typename Scalar, typename Derived>
+Derived&
+MatrixBase<Scalar, Derived>::operator*=(const Scalar& other)
+{
+  return *this = *this * other;
+}
+
+/** \sa ScalarMultipleOp */
+template<typename Scalar, typename Derived>
+Derived&
+MatrixBase<Scalar, Derived>::operator/=(const Scalar& other)
+{
+  return *this = *this / other;
+}
+
 #endif // EIGEN_CWISE_UNARY_OP_H
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index 652514fa3..d2ce1803a 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -122,7 +122,7 @@ typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm() const
   * \sa norm()
   */
 template<typename Scalar, typename Derived>
-const ScalarMultiple<Derived>
+const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>
 MatrixBase<Scalar, Derived>::normalized() const
 {
   return (*this) / norm();
diff --git a/Eigen/src/Core/ForwardDeclarations.h b/Eigen/src/Core/ForwardDeclarations.h
index 6ed9dacfe..aa9c4a052 100644
--- a/Eigen/src/Core/ForwardDeclarations.h
+++ b/Eigen/src/Core/ForwardDeclarations.h
@@ -27,7 +27,6 @@
 
 template<typename _Scalar, int _Rows, int _Cols, int _StorageOrder, int _MaxRows, int _MaxCols> class Matrix;
 template<typename MatrixType> class MatrixRef;
-template<typename NewScalar, typename MatrixType> class Cast;
 template<typename MatrixType> class Row;
 template<typename MatrixType> class Column;
 template<typename MatrixType> class Minor;
@@ -37,7 +36,6 @@ template<typename MatrixType> class Conjugate;
 template<typename BinaryOp, typename Lhs, typename Rhs> class CwiseBinaryOp;
 template<typename UnaryOp, typename MatrixType> class CwiseUnaryOp;
 template<typename Lhs, typename Rhs> class Product;
-template<typename MatrixType> class ScalarMultiple;
 template<typename MatrixType> class Random;
 template<typename MatrixType> class Zero;
 template<typename MatrixType> class Ones;
@@ -52,6 +50,8 @@ struct ScalarQuotientOp;
 struct ScalarOppositeOp;
 struct ScalarConjugateOp;
 struct ScalarAbsOp;
+template<typename NewType> struct ScalarCastOp;
+template<typename Scalar>  struct ScalarMultipleOp;
 
 template<typename T> struct Reference
 {
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f5a2ad178..c89d3784c 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -227,7 +227,8 @@ template<typename Scalar, typename Derived> class MatrixBase
 
     /// \name matrix transformation
     //@{
-    template<typename NewScalar> const Cast<NewScalar, Derived> cast() const;
+    template<typename NewType>
+    const CwiseUnaryOp<ScalarCastOp<NewType>, Derived> cast() const;
 
     const DiagonalMatrix<Derived> asDiagonal() const;
 
@@ -237,7 +238,7 @@ template<typename Scalar, typename Derived> class MatrixBase
     const CwiseUnaryOp<ScalarConjugateOp, Derived> conjugate() const;
     const Transpose<CwiseUnaryOp<ScalarConjugateOp, Derived> > adjoint() const;
 
-    const ScalarMultiple<Derived> normalized() const;
+    const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived> normalized() const;
     //@}
 
     // FIXME not sure about the following name
@@ -304,12 +305,11 @@ template<typename Scalar, typename Derived> class MatrixBase
     Derived& operator*=(const Scalar& other);
     Derived& operator/=(const Scalar& other);
 
-    const ScalarMultiple<Derived> operator*(const Scalar& scalar) const;
-    const ScalarMultiple<Derived> operator/(const Scalar& scalar) const;
+    const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived> operator*(const Scalar& scalar) const;
+    const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived> operator/(const Scalar& scalar) const;
 
-    friend
-    const ScalarMultiple<Derived> operator*(const Scalar& scalar,
-                                            const MatrixBase& matrix)
+    friend const CwiseUnaryOp<ScalarMultipleOp<Scalar>, Derived>
+    operator*(const Scalar& scalar, const MatrixBase& matrix)
     { return matrix*scalar; }
 
     template<typename OtherDerived>
@@ -356,11 +356,11 @@ template<typename Scalar, typename Derived> class MatrixBase
     const Eval<Derived> eval() const EIGEN_ALWAYS_INLINE;
 
     template<typename CustomUnaryOp>
-    const CwiseUnaryOp<CustomUnaryOp, Derived> cwise() const;
+    const CwiseUnaryOp<CustomUnaryOp, Derived> cwise(const CustomUnaryOp& func = CustomUnaryOp()) const;
 
     template<typename CustomBinaryOp, typename OtherDerived>
     const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
-    cwise(const MatrixBase<Scalar, OtherDerived> &other) const;
+    cwise(const MatrixBase<Scalar, OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const;
     //@}
 
     /** puts in *row and *col the location of the coefficient of *this
diff --git a/Eigen/src/Core/ScalarMultiple.h b/Eigen/src/Core/ScalarMultiple.h
deleted file mode 100644
index e37865790..000000000
--- a/Eigen/src/Core/ScalarMultiple.h
+++ /dev/null
@@ -1,107 +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) 2006-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_SCALARMULTIPLE_H
-#define EIGEN_SCALARMULTIPLE_H
-
-/** \class ScalarMultiple
-  *
-  * \brief Expression of the product of a matrix or vector by a scalar
-  *
-  * \param FactorTye the type of scalar by which to multiply
-  * \param MatrixType the type of the matrix or vector to multiply
-  *
-  * This class represents an expression of the product of a matrix or vector by a scalar.
-  * It is the return type of the operator* between a matrix or vector and a scalar, and most
-  * of the time this is the only way it is used.
-  */
-template<typename MatrixType> class ScalarMultiple : NoOperatorEquals,
-  public MatrixBase<typename MatrixType::Scalar, ScalarMultiple<MatrixType> >
-{
-  public:
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::AsArg MatRef;
-    friend class MatrixBase<Scalar, ScalarMultiple>;
-    friend class MatrixBase<Scalar, ScalarMultiple>::Traits;
-    typedef MatrixBase<Scalar, ScalarMultiple> Base;
-
-    ScalarMultiple(const MatRef& matrix, Scalar factor)
-      : m_matrix(matrix), m_factor(factor) {}
-
-  private:
-    enum {
-      RowsAtCompileTime = MatrixType::Traits::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::Traits::ColsAtCompileTime,
-      MaxRowsAtCompileTime = MatrixType::Traits::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::Traits::MaxColsAtCompileTime
-    };
-
-    const ScalarMultiple& _asArg() const { return *this; }
-    int _rows() const { return m_matrix.rows(); }
-    int _cols() const { return m_matrix.cols(); }
-
-    Scalar _coeff(int row, int col) const
-    {
-      return m_factor * m_matrix.coeff(row, col);
-    }
-
-  protected:
-    const MatRef m_matrix;
-    const Scalar m_factor;
-};
-
-/** relates MatrixBase sa class ScalarMultiple */
-template<typename Scalar, typename Derived>
-const ScalarMultiple<Derived>
-MatrixBase<Scalar, Derived>::operator*(const Scalar& scalar) const
-{
-  return ScalarMultiple<Derived>(asArg(), scalar);
-}
-
-/** \relates MatrixBase \sa class ScalarMultiple */
-template<typename Scalar, typename Derived>
-const ScalarMultiple<Derived>
-MatrixBase<Scalar, Derived>::operator/(const Scalar& scalar) const
-{
-  assert(NumTraits<Scalar>::HasFloatingPoint);
-  return ScalarMultiple<Derived>(asArg(), static_cast<Scalar>(1) / scalar);
-}
-
-/** \sa ScalarMultiple */
-template<typename Scalar, typename Derived>
-Derived&
-MatrixBase<Scalar, Derived>::operator*=(const Scalar& other)
-{
-  return *this = *this * other;
-}
-
-/** \sa ScalarMultiple */
-template<typename Scalar, typename Derived>
-Derived&
-MatrixBase<Scalar, Derived>::operator/=(const Scalar& other)
-{
-  return *this = *this / other;
-}
-
-#endif // EIGEN_SCALARMULTIPLE_H
diff --git a/Eigen/src/Core/Util.h b/Eigen/src/Core/Util.h
index 9d0aa07a0..962a942d1 100644
--- a/Eigen/src/Core/Util.h
+++ b/Eigen/src/Core/Util.h
@@ -113,6 +113,14 @@ const int RowMajor = 1;
 
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 
+// just a workaround because GCC seems to not really like empty structs
+#ifdef __GNUG__
+  struct EiEmptyStruct{char _ei_dummy_;};
+  #define EIGEN_EMPTY_STRUCT : Eigen::EiEmptyStruct
+#else
+  #define EIGEN_EMPTY_STRUCT
+#endif
+
 //classes inheriting NoOperatorEquals don't generate a default operator=.
 class NoOperatorEquals
 {
@@ -120,7 +128,7 @@ class NoOperatorEquals
     NoOperatorEquals& operator=(const NoOperatorEquals&);
 };
 
-template<int Value> class IntAtRunTimeIfDynamic
+template<int Value> class IntAtRunTimeIfDynamic EIGEN_EMPTY_STRUCT
 {
   public:
     IntAtRunTimeIfDynamic() {}
@@ -139,4 +147,58 @@ template<> class IntAtRunTimeIfDynamic<Dynamic>
     void setValue(int value) { m_value = value; }
 };
 
+struct ei_has_nothing {int a[1];};
+struct ei_has_std_result_type {int a[2];};
+struct ei_has_tr1_result {int a[3];};
+
+/** Convenient struct to get the result type of a unary or binary functor.
+  * 
+  * It supports both the current STL mechanism (using the result_type member) as well as
+  * upcoming next STL generation (using a templated result member).
+  * If none of these member is provided, then the type of the first argument is returned.
+  */
+template<typename T> struct ei_result_of {};
+
+template<typename Func, typename ArgType, int SizeOf=sizeof(ei_has_nothing)>
+struct ei_unary_result_of_select {typedef ArgType type;};
+
+template<typename Func, typename ArgType>
+struct ei_unary_result_of_select<Func, ArgType, sizeof(ei_has_std_result_type)> {typedef typename Func::result_type type;};
+
+template<typename Func, typename ArgType>
+struct ei_unary_result_of_select<Func, ArgType, sizeof(ei_has_tr1_result)> {typedef typename Func::template result<Func(ArgType)>::type type;};
+
+template<typename Func, typename ArgType>
+struct ei_result_of<Func(ArgType)> {
+    template<typename T>
+    static ei_has_std_result_type testFunctor(T const *, typename T::result_type const * = 0);
+    template<typename T>
+    static ei_has_tr1_result      testFunctor(T const *, typename T::template result<T(ArgType)>::type const * = 0);
+    static ei_has_nothing         testFunctor(...);
+    
+    typedef typename ei_unary_result_of_select<Func, ArgType, sizeof(testFunctor(static_cast<Func*>(0)))>::type type;
+};
+
+template<typename Func, typename ArgType0, typename ArgType1, int SizeOf=sizeof(ei_has_nothing)>
+struct ei_binary_result_of_select {typedef ArgType0 type;};
+
+template<typename Func, typename ArgType0, typename ArgType1>
+struct ei_binary_result_of_select<Func, ArgType0, ArgType1, sizeof(ei_has_std_result_type)>
+{typedef typename Func::result_type type;};
+
+template<typename Func, typename ArgType0, typename ArgType1>
+struct ei_binary_result_of_select<Func, ArgType0, ArgType1, sizeof(ei_has_tr1_result)>
+{typedef typename Func::template result<Func(ArgType0,ArgType1)>::type type;};
+
+template<typename Func, typename ArgType0, typename ArgType1>
+struct ei_result_of<Func(ArgType0,ArgType1)> {
+    template<typename T>
+    static ei_has_std_result_type testFunctor(T const *, typename T::result_type const * = 0);
+    template<typename T>
+    static ei_has_tr1_result      testFunctor(T const *, typename T::template result<T(ArgType0,ArgType1)>::type const * = 0);
+    static ei_has_nothing         testFunctor(...);
+    
+    typedef typename ei_binary_result_of_select<Func, ArgType0, ArgType1, sizeof(testFunctor(static_cast<Func*>(0)))>::type type;
+};
+
 #endif // EIGEN_UTIL_H
diff --git a/doc/bench_unrolling b/doc/bench_unrolling
new file mode 100755
index 000000000..4af791412
--- /dev/null
+++ b/doc/bench_unrolling
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+# gcc : CXX="g++  -finline-limit=10000 -ftemplate-depth-2000 --param max-inline-recursive-depth=2000"
+# icc : CXX="icpc -fast -no-inline-max-size -fno-exceptions"
+
+for ((i=1; i<16; ++i)); do
+    echo "Matrix size: $i x $i :"
+    $CXX -O3 -I.. -DNDEBUG  benchmark.cpp -DMATSIZE=$i -DEIGEN_UNROLLING_LIMIT_OPEQUAL=1024 -DEIGEN_UNROLLING_LIMIT_PRODUCT=25 -o benchmark && time ./benchmark >/dev/null
+    $CXX -O3 -I.. -DNDEBUG -finline-limit=10000 benchmark.cpp -DMATSIZE=$i -DEIGEN_DONT_USE_UNROLLED_LOOPS=1 -o benchmark && time ./benchmark >/dev/null
+    echo " "
+done
diff --git a/doc/benchmark.cpp b/doc/benchmark.cpp
index ad515e04e..0d95a5043 100644
--- a/doc/benchmark.cpp
+++ b/doc/benchmark.cpp
@@ -1,24 +1,31 @@
-// g++ -O3 -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark
+// g++ -O3 -DNDEBUG -DMATSIZE=<x> benchmark.cpp -o benchmark && time ./benchmark
 #include <cstdlib>
 #include <cmath>
 #include <Eigen/Core>
 
-//using namespace std;
+#ifndef MATSIZE
+#define MATSIZE 3
+#endif
+
+using namespace std;
 USING_PART_OF_NAMESPACE_EIGEN
 
 int main(int argc, char *argv[])
 {
-	Matrix3d I;
-	Matrix3d m;
-	for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
-	{
-		I(i,j) = (i==j);
-		m(i,j) = (i+3*j);
-	}
-	for(int a = 0; a < 400000000; a++)
-	{
-		m = I + 0.00005 * (m + m*m);
-	}
-	std::cout << m << std::endl;
-	return 0;
+    Matrix<double,MATSIZE,MATSIZE> I;
+    Matrix<double,MATSIZE,MATSIZE> m;
+    for(int i = 0; i < MATSIZE; i++)
+        for(int j = 0; j < MATSIZE; j++)
+        {
+            I(i,j) = (i==j);
+            m(i,j) = (i+MATSIZE*j);
+        }
+    asm("#begin");
+    for(int a = 0; a < 40000000; a++)
+    {
+        m = I + 0.00005 * (m + m*m);
+    }
+    asm("#end");
+    cout << m << endl;
+    return 0;
 }
diff --git a/doc/benchmarkX.cpp b/doc/benchmarkX.cpp
index 1670341e9..1c8b6b803 100644
--- a/doc/benchmarkX.cpp
+++ b/doc/benchmarkX.cpp
@@ -7,7 +7,7 @@ USING_PART_OF_NAMESPACE_EIGEN
 
 int main(int argc, char *argv[])
 {
-	MatrixXd I = MatrixXd::identity(20);
+	MatrixXd I = MatrixXd::identity(20,20);
 	MatrixXd m(20,20);
 	for(int i = 0; i < 20; i++) for(int j = 0; j < 20; j++)
 	{
diff --git a/doc/benchmark_suite b/doc/benchmark_suite
index e0969d0eb..9ddfccbf6 100755
--- a/doc/benchmark_suite
+++ b/doc/benchmark_suite
@@ -1,16 +1,17 @@
+#!/bin/bash
 echo "Fixed size 3x3, ColumnMajor, -DNDEBUG"
-g++ -O3 -I .. -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+$CXX -O3 -I .. -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
 echo "Fixed size 3x3, ColumnMajor, with asserts"
-g++ -O3 -I .. benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+$CXX -O3 -I .. benchmark.cpp -o benchmark && time ./benchmark >/dev/null
 echo "Fixed size 3x3, RowMajor, -DNDEBUG"
-g++ -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+$CXX -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
 echo "Fixed size 3x3, RowMajor, with asserts"
-g++ -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+$CXX -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR benchmark.cpp -o benchmark && time ./benchmark >/dev/null
 echo "Dynamic size 20x20, ColumnMajor, -DNDEBUG"
-g++ -O3 -I .. -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+$CXX -O3 -I .. -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
 echo "Dynamic size 20x20, ColumnMajor, with asserts"
-g++ -O3 -I .. benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+$CXX -O3 -I .. benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
 echo "Dynamic size 20x20, RowMajor, -DNDEBUG"
-g++ -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+$CXX -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
 echo "Dynamic size 20x20, RowMajor, with asserts"
-g++ -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+$CXX -O3 -I .. -DEIGEN_DEFAULT_TO_ROW_MAJOR benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
diff --git a/doc/examples/class_Cast.cpp b/doc/examples/class_Cast.cpp
index 9751b05c1..88549732f 100644
--- a/doc/examples/class_Cast.cpp
+++ b/doc/examples/class_Cast.cpp
@@ -3,16 +3,13 @@ USING_PART_OF_NAMESPACE_EIGEN
 using namespace std;
 
 template<typename Scalar, typename Derived>
-const Eigen::Cast<
-  typename Eigen::NumTraits<Scalar>::FloatingPoint,
+const Eigen::CwiseUnaryOp<
+  Eigen::ScalarCastOp<typename Eigen::NumTraits<Scalar>::FloatingPoint>,
   Derived
 >
 castToFloatingPoint(const MatrixBase<Scalar, Derived>& m)
 {
-  return Eigen::Cast<
-    typename Eigen::NumTraits<Scalar>::FloatingPoint,
-    Derived
-  >(m.asArg());
+  return m.template cast<typename Eigen::NumTraits<Scalar>::FloatingPoint>();
 }
 
 int main(int, char**)
diff --git a/doc/examples/class_CwiseBinaryOp.cpp b/doc/examples/class_CwiseBinaryOp.cpp
index 9bb10cdcb..f5439a434 100644
--- a/doc/examples/class_CwiseBinaryOp.cpp
+++ b/doc/examples/class_CwiseBinaryOp.cpp
@@ -3,9 +3,9 @@ USING_PART_OF_NAMESPACE_EIGEN
 using namespace std;
 
 // define a custom template binary functor
-struct CwiseMinOp {
+struct CwiseMinOp EIGEN_EMPTY_STRUCT {
     template<typename Scalar>
-    static Scalar op(const Scalar& a, const Scalar& b) { return std::min(a,b); }
+    Scalar operator()(const Scalar& a, const Scalar& b) const { return std::min(a,b); }
 };
 
 // define a custom binary operator between two matrices
@@ -14,14 +14,13 @@ const Eigen::CwiseBinaryOp<CwiseMinOp, Derived1, Derived2>
 cwiseMin(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
 {
   return Eigen::CwiseBinaryOp<CwiseMinOp, Derived1, Derived2>(mat1.asArg(), mat2.asArg());
-  // Note that the above is equivalent to:
-  // return mat1.template cwise<CwiseMinOp>(mat2);
 }
 
 int main(int, char**)
 {
   Matrix4d m1 = Matrix4d::random(), m2 = Matrix4d::random();
-  cout << cwiseMin(m1,m2) << endl;          // use our new global operator
-  cout << m1.cwise<CwiseMinOp>(m2) << endl; // directly use the generic expression member
+  cout << cwiseMin(m1,m2) << endl;            // use our new global operator
+  cout << m1.cwise<CwiseMinOp>(m2) << endl;   // directly use the generic expression member
+  cout << m1.cwise(m2, CwiseMinOp()) << endl; // directly use the generic expression member (variant)
   return 0;
 }
diff --git a/doc/examples/class_CwiseUnaryOp.cpp b/doc/examples/class_CwiseUnaryOp.cpp
new file mode 100644
index 000000000..e10037352
--- /dev/null
+++ b/doc/examples/class_CwiseUnaryOp.cpp
@@ -0,0 +1,17 @@
+#include <Eigen/Core>
+USING_PART_OF_NAMESPACE_EIGEN
+using namespace std;
+
+// define a custom template binary functor
+template<typename Scalar>
+struct CwiseClampOp EIGEN_EMPTY_STRUCT {
+    CwiseClampOp(const Scalar& inf, const Scalar& sup) : m_inf(inf), m_sup(sup) {}
+    Scalar operator()(const Scalar& x) const { return x<m_inf ? m_inf : (x>m_sup : m_sup : x); }
+};
+
+int main(int, char**)
+{
+  Matrix4d m1 = Matrix4d::random(), m2 = Matrix4d::random();
+  cout << m1.cwise(m2, CwiseClampOp<Matrix4d::Scalar>(-0.5,0.5)) << endl;
+  return 0;
+}
diff --git a/test/cwiseop.cpp b/test/cwiseop.cpp
index 7c4a9a364..02a0a1afd 100644
--- a/test/cwiseop.cpp
+++ b/test/cwiseop.cpp
@@ -31,7 +31,7 @@
 namespace Eigen {
 
 struct AddIfNull {
-    template<typename Scalar> static Scalar op(const Scalar a, const Scalar b) {return a<=1e-3 ? b : a;}
+    template<typename Scalar> Scalar operator() (const Scalar a, const Scalar b) const {return a<=1e-3 ? b : a;}
 };
 
 template<typename MatrixType> void cwiseops(const MatrixType& m)