* add a HouseholderSequence class (not good enough yet for Triadiagonalization and HessenbergDecomposition)

* rework a bit AnyMatrixBase, and mobe it to a separate file

5 files changed, 169 insertions, 69 deletions

diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/AnyMatrixBase.h
new file mode 100644
index 000000000..cd354d7b1
--- /dev/null
+++ b/Eigen/src/Core/AnyMatrixBase.h
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_ANYMATRIXBASE_H
+#define EIGEN_ANYMATRIXBASE_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
+{
+  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); }
+
+  /** \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(); }
+
+  /** \internal Don't use it, but do the equivalent: \code dst = *this; \endcode */
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  { derived().evalTo(dst); }
+
+  /** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */
+  template<typename Dest> inline void addToDense(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalTo(res);
+    dst += res;
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */
+  template<typename Dest> inline void subToDense(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalTo(res);
+    dst -= res;
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheRight(*this); \endcode */
+  template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    dst = dst * this->derived();
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheLeft(*this); \endcode */
+  template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    dst = this->derived() * dst;
+  }
+
+};
+
+/***************************************************************************
+* Implementation of matrix base methods
+***************************************************************************/
+
+/** Copies the generic expression \a other into *this. \returns a reference to *this.
+  * The expression must provide a (templated) evalToDense(Derived& dst) const function
+  * which does the actual job. In practice, this allows any user to write its own
+  * special matrix without having to modify MatrixBase */
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().evalTo(derived());
+  return derived();
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().addToDense(derived());
+  return derived();
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().subToDense(derived());
+  return 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 AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheRight(derived());
+  return derived();
+}
+
+/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */
+template<typename Derived>
+template<typename OtherDerived>
+inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheRight(derived());
+}
+
+/** replaces \c *this by \c *this * \a other. */
+template<typename Derived>
+template<typename OtherDerived>
+inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheLeft(derived());
+}
+
+#endif // EIGEN_ANYMATRIXBASE_H
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f267aa34d..4835f167c 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -26,44 +26,6 @@
 #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
-{
-  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); }
-  /** \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
   *
   * \brief Base class for all matrices, vectors, and expressions
@@ -96,7 +58,6 @@ template<typename Derived> class MatrixBase
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 {
   public:
-
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
                 typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*;
@@ -301,21 +262,14 @@ template<typename Derived> class MatrixBase
       */
     Derived& operator=(const MatrixBase& other);
 
-    /** Copies the generic expression \a other into *this. \returns a reference to *this.
-      * The expression must provide a (templated) evalToDense(Derived& dst) const function
-      * which does the actual job. In practice, this allows any user to write its own
-      * special matrix without having to modify MatrixBase */
     template<typename OtherDerived>
-    Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().evalToDense(derived()); return derived(); }
+    Derived& operator=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived>
-    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().addToDense(derived()); return derived(); }
+    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived>
-    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().subToDense(derived()); return derived(); }
+    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
@@ -436,6 +390,12 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
 
+    template<typename OtherDerived>
+    void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other);
+
+    template<typename OtherDerived>
+    void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other);
+
     template<typename DiagonalDerived>
     const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
     operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index e7227d4f6..7f0c2df6e 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -434,18 +434,4 @@ 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 AnyMatrixBase<OtherDerived> &other)
-{
-  return derived() = derived() * other.derived();
-}
-
 #endif // EIGEN_PRODUCT_H
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index b0362f20c..17726bca3 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -91,9 +91,9 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
     #endif // not EIGEN_PARSED_BY_DOXYGEN
 
     template<typename DenseDerived>
-    void evalToDense(MatrixBase<DenseDerived> &other) const;
+    void evalTo(MatrixBase<DenseDerived> &other) const;
     template<typename DenseDerived>
-    void evalToDenseLazy(MatrixBase<DenseDerived> &other) const;
+    void evalToLazy(MatrixBase<DenseDerived> &other) const;
 
   protected:
 
@@ -546,23 +546,23 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
 template<typename DenseDerived>
-void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
+void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
 {
   if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit)
   {
     typename Derived::PlainMatrixType other_evaluated(rows(), cols());
-    evalToDenseLazy(other_evaluated);
+    evalToLazy(other_evaluated);
     other.derived().swap(other_evaluated);
   }
   else
-    evalToDenseLazy(other.derived());
+    evalToLazy(other.derived());
 }
 
 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
 template<typename DenseDerived>
-void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const
+void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
 {
   const bool unroll =   DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2
                      <= EIGEN_UNROLLING_LIMIT;
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index c5f27d80b..3f66738f0 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -123,6 +123,7 @@ template<typename MatrixType> class SVD;
 template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
 template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
 template<typename MatrixType> class LDLT;
+template<typename VectorsType, typename CoeffsType> class HouseholderSequence;
 template<typename Scalar>     class PlanarRotation;
 
 // Geometry module: