* Bye bye MultiplierBase, extend a bit AnyMatrixBase to allow =, +=, and -=
* This probably makes ReturnByValue needless

Eigen/src/Core/BandMatrix.h

diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h
index 2da463a..c226969 100644
--- a/Eigen/src/Core/BandMatrix.h
+++ b/Eigen/src/Core/BandMatrix.h
@@ -57,7 +57,7 @@
 };
 
 template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
-class BandMatrix : public MultiplierBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
+class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
 {
   public:
 

Eigen/src/Core/DiagonalMatrix.h

diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 5fc80c9..ebbed15 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -27,7 +27,7 @@
 #define EIGEN_DIAGONALMATRIX_H
 
 template<typename Derived>
-class DiagonalBase : public MultiplierBase<Derived>
+class DiagonalBase : public AnyMatrixBase<Derived>
 {
   public:
     typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType;
@@ -52,6 +52,12 @@
     DenseMatrixType toDenseMatrix() const { return derived(); }
     template<typename DenseDerived>
     void evalToDense(MatrixBase<DenseDerived> &other) const;
+    template<typename DenseDerived>
+    void addToDense(MatrixBase<DenseDerived> &other) const
+    { other.diagonal() += diagonal(); }
+    template<typename DenseDerived>
+    void subToDense(MatrixBase<DenseDerived> &other) const
+    { other.diagonal() -= diagonal(); }
 
     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
@@ -84,6 +90,7 @@
   */
 template<typename _Scalar, int _Size>
 struct ei_traits<DiagonalMatrix<_Scalar,_Size> >
+ : ei_traits<Matrix<_Scalar,_Size,_Size> >
 {
   typedef Matrix<_Scalar,_Size,1> DiagonalVectorType;
 };
@@ -170,6 +177,14 @@
 struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
 {
   typedef _DiagonalVectorType DiagonalVectorType;
+  typedef typename DiagonalVectorType::Scalar Scalar;
+  enum {
+    RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    MaxRowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    MaxColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
+    Flags = 0
+  };
 };
 
 template<typename _DiagonalVectorType>

Eigen/src/Core/Matrix.h

diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 8937596..c31acab 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -462,6 +462,7 @@
       : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols())
     {
       _check_template_params();
+      resize(other.rows(), other.cols());
       *this = other;
     }
 

Eigen/src/Core/MatrixBase.h

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index b881c09..30cfbb1 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -37,19 +37,31 @@
   */
 template<typename Derived> struct AnyMatrixBase
 {
+  typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
+  
   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;
+  /** \returns the number of rows. \sa cols(), RowsAtCompileTime */
+  inline int rows() const { return derived().rows(); }
+  /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
+  inline int cols() const { return derived().cols(); }
+
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  { derived().evalTo(dst); }
+
+  template<typename Dest> inline void addToDense(Dest& dst) const
+  {
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalToDense(res);
+    dst += res;
+  }
+  
+  template<typename Dest> inline void subToDense(Dest& dst) const
+  {
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalToDense(res);
+    dst -= res;
+  }
 };
 
 /** \class MatrixBase
@@ -79,7 +91,7 @@
   */
 template<typename Derived> class MatrixBase
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-  : public MultiplierBase<Derived>,
+  : public AnyMatrixBase<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
@@ -298,12 +310,16 @@
     Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
     { other.derived().evalToDense(derived()); return derived(); }
 
+    template<typename OtherDerived>
+    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
+    { other.derived().addToDense(derived()); return derived(); }
+
+    template<typename OtherDerived>
+    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
+    { other.derived().subToDense(derived()); return derived(); }
+
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
-    template<typename OtherDerived,typename OtherEvalType>
-    Derived& operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
-    template<typename OtherDerived,typename OtherEvalType>
-    Derived& operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** Copies \a other into *this without evaluating other. \returns a reference to *this. */
@@ -410,7 +426,7 @@
     operator*(const MatrixBase<OtherDerived> &other) const;
 
     template<typename OtherDerived>
-    Derived& operator*=(const MultiplierBase<OtherDerived>& other);
+    Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
 
     template<typename DiagonalDerived>
     const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
@@ -645,7 +661,7 @@
     void visit(Visitor& func) const;
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-    using MultiplierBase<Derived>::derived;
+    using AnyMatrixBase<Derived>::derived;
     inline Derived& const_cast_derived() const
     { return *static_cast<Derived*>(const_cast<MatrixBase*>(this)); }
 #endif // not EIGEN_PARSED_BY_DOXYGEN

Eigen/src/Core/Product.h

diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index ff45cba..1a32eb5 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -294,7 +294,7 @@
 template<typename Derived>
 template<typename OtherDerived>
 inline Derived &
-MatrixBase<Derived>::operator*=(const MultiplierBase<OtherDerived> &other)
+MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
 {
   return derived() = derived() * other.derived();
 }

Eigen/src/Core/ReturnByValue.h

diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h
index 58b205e..3f2b478 100644
--- a/Eigen/src/Core/ReturnByValue.h
+++ b/Eigen/src/Core/ReturnByValue.h
@@ -48,14 +48,6 @@
   public:
     template<typename Dest> inline void evalTo(Dest& dst) const
     { static_cast<const Functor*>(this)->evalTo(dst); }
-    template<typename Dest> inline void addTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_addTo(dst); }
-    template<typename Dest> inline void subTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_subTo(dst); }
-    template<typename Dest> inline void _addTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst += res; }
-    template<typename Dest> inline void _subTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst -= res; }
 };
 
 template<typename Functor, typename _Scalar,int _Rows,int _Cols,int _Options,int _MaxRows,int _MaxCols>
@@ -68,14 +60,6 @@
     template<typename Dest>
     inline void evalTo(Dest& dst) const
     { static_cast<const Functor* const>(this)->evalTo(dst); }
-    template<typename Dest> inline void addTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_addTo(dst); }
-    template<typename Dest> inline void subTo(Dest& dst) const
-    { static_cast<const Functor*>(this)->_subTo(dst); }
-    template<typename Dest> inline void _addTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst += res; }
-    template<typename Dest> inline void _subTo(Dest& dst) const
-    { EvalType res; evalTo(res); dst -= res; }
     inline int rows() const { return static_cast<const Functor* const>(this)->rows(); }
     inline int cols() const { return static_cast<const Functor* const>(this)->cols(); }
 };
@@ -88,20 +72,4 @@
   return derived();
 }
 
-template<typename Derived>
-template<typename OtherDerived,typename OtherEvalType>
-Derived& MatrixBase<Derived>::operator+=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
-{
-  other.addTo(derived());
-  return derived();
-}
-
-template<typename Derived>
-template<typename OtherDerived,typename OtherEvalType>
-Derived& MatrixBase<Derived>::operator-=(const ReturnByValue<OtherDerived,OtherEvalType>& other)
-{
-  other.subTo(derived());
-  return derived();
-}
-
 #endif // EIGEN_RETURNBYVALUE_H

Eigen/src/Core/SelfAdjointView.h

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index c21f3a3..0a4ba17 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -55,7 +55,7 @@
 
 template <typename Lhs, int LhsMode, bool LhsIsVector,
           typename Rhs, int RhsMode, bool RhsIsVector>
-struct ei_selfadjoint_product_returntype;
+struct SelfadjointProductMatrix;
 
 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
@@ -100,20 +100,20 @@
 
     /** Efficient self-adjoint matrix times vector/matrix product */
     template<typename OtherDerived>
-    ei_selfadjoint_product_returntype<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
+    SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
     operator*(const MatrixBase<OtherDerived>& rhs) const
     {
-      return ei_selfadjoint_product_returntype
+      return SelfadjointProductMatrix
               <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
               (m_matrix, rhs.derived());
     }
 
     /** Efficient vector/matrix times self-adjoint matrix product */
     template<typename OtherDerived> friend
-    ei_selfadjoint_product_returntype<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
+    SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
     operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
     {
-      return ei_selfadjoint_product_returntype
+      return SelfadjointProductMatrix
               <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
               (lhs.derived(),rhs.m_matrix);
     }
@@ -201,10 +201,13 @@
 ***************************************************************************/
 
 template<typename Lhs, int LhsMode, typename Rhs>
-struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>
-  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,0,true>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<typename Lhs, int LhsMode, typename Rhs>
+struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
+  : public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> >
 {
   typedef typename Lhs::Scalar Scalar;
 
@@ -224,19 +227,19 @@
     LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit)
   };
 
-  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
   inline int rows() const { return m_lhs.rows(); }
   inline int cols() const { return m_rhs.cols(); }
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);
@@ -272,12 +275,15 @@
 ***************************************************************************/
 
 template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
-struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>
-  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
+struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>
+  : public AnyMatrixBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> >
 {
-  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  SelfadjointProductMatrix(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
@@ -305,12 +311,12 @@
     RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit
   };
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);

Eigen/src/Core/TriangularMatrix.h

diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 8b6c9a2..c262ea7 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -43,7 +43,7 @@
   *
   * \sa MatrixBase::part()
   */
-template<typename Derived> class TriangularBase : public MultiplierBase<Derived>
+template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
 {
   public:
 
@@ -145,7 +145,7 @@
 template<int Mode, bool LhsIsTriangular,
          typename Lhs, bool LhsIsVector,
          typename Rhs, bool RhsIsVector>
-struct ei_triangular_product_returntype;
+struct TriangularProduct;
 
 template<typename _MatrixType, unsigned int _Mode> class TriangularView
   : public TriangularBase<TriangularView<_MatrixType, _Mode> >
@@ -253,20 +253,20 @@
 
     /** Efficient triangular matrix times vector/matrix product */
     template<typename OtherDerived>
-    ei_triangular_product_returntype<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
+    TriangularProduct<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
     operator*(const MatrixBase<OtherDerived>& rhs) const
     {
-      return ei_triangular_product_returntype
+      return TriangularProduct
               <Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
               (m_matrix, rhs.derived());
     }
 
     /** Efficient vector/matrix times triangular matrix product */
     template<typename OtherDerived> friend
-    ei_triangular_product_returntype<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
+    TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
     {
-      return ei_triangular_product_returntype
+      return TriangularProduct
               <Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
               (lhs.derived(),rhs.m_matrix);
     }

Eigen/src/Core/products/TriangularMatrixMatrix.h

diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index ce18941..f69c043 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -321,12 +321,15 @@
 ***************************************************************************/
 
 template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
-struct ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
-  : public ReturnByValue<ei_triangular_product_returntype<Mode,LhsIsTriangular,Lhs,false,Rhs,false>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
+struct TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
+  : public AnyMatrixBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,false> >
 {
-  ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  TriangularProduct(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
@@ -347,12 +350,12 @@
   typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
   typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.resize(m_lhs.rows(), m_rhs.cols());
     dst.setZero();

Eigen/src/Core/products/TriangularMatrixVector.h

diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h
index 18d76b9..42239fa 100644
--- a/Eigen/src/Core/products/TriangularMatrixVector.h
+++ b/Eigen/src/Core/products/TriangularMatrixVector.h
@@ -118,10 +118,13 @@
 ***************************************************************************/
 
 template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
-struct ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>
-  : public ReturnByValue<ei_triangular_product_returntype<Mode,true,Lhs,false,Rhs,true>,
-                         Matrix<typename ei_traits<Rhs>::Scalar,
-                                Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+struct ei_traits<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
+ : ei_traits<Matrix<typename ei_traits<Rhs>::Scalar,Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+{};
+
+template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs>
+struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
+  : public AnyMatrixBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true> >
 {
   typedef typename Lhs::Scalar Scalar;
 
@@ -137,19 +140,19 @@
   typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
   typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
 
-  ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  TriangularProduct(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
   inline int rows() const { return m_lhs.rows(); }
   inline int cols() const { return m_rhs.cols(); }
 
-  template<typename Dest> inline void _addTo(Dest& dst) const
+  template<typename Dest> inline void addToDense(Dest& dst) const
   { evalTo(dst,1); }
-  template<typename Dest> inline void _subTo(Dest& dst) const
+  template<typename Dest> inline void subToDense(Dest& dst) const
   { evalTo(dst,-1); }
 
-  template<typename Dest> void evalTo(Dest& dst) const
+  template<typename Dest> void evalToDense(Dest& dst) const
   {
     dst.setZero();
     evalTo(dst,1);

Eigen/src/Core/util/ForwardDeclarations.h

diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index b4fbae2..310d0fb 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -29,7 +29,6 @@
 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,

bench/btl/libs/C_BLAS/C_BLAS_interface.hh

diff --git a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
index e08c4ae..c8b6988 100644
--- a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
+++ b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
@@ -237,6 +237,15 @@
     spotrf_(&uplo, &N, C, &N, &info);
   }
 
+  static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
+    int N2 = N*N;
+    scopy_(&N2, X, &intone, C, &intone);
+    char uplo = 'L';
+    int info = 0;
+    int * ipiv = (int*)alloca(sizeof(int)*N);
+    sgetrf_(&N, &N, C, &N, ipiv, &info);
+  }
+
   #ifdef HAS_LAPACK
 
   static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
@@ -249,14 +258,7 @@
     sgetc2_(&N, C, &N, ipiv, jpiv, &info);
   }
 
-  static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
-    int N2 = N*N;
-    scopy_(&N2, X, &intone, C, &intone);
-    char uplo = 'L';
-    int info = 0;
-    int * ipiv = (int*)alloca(sizeof(int)*N);
-    sgetrf_(&N, &N, C, &N, ipiv, &info);
-  }
+  
 
   static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
 #ifdef PUREBLAS

bench/btl/libs/C_BLAS/main.cpp

diff --git a/bench/btl/libs/C_BLAS/main.cpp b/bench/btl/libs/C_BLAS/main.cpp
index 5c2f8a3..bb3022d 100644
--- a/bench/btl/libs/C_BLAS/main.cpp
+++ b/bench/btl/libs/C_BLAS/main.cpp
@@ -24,7 +24,7 @@
 
 #include "action_cholesky.hh"
 #include "action_lu_decomp.hh"
-// #include "action_partial_lu_decomp.hh"
+#include "action_partial_lu.hh"
 #include "action_trisolve_matrix.hh"
 
 #ifdef HAS_LAPACK
@@ -52,10 +52,10 @@
   bench<Action_trisolve_matrix<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   bench<Action_cholesky<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
+  bench<Action_partial_lu<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   #ifdef HAS_LAPACK
   bench<Action_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
-//   bench<Action_partial_lu_decomp<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_hessenberg<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_tridiagonalization<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   #endif

bench/btl/libs/eigen2/eigen2_interface.hh

diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh
index ffc800a..7a01472 100644
--- a/bench/btl/libs/eigen2/eigen2_interface.hh
+++ b/bench/btl/libs/eigen2/eigen2_interface.hh
@@ -190,10 +190,10 @@
   }
 
   static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
-//     X = L.template triangularView<LowerTriangular>().solve(B);
-    X = B;
-    ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
-      ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
+    X = L.template triangularView<LowerTriangular>().solve(B);
+//     
+//     ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular>
+//       ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride());
   }
 
   static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){

bench/btl/libs/eigen2/main_adv.cpp

diff --git a/bench/btl/libs/eigen2/main_adv.cpp b/bench/btl/libs/eigen2/main_adv.cpp
index c3e5fac..d98c0cd 100644
--- a/bench/btl/libs/eigen2/main_adv.cpp
+++ b/bench/btl/libs/eigen2/main_adv.cpp
@@ -23,7 +23,7 @@
 #include "action_cholesky.hh"
 #include "action_hessenberg.hh"
 #include "action_lu_decomp.hh"
-// #include "action_partial_lu_decomp.hh"
+#include "action_partial_lu.hh"
 
 BTL_MAIN;
 
@@ -33,7 +33,7 @@
   bench<Action_trisolve_matrix<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_cholesky<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
-//   bench<Action_partial_lu_decomp<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
+  bench<Action_partial_lu<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
 
   bench<Action_hessenberg<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
   bench<Action_tridiagonalization<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);