diff --git a/Eigen/Core b/Eigen/Core
index ace3a204e..792fc0135 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -159,7 +159,6 @@ namespace Eigen {
 #include "src/Core/CwiseNullaryOp.h"
 #include "src/Core/CwiseUnaryView.h"
 #include "src/Core/Dot.h"
-#include "src/Core/DiagonalProduct.h"
 #include "src/Core/SolveTriangular.h"
 #include "src/Core/MapBase.h"
 #include "src/Core/Map.h"
@@ -168,6 +167,7 @@ namespace Eigen {
 #include "src/Core/Transpose.h"
 #include "src/Core/DiagonalMatrix.h"
 #include "src/Core/Diagonal.h"
+#include "src/Core/DiagonalProduct.h"
 #include "src/Core/Redux.h"
 #include "src/Core/Visitor.h"
 #include "src/Core/Fuzzy.h"
diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h
index 1258ae051..87e7a1337 100644
--- a/Eigen/src/Core/Assign.h
+++ b/Eigen/src/Core/Assign.h
@@ -436,10 +436,16 @@ struct ei_assign_selector<Derived,OtherDerived,true,true> {
 
 template<typename Derived>
 template<typename OtherDerived>
-EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>
-  ::operator=(const MatrixBase<OtherDerived>& other)
+EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase<OtherDerived>& other)
 {
   return ei_assign_selector<Derived,OtherDerived>::run(derived(), other.derived());
 }
 
+template<typename Derived>
+EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase& other)
+{
+  return ei_assign_selector<Derived,Derived>::run(derived(), other.derived());
+}
+
+
 #endif // EIGEN_ASSIGN_H
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 2094913bc..413ad5c9f 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -26,250 +26,164 @@
 #ifndef EIGEN_DIAGONALMATRIX_H
 #define EIGEN_DIAGONALMATRIX_H
 
-
-template<typename CoeffsVectorType, typename Derived>
-class DiagonalMatrixBase : ei_no_assignment_operator,
-   public MatrixBase<Derived>
+template<typename Derived>
+class DiagonalBase : public MultiplierBase<Derived>
 {
   public:
-    typedef MatrixBase<Derived> Base;
-    typedef typename ei_traits<Derived>::Scalar Scalar;
-    typedef typename Base::PacketScalar PacketScalar;
-    using Base::derived;
-    typedef typename ei_cleantype<CoeffsVectorType>::type _CoeffsVectorType;
-
-  protected:
-
-    // MSVC gets crazy if we define default parameters
-    template<typename OtherDerived, bool IsVector, bool IsDiagonal> struct construct_from_expression;
-
-    // = vector
-    template<typename OtherDerived>
-    struct construct_from_expression<OtherDerived,true,false>
-    {
-      static void run(Derived& dst, const OtherDerived& src)
-      { dst.diagonal() = src; }
-    };
-
-    // = diagonal expression
-    template<typename OtherDerived, bool IsVector>
-    struct construct_from_expression<OtherDerived,IsVector,true>
-    {
-      static void run(Derived& dst, const OtherDerived& src)
-      { dst.diagonal() = src.diagonal(); }
-    };
-
-    /** Default constructor without initialization */
-    inline DiagonalMatrixBase() {}
-    /** Constructs a diagonal matrix with given dimension */
-    inline DiagonalMatrixBase(int dim) : m_coeffs(dim) {}
-    /** Generic constructor from an expression */
-    template<typename OtherDerived>
-    inline DiagonalMatrixBase(const MatrixBase<OtherDerived>& other)
-    {
-      construct_from_expression<OtherDerived,OtherDerived::IsVectorAtCompileTime,ei_is_diagonal<OtherDerived>::ret>
-        ::run(derived(),other.derived());
-    }
+    typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
+    typedef typename DiagonalVectorType::Scalar Scalar;
     
-    template<typename NewType,int dummy=0>
-    struct cast_selector {
-      typedef const DiagonalMatrixWrapper<NestByValue<CwiseUnaryOp<ei_scalar_cast_op<Scalar, NewType>, _CoeffsVectorType> > > return_type;
-      inline static return_type run(const DiagonalMatrixBase& d) {
-        return d.m_coeffs.template cast<NewType>().nestByValue().asDiagonal();
-      }
+    enum {
+      RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+      ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+      MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
+      MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
+      IsVectorAtCompileTime = 0
     };
     
-    template<int dummy>
-    struct cast_selector<Scalar,dummy> {
-      typedef const Derived& return_type;
-      inline static return_type run(const DiagonalMatrixBase& d) {
-        return d.derived();
-      }
-    };
+    typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
+        
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
+    inline Derived& derived() { return *static_cast<Derived*>(this); }
+    #endif // not EIGEN_PARSED_BY_DOXYGEN
+    
+    DenseMatrixType toDenseMatrix() const { return derived(); }
 
-  public:
-
-    inline DiagonalMatrixBase(const _CoeffsVectorType& coeffs) : m_coeffs(coeffs)
-    {
-      EIGEN_STATIC_ASSERT_VECTOR_ONLY(_CoeffsVectorType);
-      ei_assert(coeffs.size() > 0);
-    }
-
-    template<typename NewType>
-    inline typename cast_selector<NewType,0>::return_type
-    cast() const
-    {
-      return cast_selector<NewType,0>::run(*this);
-    }
-
-    /** Assignment operator.
-      * The right-hand-side \a other must be either a vector representing the diagonal
-      * coefficients or a diagonal matrix expression.
-      */
-    template<typename OtherDerived>
-    inline Derived& operator=(const MatrixBase<OtherDerived>& other)
-    {
-      construct_from_expression<OtherDerived,OtherDerived::IsVectorAtCompileTime,ei_is_diagonal<OtherDerived>::ret>
-        ::run(derived(),other);
-      return derived();
-    }
-
-    inline int rows() const { return m_coeffs.size(); }
-    inline int cols() const { return m_coeffs.size(); }
-
-    inline const Scalar coeff(int row, int col) const
-    {
-      return row == col ? m_coeffs.coeff(row) : static_cast<Scalar>(0);
-    }
-
-    inline Scalar& coeffRef(int row, int col)
-    {
-      ei_assert(row==col);
-      return m_coeffs.coeffRef(row);
-    }
-
-    inline _CoeffsVectorType& diagonal() { return m_coeffs; }
-    inline const _CoeffsVectorType& diagonal() const { return m_coeffs; }
-
-  protected:
-    CoeffsVectorType m_coeffs;
+    inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
+    inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
+    
+    template<typename MatrixDerived>
+    const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft>
+    operator*(const MatrixBase<MatrixDerived> &matrix) const;
 };
 
+template<typename Derived>
+template<typename DiagonalDerived>
+Derived& MatrixBase<Derived>::operator=(const DiagonalBase<DiagonalDerived> &other)
+{
+  setZero();
+  diagonal() = other.diagonal();
+  return derived();
+}
+
 /** \class DiagonalMatrix
   * \nonstableyet
   *
-  * \brief Represent a diagonal matrix with its storage
+  * \brief Represents a diagonal matrix with its storage
   *
   * \param _Scalar the type of coefficients
-  * \param _Size the dimension of the matrix
+  * \param _Size the dimension of the matrix, or Dynamic
   *
   * \sa class Matrix
   */
-template<typename _Scalar,int _Size>
-struct ei_traits<DiagonalMatrix<_Scalar,_Size> > : ei_traits<Matrix<_Scalar,_Size,_Size> >
+template<typename _Scalar, int _Size>
+struct ei_traits<DiagonalMatrix<_Scalar,_Size> >
 {
-  enum {
-    Flags = (ei_traits<Matrix<_Scalar,_Size,_Size> >::Flags & HereditaryBits) | DiagonalBits
-  };
+  typedef Matrix<_Scalar,_Size,1> DiagonalVectorType;
 };
 
 template<typename _Scalar, int _Size>
 class DiagonalMatrix
-  : public DiagonalMatrixBase<Matrix<_Scalar,_Size,1>, DiagonalMatrix<_Scalar,_Size> >
+  : public DiagonalBase<DiagonalMatrix<_Scalar,_Size> >
 {
   public:
-    EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrix)
-    typedef DiagonalMatrixBase<Matrix<_Scalar,_Size,1>, DiagonalMatrix<_Scalar,_Size> > DiagonalBase;
-
+    
+    typedef typename ei_traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType;
+    typedef const DiagonalMatrix& Nested;
+    typedef _Scalar Scalar;
+  
   protected:
-    typedef Matrix<_Scalar,_Size,1> CoeffVectorType;
-    using DiagonalBase::m_coeffs;
 
+    DiagonalVectorType m_diagonal;
+    
   public:
 
+    inline const DiagonalVectorType& diagonal() const { return m_diagonal; }
+    inline DiagonalVectorType& diagonal() { return m_diagonal; }
+    
     /** Default constructor without initialization */
-    inline DiagonalMatrix() : DiagonalBase()
-    {}
+    inline DiagonalMatrix() {}
 
     /** Constructs a diagonal matrix with given dimension  */
-    inline DiagonalMatrix(int dim) : DiagonalBase(dim)
-    {}
+    inline DiagonalMatrix(int dim) : m_diagonal(dim) {}
 
     /** 2D only */
-    inline DiagonalMatrix(const Scalar& sx, const Scalar& sy)
-    {
-      EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DiagonalMatrix,2,2);
-      m_coeffs.x() = sx;
-      m_coeffs.y() = sy;
-    }
+    inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {}
+
     /** 3D only */
-    inline DiagonalMatrix(const Scalar& sx, const Scalar& sy, const Scalar& sz)
+    inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {}
+
+    template<typename OtherDerived>
+    inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {}
+
+    /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
+    inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {}
+
+    /** generic constructor from expression of the diagonal coefficients */
+    template<typename OtherDerived>
+    explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other)
+    {}
+
+    template<typename OtherDerived>
+    DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other)
     {
-      EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DiagonalMatrix,3,3);
-      m_coeffs.x() = sx;
-      m_coeffs.y() = sy;
-      m_coeffs.z() = sz;
+      m_diagonal = other.diagonal();
+      return *this;
     }
 
-    /** copy constructor */
-    inline DiagonalMatrix(const DiagonalMatrix& other) : DiagonalBase(other.m_coeffs)
-    {}
-
-    /** generic constructor from expression */
-    template<typename OtherDerived>
-    explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : DiagonalBase(other)
-    {}
-
+    /** This is a special case of the templated operator=. Its purpose is to
+      * prevent a default operator= from hiding the templated operator=.
+      */
     DiagonalMatrix& operator=(const DiagonalMatrix& other)
     {
-      m_coeffs = other.m_coeffs;
-      return *this;
-    }
-
-    template<typename OtherDerived>
-    DiagonalMatrix& operator=(const MatrixBase<OtherDerived>& other)
-    {
-      EIGEN_STATIC_ASSERT(ei_is_diagonal<OtherDerived>::ret, THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX);
-      m_coeffs = other.diagonal();
+      m_diagonal = other.m_diagonal();
       return *this;
     }
     
-    inline void resize(int size)
-    {
-      m_coeffs.resize(size);
-    }
-    
-    inline void resize(int rows, int cols)
-    {
-      ei_assert(rows==cols && "a diagonal matrix must be square");
-      m_coeffs.resize(rows);
-    }
-    
-    inline void setZero() { m_coeffs.setZero(); }
+    inline void resize(int size) { m_diagonal.resize(size); }
+    inline void setZero() { m_diagonal.setZero(); }
+    inline void setZero(int size) { m_diagonal.setZero(size); }
+    inline void setIdentity() { m_diagonal.setIdentity(); }
+    inline void setIdentity(int size) { m_diagonal.setIdentity(size); }
 };
 
-/** \class DiagonalMatrixWrapper
+/** \class DiagonalWrapper
   * \nonstableyet
   *
   * \brief Expression of a diagonal matrix
   *
-  * \param CoeffsVectorType the type of the vector of diagonal coefficients
+  * \param _DiagonalVectorType the type of the vector of diagonal coefficients
   *
   * This class is an expression of a diagonal matrix with given vector of diagonal
-  * coefficients. It is the return type of MatrixBase::diagonal(const OtherDerived&)
-  * and most of the time this is the only way it is used.
+  * coefficients. It is the return type of MatrixBase::asDiagonal()
+  * and most of the time this is the only way that it is used.
   *
-  * \sa class DiagonalMatrixBase, class DiagonalMatrix, MatrixBase::asDiagonal()
+  * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal()
   */
-template<typename CoeffsVectorType>
-struct ei_traits<DiagonalMatrixWrapper<CoeffsVectorType> >
+template<typename _DiagonalVectorType>
+struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
 {
-  typedef typename CoeffsVectorType::Scalar Scalar;
-  typedef typename ei_nested<CoeffsVectorType>::type CoeffsVectorTypeNested;
-  typedef typename ei_unref<CoeffsVectorTypeNested>::type _CoeffsVectorTypeNested;
-  enum {
-    RowsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
-    ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
-    MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
-    MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
-    Flags = (_CoeffsVectorTypeNested::Flags & HereditaryBits) | DiagonalBits,
-    CoeffReadCost = _CoeffsVectorTypeNested::CoeffReadCost
-  };
+  typedef _DiagonalVectorType DiagonalVectorType;
 };
-template<typename CoeffsVectorType>
-class DiagonalMatrixWrapper
-  : public DiagonalMatrixBase<typename CoeffsVectorType::Nested, DiagonalMatrixWrapper<CoeffsVectorType> >
+
+template<typename _DiagonalVectorType>
+class DiagonalWrapper
+  : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, ei_no_assignment_operator
 {
-    typedef typename CoeffsVectorType::Nested CoeffsVectorTypeNested;
-    typedef DiagonalMatrixBase<CoeffsVectorTypeNested, DiagonalMatrixWrapper<CoeffsVectorType> > DiagonalBase;
   public:
-    EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrixWrapper)
-    inline DiagonalMatrixWrapper(const CoeffsVectorType& coeffs) : DiagonalBase(coeffs)
-    {}
+    typedef _DiagonalVectorType DiagonalVectorType;
+    typedef DiagonalWrapper Nested;
+    
+    inline DiagonalWrapper(const DiagonalVectorType& diagonal) : m_diagonal(diagonal) {}
+    const DiagonalVectorType& diagonal() const { return m_diagonal; }
+    
+  protected:
+    const typename DiagonalVectorType::Nested m_diagonal;
 };
 
 /** \nonstableyet
-  * \returns an expression of a diagonal matrix with *this as vector of diagonal coefficients
+  * \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients
   *
   * \only_for_vectors
   *
@@ -278,10 +192,10 @@ class DiagonalMatrixWrapper
   * Example: \include MatrixBase_asDiagonal.cpp
   * Output: \verbinclude MatrixBase_asDiagonal.out
   *
-  * \sa class DiagonalMatrix, isDiagonal()
+  * \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal()
   **/
 template<typename Derived>
-inline const DiagonalMatrixWrapper<Derived>
+inline const DiagonalWrapper<Derived>
 MatrixBase<Derived>::asDiagonal() const
 {
   return derived();
diff --git a/Eigen/src/Core/DiagonalProduct.h b/Eigen/src/Core/DiagonalProduct.h
index f0e0d142b..5948111c6 100644
--- a/Eigen/src/Core/DiagonalProduct.h
+++ b/Eigen/src/Core/DiagonalProduct.h
@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -26,105 +26,66 @@
 #ifndef EIGEN_DIAGONALPRODUCT_H
 #define EIGEN_DIAGONALPRODUCT_H
 
-/** \internal Specialization of ei_nested for DiagonalMatrix.
- *  Unlike ei_nested, if the argument is a DiagonalMatrix and if it must be evaluated,
- *  then it evaluated to a DiagonalMatrix having its own argument evaluated.
- */
-template<typename T, int N, bool IsDiagonal = ei_is_diagonal<T>::ret> struct ei_nested_diagonal : ei_nested<T,N> {};
-template<typename T, int N> struct ei_nested_diagonal<T,N,true>
- : ei_nested<T, N, DiagonalMatrix<typename T::Scalar, EIGEN_ENUM_MIN(T::RowsAtCompileTime,T::ColsAtCompileTime)> >
-{};
-
-// specialization of ProductReturnType
-template<typename Lhs, typename Rhs>
-struct ProductReturnType<Lhs,Rhs,DiagonalProduct>
+template<typename MatrixType, typename DiagonalType, int Order>
+struct ei_traits<DiagonalProduct<MatrixType, DiagonalType, Order> >
 {
-  typedef typename ei_nested_diagonal<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
-  typedef typename ei_nested_diagonal<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
-
-  typedef Product<LhsNested, RhsNested, DiagonalProduct> Type;
-};
-
-template<typename LhsNested, typename RhsNested>
-struct ei_traits<Product<LhsNested, RhsNested, DiagonalProduct> >
-{
-  // clean the nested types:
-  typedef typename ei_cleantype<LhsNested>::type _LhsNested;
-  typedef typename ei_cleantype<RhsNested>::type _RhsNested;
-  typedef typename _LhsNested::Scalar Scalar;
-
+  typedef typename MatrixType::Scalar Scalar;
   enum {
-    LhsFlags = _LhsNested::Flags,
-    RhsFlags = _RhsNested::Flags,
-    RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
-    ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
-    MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
-    MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
-
-    LhsIsDiagonal = ei_is_diagonal<_LhsNested>::ret,
-    RhsIsDiagonal = ei_is_diagonal<_RhsNested>::ret,
-
-    CanVectorizeRhs =  (!RhsIsDiagonal) && (RhsFlags & RowMajorBit) && (RhsFlags & PacketAccessBit)
-                     && (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
-
-    CanVectorizeLhs =  (!LhsIsDiagonal) && (!(LhsFlags & RowMajorBit)) && (LhsFlags & PacketAccessBit)
-                     && (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
-
-    RemovedBits = ~((RhsFlags & RowMajorBit) && (!CanVectorizeLhs) ? 0 : RowMajorBit),
-
-    Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
-          | (((CanVectorizeLhs&&RhsIsDiagonal) || (CanVectorizeRhs&&LhsIsDiagonal)) ? PacketAccessBit : 0),
-
-    CoeffReadCost = NumTraits<Scalar>::MulCost + _LhsNested::CoeffReadCost + _RhsNested::CoeffReadCost
+    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
+    Flags = (unsigned int)(MatrixType::Flags) & HereditaryBits,
+    CoeffReadCost = NumTraits<Scalar>::MulCost + MatrixType::CoeffReadCost + DiagonalType::DiagonalVectorType::CoeffReadCost
   };
 };
 
-template<typename LhsNested, typename RhsNested> class Product<LhsNested, RhsNested, DiagonalProduct> : ei_no_assignment_operator,
-  public MatrixBase<Product<LhsNested, RhsNested, DiagonalProduct> >
+template<typename MatrixType, typename DiagonalType, int Order>
+class DiagonalProduct : ei_no_assignment_operator,
+                        public MatrixBase<DiagonalProduct<MatrixType, DiagonalType, Order> >
 {
-    typedef typename ei_traits<Product>::_LhsNested _LhsNested;
-    typedef typename ei_traits<Product>::_RhsNested _RhsNested;
-
-    enum {
-      RhsIsDiagonal = ei_is_diagonal<_RhsNested>::ret
-    };
-
   public:
 
-    EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
+    EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalProduct)
 
-    template<typename Lhs, typename Rhs>
-    inline Product(const Lhs& lhs, const Rhs& rhs)
-      : m_lhs(lhs), m_rhs(rhs)
+    inline DiagonalProduct(const MatrixType& matrix, const DiagonalType& diagonal)
+      : m_matrix(matrix), m_diagonal(diagonal)
     {
-      ei_assert(lhs.cols() == rhs.rows());
+      ei_assert(diagonal.diagonal().size() == (Order == DiagonalOnTheLeft ? matrix.rows() : matrix.cols()));
     }
 
-    inline int rows() const { return m_lhs.rows(); }
-    inline int cols() const { return m_rhs.cols(); }
+    inline int rows() const { return m_matrix.rows(); }
+    inline int cols() const { return m_matrix.cols(); }
 
     const Scalar coeff(int row, int col) const
     {
-      const int unique = RhsIsDiagonal ? col : row;
-      return m_lhs.coeff(row, unique) * m_rhs.coeff(unique, col);
-    }
-
-    template<int LoadMode>
-    const PacketScalar packet(int row, int col) const
-    {
-      if (RhsIsDiagonal)
-      {
-        return ei_pmul(m_lhs.template packet<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
-      }
-      else
-      {
-        return ei_pmul(ei_pset1(m_lhs.coeff(row, row)), m_rhs.template packet<LoadMode>(row, col));
-      }
+      return m_diagonal.diagonal().coeff(Order == DiagonalOnTheLeft ? row : col) * m_matrix.coeff(row, col);
     }
 
   protected:
-    const LhsNested m_lhs;
-    const RhsNested m_rhs;
+    const typename MatrixType::Nested m_matrix;
+    const typename DiagonalType::Nested m_diagonal;
 };
 
+/** \returns the diagonal matrix product of \c *this by the diagonal matrix \a diagonal.
+  */
+template<typename Derived>
+template<typename DiagonalDerived>
+inline const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
+MatrixBase<Derived>::operator*(const DiagonalBase<DiagonalDerived> &diagonal) const
+{
+  return DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>(derived(), diagonal.derived());
+}
+
+/** \returns the diagonal matrix product of \c *this by the matrix \a matrix.
+  */
+template<typename DiagonalDerived>
+template<typename MatrixDerived>
+inline const DiagonalProduct<MatrixDerived, DiagonalDerived, DiagonalOnTheLeft>
+DiagonalBase<DiagonalDerived>::operator*(const MatrixBase<MatrixDerived> &matrix) const
+{
+  return DiagonalProduct<MatrixDerived, DiagonalDerived, DiagonalOnTheLeft>(matrix.derived(), derived());
+}
+
+
 #endif // EIGEN_DIAGONALPRODUCT_H
diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h
index 126974918..ec9ab3c5e 100644
--- a/Eigen/src/Core/MapBase.h
+++ b/Eigen/src/Core/MapBase.h
@@ -171,11 +171,7 @@ template<typename Derived> class MapBase
       return Base::operator=(other);
     }
 
-    template<typename OtherDerived>
-    Derived& operator=(const MatrixBase<OtherDerived>& other)
-    {
-      return Base::operator=(other);
-    }
+    using Base::operator=;
 
     using Base::operator*=;
 
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index e099ba1e7..8cf5bc64b 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -124,6 +124,7 @@ class Matrix
 {
   public:
     EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix)
+    
     enum { Options = _Options };
     friend class Eigen::Map<Matrix, Unaligned>;
     typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType;
@@ -335,11 +336,11 @@ class Matrix
     EIGEN_STRONG_INLINE Matrix& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func)
     { return Base::operator=(func); }
 
-    EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Matrix, +=)
-    EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Matrix, -=)
-    EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, *=)
-    EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, /=)
-
+    using Base::operator +=;
+    using Base::operator -=;
+    using Base::operator *=;
+    using Base::operator /=;
+    
     /** Default constructor.
       *
       * For fixed-size matrices, does nothing.
@@ -438,6 +439,22 @@ class Matrix
     { other.evalTo(*this); }
     /** Destructor */
     inline ~Matrix() {}
+    
+
+    template<typename DiagonalDerived>
+    EIGEN_STRONG_INLINE Matrix& operator=(const DiagonalBase<DiagonalDerived> &other)
+    {
+      resize(other.diagonal().size(), other.diagonal().size());
+      Base::operator=(other);
+      return *this;
+    }
+
+    template<typename DiagonalDerived>
+    EIGEN_STRONG_INLINE Matrix(const DiagonalBase<DiagonalDerived> &other)
+      : m_storage(other.diagonal().size() * other.diagonal().size(), other.diagonal().size(), other.diagonal().size())
+    {
+      *this = other;
+    }
 
     /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
       * data pointers.
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index addf59570..b0ff09531 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -26,6 +26,32 @@
 #ifndef EIGEN_MATRIXBASE_H
 #define EIGEN_MATRIXBASE_H
 
+
+/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T).
+  *
+  * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase.
+  *
+  * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc.
+  *
+  * Notice that this class is trivial, it is only used to disambiguate overloaded functions.
+  */
+template<typename Derived> struct AnyMatrixBase
+{
+  Derived& derived() { return *static_cast<Derived*>(this); }
+  const Derived& derived() const { return *static_cast<const Derived*>(this); }
+};
+/** Common base class for all classes T such that there are overloaded operator* allowing to
+  * multiply a MatrixBase by a T on both sides.
+  *
+  * In other words, an AnyMatrixBase object is an object that can be multiplied a MatrixBase, the result being again a MatrixBase.
+  *
+  * Besides MatrixBase-derived classes, this also includes certain special matrix classes, such as diagonal matrices.
+  */
+template<typename Derived> struct MultiplierBase : public AnyMatrixBase<Derived>
+{
+  using AnyMatrixBase<Derived>::derived;
+};
+
 /** \class MatrixBase
   *
   * \brief Base class for all matrices, vectors, and expressions
@@ -50,11 +76,11 @@
       cout << x.row(0) << endl;
     }
   * \endcode
-  *
   */
 template<typename Derived> class MatrixBase
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-  : public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
+  : public MultiplierBase<Derived>,
+    public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
                 typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 {
@@ -256,10 +282,7 @@ template<typename Derived> class MatrixBase
     /** Special case of the template operator=, in order to prevent the compiler
       * from generating a default operator= (issue hit with g++ 4.1)
       */
-    inline Derived& operator=(const MatrixBase& other)
-    {
-      return this->operator=<Derived>(other);
-    }
+    Derived& operator=(const MatrixBase& other);
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** Copies \a other into *this without evaluating other. \returns a reference to *this. */
@@ -358,17 +381,30 @@ template<typename Derived> class MatrixBase
     operator*(const Scalar& scalar, const MatrixBase& matrix)
     { return matrix*scalar; }
 
-
     template<typename OtherDerived>
     const typename ProductReturnType<Derived,OtherDerived>::Type
     operator*(const MatrixBase<OtherDerived> &other) const;
 
+    /** replaces \c *this by \c *this * \a other.
+      *
+      * \returns a reference to \c *this
+      */
     template<typename OtherDerived>
-    Derived& operator*=(const MatrixBase<OtherDerived>& other);
+    Derived& operator*=(const MultiplierBase<OtherDerived>& other)
+    {
+      return *this = *this * other.derived();
+    }
+    
+    template<typename DiagonalDerived>
+    const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
+    operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
+
+    template<typename DiagonalDerived>
+    Derived& operator=(const DiagonalBase<DiagonalDerived> &other);
 
     template<typename OtherDerived>
     typename ei_plain_matrix_type_column_major<OtherDerived>::type
-		solveTriangular(const MatrixBase<OtherDerived>& other) const;
+    solveTriangular(const MatrixBase<OtherDerived>& other) const;
 
     template<typename OtherDerived>
     void solveTriangularInPlace(const MatrixBase<OtherDerived>& other) const;
@@ -477,7 +513,7 @@ template<typename Derived> class MatrixBase
     static const BasisReturnType UnitZ();
     static const BasisReturnType UnitW();
 
-    const DiagonalMatrixWrapper<Derived> asDiagonal() const;
+    const DiagonalWrapper<Derived> asDiagonal() const;
 
     void fill(const Scalar& value);
     Derived& setConstant(const Scalar& value);
@@ -588,8 +624,7 @@ template<typename Derived> class MatrixBase
     void visit(Visitor& func) const;
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-    inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
-    inline Derived& derived() { return *static_cast<Derived*>(this); }
+    using MultiplierBase<Derived>::derived;
     inline Derived& const_cast_derived() const
     { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
 #endif // not EIGEN_PARSED_BY_DOXYGEN
@@ -685,16 +720,16 @@ template<typename Derived> class MatrixBase
 
 /////////// Sparse module ///////////
 
-    // dense = spasre * dense
+    // dense = sparse * dense
     template<typename Derived1, typename Derived2>
     Derived& lazyAssign(const SparseProduct<Derived1,Derived2,SparseTimeDenseProduct>& product);
-    // dense = dense * spasre
+    // dense = dense * sparse
     template<typename Derived1, typename Derived2>
     Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
 
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
-
+    
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN
     #endif
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 353aec48e..fe6c5027d 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -47,10 +47,10 @@ struct ei_product_packet_impl;
   * This class defines the typename Type representing the optimized product expression
   * between two matrix expressions. In practice, using ProductReturnType<Lhs,Rhs>::Type
   * is the recommended way to define the result type of a function returning an expression
-  * which involve a matrix product. The class Product or DiagonalProduct should never be
+  * which involve a matrix product. The class Product should never be
   * used directly.
   *
-  * \sa class Product, class DiagonalProduct, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
+  * \sa class Product, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
   */
 template<typename Lhs, typename Rhs, int ProductMode>
 struct ProductReturnType
@@ -62,7 +62,6 @@ struct ProductReturnType
 };
 
 // cache friendly specialization
-// note that there is a DiagonalProduct specialization in DiagonalProduct.h
 template<typename Lhs, typename Rhs>
 struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
 {
@@ -78,15 +77,12 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
 /*  Helper class to determine the type of the product, can be either:
  *    - NormalProduct
  *    - CacheFriendlyProduct
- *    - DiagonalProduct
  */
 template<typename Lhs, typename Rhs> struct ei_product_mode
 {
   enum{
 
-    value = ei_is_diagonal<Rhs>::ret || ei_is_diagonal<Lhs>::ret
-          ? DiagonalProduct
-          : Lhs::MaxColsAtCompileTime == Dynamic
+    value = Lhs::MaxColsAtCompileTime == Dynamic
             && ( Lhs::MaxRowsAtCompileTime == Dynamic
               || Rhs::MaxColsAtCompileTime == Dynamic )
             && (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit)  && (!(Lhs::Flags&DirectAccessBit))))
@@ -290,18 +286,6 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
 }
 
-/** replaces \c *this by \c *this * \a other.
-  *
-  * \returns a reference to \c *this
-  */
-template<typename Derived>
-template<typename OtherDerived>
-inline Derived &
-MatrixBase<Derived>::operator*=(const MatrixBase<OtherDerived> &other)
-{
-  return *this = *this * other;
-}
-
 /***************************************************************************
 * Normal product .coeff() implementation (with meta-unrolling)
 ***************************************************************************/
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index c2641537a..35e4767b8 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -205,12 +205,14 @@ template<typename T> struct ei_is_diagonal
   };
 };
 
+enum { DiagonalOnTheLeft, DiagonalOnTheRight };
+
 enum { Aligned, Unaligned };
 enum { ForceAligned, AsRequested };
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 enum DirectionType { Vertical, Horizontal, BothDirections };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
 
 enum {
   /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 9606771a1..3e4c1bf0b 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -28,6 +28,9 @@
 template<typename T> struct ei_traits;
 template<typename T> struct NumTraits;
 
+template<typename Derived> struct AnyMatrixBase;
+template<typename Derived> struct MultiplierBase;
+
 template<typename _Scalar, int _Rows, int _Cols,
          int _Options = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION | AutoAlign,
          int _MaxRows = _Rows, int _MaxCols = _Cols> class Matrix;
@@ -46,10 +49,13 @@ template<typename UnaryOp,   typename MatrixType>         class CwiseUnaryOp;
 template<typename ViewOp,    typename MatrixType>         class CwiseUnaryView;
 template<typename BinaryOp,  typename Lhs, typename Rhs>  class CwiseBinaryOp;
 template<typename Lhs, typename Rhs, int ProductMode> class Product;
-template<typename CoeffsVectorType, typename Derived> class DiagonalMatrixBase;
-template<typename CoeffsVectorType> class DiagonalMatrixWrapper;
+
+template<typename Derived> class DiagonalBase;
+template<typename _DiagonalVectorType> class DiagonalWrapper;
 template<typename _Scalar, int _Size> class DiagonalMatrix;
+template<typename MatrixType, typename DiagonalType, int Order> class DiagonalProduct;
 template<typename MatrixType, int Index> class Diagonal;
+
 template<typename MatrixType, int PacketAccess = AsRequested> class Map;
 template<typename MatrixType, unsigned int Mode> class Part;
 template<typename MatrixType, unsigned int Mode> class Extract;
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 8f62f9b69..38af7caf3 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -233,30 +233,16 @@ using Eigen::ei_cos;
 // needed to define it here as escaping characters in CMake add_definition's argument seems very problematic.
 #define EIGEN_DOCS_IO_FORMAT IOFormat(3, AlignCols, " ", "\n", "", "")
 
-#define EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \
-template<typename OtherDerived> \
-EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::MatrixBase<OtherDerived>& other) \
-{ \
-  return Base::operator Op(other.derived()); \
-} \
-EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \
-{ \
-  return Base::operator Op(other); \
-}
-
-#define EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \
-template<typename Other> \
-EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \
-{ \
-  return Base::operator Op(scalar); \
-}
-
 #define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
-EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \
-EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \
-EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \
-EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
-EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
+using Base::operator =; \
+using Base::operator +=; \
+using Base::operator -=; \
+using Base::operator *=; \
+using Base::operator /=; \
+EIGEN_STRONG_INLINE Derived& operator=(const Derived& other) \
+{ \
+  return Base::operator=(other); \
+}
 
 #define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \
 typedef BaseClass Base; \
diff --git a/Eigen/src/Geometry/RotationBase.h b/Eigen/src/Geometry/RotationBase.h
index 632ea3991..98d68c294 100644
--- a/Eigen/src/Geometry/RotationBase.h
+++ b/Eigen/src/Geometry/RotationBase.h
@@ -91,12 +91,12 @@ class RotationBase
       */
     template<typename OtherDerived>
     inline typename generic_product_selector<OtherDerived,OtherDerived::IsVectorAtCompileTime>::ReturnType
-    operator*(const MatrixBase<OtherDerived>& e) const
+    operator*(const MultiplierBase<OtherDerived>& e) const
     { return generic_product_selector<OtherDerived>::run(derived(), e.derived()); }
 
     /** \returns the concatenation of a linear transformation \a l with the rotation \a r */
     template<typename OtherDerived> friend
-    inline RotationMatrixType operator*(const MatrixBase<OtherDerived>& l, const Derived& r)
+    inline RotationMatrixType operator*(const MultiplierBase<OtherDerived>& l, const Derived& r)
     { return l.derived() * r.toRotationMatrix(); }
 
     /** \returns the concatenation of the rotation \c *this with a transformation \a t */
diff --git a/Eigen/src/Geometry/Scaling.h b/Eigen/src/Geometry/Scaling.h
index 785a5f40c..ce191b5da 100644
--- a/Eigen/src/Geometry/Scaling.h
+++ b/Eigen/src/Geometry/Scaling.h
@@ -141,7 +141,7 @@ static inline DiagonalMatrix<Scalar,3> Scaling(Scalar sx, Scalar sy, Scalar sz)
   * This is an alias for coeffs.asDiagonal()
   */
 template<typename Derived>
-static inline const DiagonalMatrixWrapper<Derived> Scaling(const MatrixBase<Derived>& coeffs)
+static inline const DiagonalWrapper<Derived> Scaling(const MatrixBase<Derived>& coeffs)
 { return coeffs.asDiagonal(); }
 
 /** \addtogroup Geometry_Module */
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index c186e5ec4..eb6048c47 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -226,14 +226,14 @@ public:
 
   /** Constructs and initializes a transformation from a Dim^2 or a (Dim+1)^2 matrix. */
   template<typename OtherDerived>
-  inline explicit Transform(const MatrixBase<OtherDerived>& other)
+  inline explicit Transform(const AnyMatrixBase<OtherDerived>& other)
   {
     ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived());
   }
 
   /** Set \c *this from a Dim^2 or (Dim+1)^2 matrix. */
   template<typename OtherDerived>
-  inline Transform& operator=(const MatrixBase<OtherDerived>& other)
+  inline Transform& operator=(const AnyMatrixBase<OtherDerived>& other)
   {
     ei_transform_construct_from_matrix<OtherDerived,Mode,Dim,HDim>::run(this, other.derived());
     return *this;
@@ -310,7 +310,7 @@ public:
   // note: this function is defined here because some compilers cannot find the respective declaration
   template<typename OtherDerived>
   inline const typename ei_transform_right_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType
-  operator * (const MatrixBase<OtherDerived> &other) const
+  operator * (const MultiplierBase<OtherDerived> &other) const
   { return ei_transform_right_product_impl<OtherDerived,Mode,Dim,HDim>::run(*this,other.derived()); }
 
   /** \returns the product expression of a transformation matrix \a a times a transform \a b
@@ -322,11 +322,11 @@ public:
     */
   template<typename OtherDerived> friend
   inline const typename ei_transform_left_product_impl<OtherDerived,Mode,_Dim,_Dim+1>::ResultType
-  operator * (const MatrixBase<OtherDerived> &a, const Transform &b)
+  operator * (const MultiplierBase<OtherDerived> &a, const Transform &b)
   { return ei_transform_left_product_impl<OtherDerived,Mode,Dim,HDim>::run(a.derived(),b); }
 
   template<typename OtherDerived>
-  inline Transform& operator*=(const MatrixBase<OtherDerived>& other) { return *this = *this * other; }
+  inline Transform& operator*=(const MultiplierBase<OtherDerived>& other) { return *this = *this * other; }
 
   /** Contatenates two transformations */
   inline const Transform operator * (const Transform& other) const
@@ -944,7 +944,7 @@ struct ei_transform_take_affine_part<Transform<Scalar,Dim,AffineCompact> > {
 template<typename Other, int Mode, int Dim, int HDim>
 struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,Dim>
 {
-  static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
+  static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
   {
     transform->linear() = other;
     transform->translation().setZero();
@@ -955,7 +955,7 @@ struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,Dim>
 template<typename Other, int Mode, int Dim, int HDim>
 struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,HDim>
 {
-  static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
+  static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
   {
     transform->affine() = other;
     transform->makeAffine();
@@ -965,20 +965,32 @@ struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, Dim,HDim>
 template<typename Other, int Mode, int Dim, int HDim>
 struct ei_transform_construct_from_matrix<Other, Mode,Dim,HDim, HDim,HDim>
 {
-  static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,Mode> *transform, const Other& other)
+  static inline void run(Transform<typename Other::Scalar,Dim,Mode> *transform, const Other& other)
   { transform->matrix() = other; }
 };
 
 template<typename Other, int Dim, int HDim>
 struct ei_transform_construct_from_matrix<Other, AffineCompact,Dim,HDim, HDim,HDim>
 {
-  static inline void run(Transform<typename ei_traits<Other>::Scalar,Dim,AffineCompact> *transform, const Other& other)
+  static inline void run(Transform<typename Other::Scalar,Dim,AffineCompact> *transform, const Other& other)
   { transform->matrix() = other.template block<Dim,HDim>(0,0); }
 };
 
-/*****************************************************
-*** Specializations of operator* with a MatrixBase ***
-*****************************************************/
+/*********************************************************
+*** Specializations of operator* with a MultiplierBase ***
+*********************************************************/
+
+// ei_general_product_return_type is a generalization of ProductReturnType, for all types (including e.g. DiagonalBase...),
+// instead of being restricted to MatrixBase.
+template<typename Lhs, typename Rhs> struct ei_general_product_return_type;
+template<typename D1, typename D2> struct ei_general_product_return_type<MatrixBase<D1>, MatrixBase<D2> >
+ : ProductReturnType<D1,D2> {};
+template<typename Lhs, typename D2> struct ei_general_product_return_type<Lhs, MatrixBase<D2> >
+{ typedef D2 Type; };
+template<typename D1, typename Rhs> struct ei_general_product_return_type<MatrixBase<D1>, Rhs >
+{ typedef D1 Type; };
+ 
+
 
 // Projective * set of homogeneous column vectors
 template<typename Other, int Dim, int HDim>
diff --git a/Eigen/src/Geometry/Translation.h b/Eigen/src/Geometry/Translation.h
index ba54c4631..a90e1b2f2 100644
--- a/Eigen/src/Geometry/Translation.h
+++ b/Eigen/src/Geometry/Translation.h
@@ -93,7 +93,7 @@ public:
 
   /** Concatenates a translation and a linear transformation */
   template<typename OtherDerived>
-  inline AffineTransformType operator* (const MatrixBase<OtherDerived>& linear) const;
+  inline AffineTransformType operator* (const MultiplierBase<OtherDerived>& linear) const;
 
   /** Concatenates a translation and a rotation */
   template<typename Derived>
@@ -103,7 +103,7 @@ public:
   /** \returns the concatenation of a linear transformation \a l with the translation \a t */
   // its a nightmare to define a templated friend function outside its declaration
   template<typename OtherDerived> friend
-  inline AffineTransformType operator*(const MatrixBase<OtherDerived>& linear, const Translation& t)
+  inline AffineTransformType operator*(const MultiplierBase<OtherDerived>& linear, const Translation& t)
   {
     AffineTransformType res;
     res.matrix().setZero();
@@ -182,7 +182,7 @@ Translation<Scalar,Dim>::operator* (const UniformScaling<Scalar>& other) const
 template<typename Scalar, int Dim>
 template<typename OtherDerived>
 inline typename Translation<Scalar,Dim>::AffineTransformType
-Translation<Scalar,Dim>::operator* (const MatrixBase<OtherDerived>& linear) const
+Translation<Scalar,Dim>::operator* (const MultiplierBase<OtherDerived>& linear) const
 {
   AffineTransformType res;
   res.matrix().setZero();
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a266ec482..68c42579e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -98,6 +98,7 @@ ei_add_test(redux)
 ei_add_test(product_small)
 ei_add_test(product_large ${EI_OFLAG})
 ei_add_test(product_selfadjoint)
+ei_add_test(diagonalmatrices)
 ei_add_test(adjoint)
 ei_add_test(submatrices)
 ei_add_test(miscmatrices)
diff --git a/test/diagonalmatrices.cpp b/test/diagonalmatrices.cpp
new file mode 100644
index 000000000..0a8b7086b
--- /dev/null
+++ b/test/diagonalmatrices.cpp
@@ -0,0 +1,95 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// 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 "main.h"
+
+template<typename MatrixType> void diagonalmatrices(const MatrixType& m)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
+  typedef Matrix<Scalar, Rows, 1> VectorType;
+  typedef Matrix<Scalar, 1, Cols> RowVectorType;
+  typedef Matrix<Scalar, Rows, Rows> SquareMatrixType;
+  typedef DiagonalMatrix<Scalar, Rows> LeftDiagonalMatrix;
+  typedef DiagonalMatrix<Scalar, Cols> RightDiagonalMatrix;
+  
+  int rows = m.rows();
+  int cols = m.cols();
+
+  MatrixType m1 = MatrixType::Random(rows, cols),
+             m2 = MatrixType::Random(rows, cols);
+  VectorType v1 = VectorType::Random(rows),
+             v2 = VectorType::Random(rows);
+  RowVectorType rv1 = RowVectorType::Random(cols),
+             rv2 = RowVectorType::Random(cols);
+  LeftDiagonalMatrix ldm1(v1), ldm2(v2);
+  RightDiagonalMatrix rdm1(rv1), rdm2(rv2);
+  
+  int i = ei_random<int>(0, rows-1);
+  int j = ei_random<int>(0, cols-1);
+  
+  VERIFY_IS_APPROX( ((ldm1 * m1)(i,j))  , ldm1.diagonal()(i) * m1(i,j) );
+  VERIFY_IS_APPROX( ((ldm1 * (m1+m2))(i,j))  , ldm1.diagonal()(i) * (m1+m2)(i,j) );
+  VERIFY_IS_APPROX( ((m1 * rdm1)(i,j))  , rdm1.diagonal()(j) * m1(i,j) );
+  VERIFY_IS_APPROX( ((v1.asDiagonal() * m1)(i,j))  , v1(i) * m1(i,j) );
+  VERIFY_IS_APPROX( ((m1 * rv1.asDiagonal())(i,j))  , rv1(j) * m1(i,j) );
+  VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * m1)(i,j))  , (v1+v2)(i) * m1(i,j) );
+  VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * (m1+m2))(i,j))  , (v1+v2)(i) * (m1+m2)(i,j) );
+  VERIFY_IS_APPROX( ((m1 * (rv1+rv2).asDiagonal())(i,j))  , (rv1+rv2)(j) * m1(i,j) );
+  VERIFY_IS_APPROX( (((m1+m2) * (rv1+rv2).asDiagonal())(i,j))  , (rv1+rv2)(j) * (m1+m2)(i,j) );
+  
+  SquareMatrixType sq_m1 (v1.asDiagonal());
+  VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
+  sq_m1 = v1.asDiagonal();
+  VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
+  SquareMatrixType sq_m2 = v1.asDiagonal();
+  VERIFY_IS_APPROX(sq_m1, sq_m2);
+  
+  ldm1 = v1.asDiagonal();
+  LeftDiagonalMatrix ldm3(v1);
+  VERIFY_IS_APPROX(ldm1.diagonal(), ldm3.diagonal());
+  LeftDiagonalMatrix ldm4 = v1.asDiagonal();
+  VERIFY_IS_APPROX(ldm1.diagonal(), ldm4.diagonal());
+  
+  sq_m1.block(0,0,rows,rows) = ldm1;
+  VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
+  sq_m1.transpose() = ldm1;
+  VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
+}
+
+void test_diagonalmatrices()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( diagonalmatrices(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( diagonalmatrices(Matrix3f()) );
+    CALL_SUBTEST( diagonalmatrices(Matrix<double,3,3,RowMajor>()) );
+    CALL_SUBTEST( diagonalmatrices(Matrix4d()) );
+    CALL_SUBTEST( diagonalmatrices(Matrix<float,4,4,RowMajor>()) );
+    CALL_SUBTEST( diagonalmatrices(MatrixXcf(3, 5)) );
+    CALL_SUBTEST( diagonalmatrices(MatrixXi(10, 8)) );
+    CALL_SUBTEST( diagonalmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
+    CALL_SUBTEST( diagonalmatrices(MatrixXf(21, 24)) );
+  }
+}
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index b7202ab34..2be49faf4 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -57,7 +57,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   EigenSolver<MatrixType> ei1(a);
   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
-                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval());
+                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
 
 }
 
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index b8e2be98a..c93953714 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -75,7 +75,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     convert(gEvec, _evec);
 
     // test gsl itself !
-    VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal().eval(), largerEps));
+    VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
 
     // compare with eigen
     VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
@@ -86,7 +86,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     convert(gEval, _eval);
     convert(gEvec, _evec);
     // test GSL itself:
-    VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal().eval()), largerEps));
+    VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
 
     // compare with eigen
 //     std::cerr << _eval.transpose() << "\n" << eiSymmGen.eigenvalues().transpose() << "\n\n";
@@ -102,11 +102,11 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   #endif
 
   VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
-          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval(), largerEps));
+          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
 
   // generalized eigen problem Ax = lBx
   VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
-          symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
+          symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
   VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp
index fcd247bc5..a4f75e384 100644
--- a/test/geo_transformations.cpp
+++ b/test/geo_transformations.cpp
@@ -332,12 +332,6 @@ template<typename Scalar, int Mode> void transformations(void)
   Translation<double,3> tr1d = tr1.template cast<double>();
   VERIFY_IS_APPROX(tr1d.template cast<Scalar>(),tr1);
 
-  AlignedScaling3 sc1(v0);
-  DiagonalMatrix<float,3> sc1f; sc1f = sc1.template cast<float>();
-  VERIFY_IS_APPROX(sc1f.template cast<Scalar>(),sc1);
-  DiagonalMatrix<double,3> sc1d; sc1d = (sc1.template cast<double>());
-  VERIFY_IS_APPROX(sc1d.template cast<Scalar>(),sc1);
-
   AngleAxis<float> aa1f = aa1.template cast<float>();
   VERIFY_IS_APPROX(aa1f.template cast<Scalar>(),aa1);
   AngleAxis<double> aa1d = aa1.template cast<double>();
diff --git a/test/miscmatrices.cpp b/test/miscmatrices.cpp
index 5b0367be6..bf885e252 100644
--- a/test/miscmatrices.cpp
+++ b/test/miscmatrices.cpp
@@ -43,7 +43,7 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m)
   VectorType v1 = VectorType::Random(rows);
   v1[0];
   Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-  square = v1.asDiagonal();
+  square(v1.asDiagonal());
   if(r==r2) VERIFY_IS_APPROX(square(r,r2), v1[r]);
   else VERIFY_IS_MUCH_SMALLER_THAN(square(r,r2), static_cast<Scalar>(1));
   square = MatrixType::Zero(rows, rows);