merge

Eigen/Core

@@ -143,6 +143,7 @@ namespace Eigen {
 
 #include "src/Core/Functors.h"
 #include "src/Core/MatrixBase.h"
+#include "src/Core/AnyMatrixBase.h"
 #include "src/Core/Coeffs.h"
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN // work around Doxygen bug triggered by Assign.h r814874

Eigen/Householder

@@ -16,6 +16,7 @@ namespace Eigen {
   */
 
 #include "src/Householder/Householder.h"
+#include "src/Householder/HouseholderSequence.h"
 
 } // namespace Eigen
 

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

Eigen/src/Core/Matrix.h

@@ -25,6 +25,7 @@
 #ifndef EIGEN_MATRIX_H
 #define EIGEN_MATRIX_H
 
+template <typename Derived, typename OtherDerived, bool IsVector = static_cast<bool>(Derived::IsVectorAtCompileTime)> struct ei_conservative_resize_like_impl;
 
 /** \class Matrix
   *
@@ -324,40 +325,28 @@ class Matrix
 
     /** Resizes \c *this to a \a rows x \a cols matrix while leaving old values of *this untouched.
     *
-    * This method is intended for dynamic-size matrices, although it is legal to call it on any
+    * This method is intended for dynamic-size matrices. If you only want to change the number
-    * matrix as long as fixed dimensions are left unchanged. If you only want to change the number
     * of rows and/or of columns, you can use conservativeResize(NoChange_t, int),
     * conservativeResize(int, NoChange_t).
     *
     * The top-left part of the resized matrix will be the same as the overlapping top-left corner
-    * of *this. In case values need to be appended to the matrix they will be uninitialized per
+    * of *this. In case values need to be appended to the matrix they will be uninitialized.
-    * default and set to zero when init_with_zero is set to true.
     */
-    inline void conservativeResize(int rows, int cols, bool init_with_zero = false)
+    EIGEN_STRONG_INLINE void conservativeResize(int rows, int cols)
     {
-      // Note: Here is space for improvement. Basically, for conservativeResize(int,int),
+      conservativeResizeLike(PlainMatrixType(rows, cols));
-      // neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the
-      // dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or
-      // conservativeResize(NoChange_t, int cols). For these methods new static asserts like
-      // EIGEN_STATIC_ASSERT_DYNAMIC_ROWS and EIGEN_STATIC_ASSERT_DYNAMIC_COLS would be good.
-      EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Matrix)
-      PlainMatrixType tmp = init_with_zero ? PlainMatrixType::Zero(rows, cols) : PlainMatrixType(rows,cols);
-      const int common_rows = std::min(rows, this->rows());
-      const int common_cols = std::min(cols, this->cols());
-      tmp.block(0,0,common_rows,common_cols) = this->block(0,0,common_rows,common_cols);
-      this->derived().swap(tmp);
     }
 
-    EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t, bool init_with_zero = false)
+    EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t)
     {
-      // Note: see the comment in conservativeResize(int,int,bool)
+      // Note: see the comment in conservativeResize(int,int)
-      conservativeResize(rows, cols(), init_with_zero);
+      conservativeResize(rows, cols());
     }
 
-    EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols, bool init_with_zero = false)
+    EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols)
     {
-      // Note: see the comment in conservativeResize(int,int,bool)
+      // Note: see the comment in conservativeResize(int,int)
-      conservativeResize(rows(), cols, init_with_zero);
+      conservativeResize(rows(), cols);
     }
 
     /** Resizes \c *this to a vector of length \a size while retaining old values of *this.
@@ -366,21 +355,17 @@ class Matrix
     * partially dynamic matrices when the static dimension is anything other
     * than 1. For example it will not work with Matrix<double, 2, Dynamic>.
     *
-    * When values are appended, they will be uninitialized per default and set
+    * When values are appended, they will be uninitialized.
-    * to zero when init_with_zero is set to true.
     */
-    inline void conservativeResize(int size, bool init_with_zero = false)
+    EIGEN_STRONG_INLINE void conservativeResize(int size)
     {
-      EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
+      conservativeResizeLike(PlainMatrixType(size));
-      EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Matrix)
-
-      if (RowsAtCompileTime == 1 || ColsAtCompileTime == 1)
-      {
-        PlainMatrixType tmp = init_with_zero ? PlainMatrixType::Zero(size) : PlainMatrixType(size);
-        const int common_size = std::min<int>(this->size(),size);
-        tmp.segment(0,common_size) = this->segment(0,common_size);
-        this->derived().swap(tmp);
     }
+
+    template<typename OtherDerived>
+    EIGEN_STRONG_INLINE void conservativeResizeLike(const MatrixBase<OtherDerived>& other)
+    {
+      ei_conservative_resize_like_impl<Matrix, OtherDerived>::run(*this, other);
     }
 
     /** Copies the value of the expression \a other into \c *this with automatic resizing.
@@ -713,13 +698,45 @@ class Matrix
       m_storage.data()[1] = y;
     }
 
-    template<typename MatrixType, typename OtherDerived, bool IsSameType, bool IsDynamicSize>
+    template<typename MatrixType, typename OtherDerived, bool SwapPointers>
     friend struct ei_matrix_swap_impl;
 };
 
-template<typename MatrixType, typename OtherDerived,
+template <typename Derived, typename OtherDerived, bool IsVector>
-         bool IsSameType = ei_is_same_type<MatrixType, OtherDerived>::ret,
+struct ei_conservative_resize_like_impl
-         bool IsDynamicSize = MatrixType::SizeAtCompileTime==Dynamic>
+{
+  static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
+  {
+    // Note: Here is space for improvement. Basically, for conservativeResize(int,int),
+    // neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the
+    // dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or
+    // conservativeResize(NoChange_t, int cols). For these methods new static asserts like
+    // EIGEN_STATIC_ASSERT_DYNAMIC_ROWS and EIGEN_STATIC_ASSERT_DYNAMIC_COLS would be good.
+    EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived)
+    EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(OtherDerived)
+
+    typename MatrixBase<Derived>::PlainMatrixType tmp(other);
+    const int common_rows = std::min(tmp.rows(), _this.rows());
+    const int common_cols = std::min(tmp.cols(), _this.cols());
+    tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
+    _this.derived().swap(tmp);
+  }
+};
+
+template <typename Derived, typename OtherDerived>
+struct ei_conservative_resize_like_impl<Derived,OtherDerived,true>
+{
+  static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other)
+  {
+    // segment(...) will check whether Derived/OtherDerived are vectors!
+    typename MatrixBase<Derived>::PlainMatrixType tmp(other);
+    const int common_size = std::min<int>(_this.size(),tmp.size());
+    tmp.segment(0,common_size) = _this.segment(0,common_size);
+    _this.derived().swap(tmp);
+  }
+};
+
+template<typename MatrixType, typename OtherDerived, bool SwapPointers>
 struct ei_matrix_swap_impl
 {
   static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
@@ -729,7 +746,7 @@ struct ei_matrix_swap_impl
 };
 
 template<typename MatrixType, typename OtherDerived>
-struct ei_matrix_swap_impl<MatrixType, OtherDerived, true, true>
+struct ei_matrix_swap_impl<MatrixType, OtherDerived, true>
 {
   static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
   {
@@ -741,7 +758,8 @@ template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int
 template<typename OtherDerived>
 inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
 {
-  ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
+  enum { SwapPointers = ei_is_same_type<Matrix, OtherDerived>::ret && Base::SizeAtCompileTime==Dynamic };
+  ei_matrix_swap_impl<Matrix, OtherDerived, bool(SwapPointers)>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
 }
 
 /** \defgroup matrixtypedefs Global matrix typedefs

Eigen/src/Core/MatrixBase.h

@@ -26,46 +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
-  : public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
-                                     typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>
-{
-  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
@@ -93,11 +53,11 @@ template<typename Derived> struct AnyMatrixBase
   */
 template<typename Derived> class MatrixBase
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-  : 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
 {
   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*;
@@ -302,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)
+    Derived& operator=(const AnyMatrixBase<OtherDerived> &other);
-    { other.derived().evalToDense(derived()); return derived(); }
 
     template<typename OtherDerived>
-    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
+    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other);
-    { other.derived().addToDense(derived()); return derived(); }
 
     template<typename OtherDerived>
-    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
+    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other);
-    { other.derived().subToDense(derived()); return derived(); }
 
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
@@ -437,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;
@@ -676,8 +635,11 @@ template<typename Derived> class MatrixBase
     typename ei_traits<Derived>::Scalar minCoeff() const;
     typename ei_traits<Derived>::Scalar maxCoeff() const;
 
-    typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const;
+    typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col) const;
-    typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const;
+    typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col) const;
+
+    typename ei_traits<Derived>::Scalar minCoeff(int* index) const;
+    typename ei_traits<Derived>::Scalar maxCoeff(int* index) const;
 
     template<typename BinaryOp>
     typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type

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

Eigen/src/Core/StableNorm.h

@@ -56,7 +56,7 @@ MatrixBase<Derived>::stableNorm() const
 {
   const int blockSize = 4096;
   RealScalar scale = 0;
-  RealScalar invScale;
+  RealScalar invScale = 1;
   RealScalar ssq = 0; // sum of square
   enum {
     Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? ForceAligned : AsRequested

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;

Eigen/src/Core/VectorBlock.h

@@ -77,11 +77,12 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock
     typedef Block<VectorType,
                   ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
                   ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
-                  PacketAccess> Base;
+                  PacketAccess> _Base;
     enum {
       IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1
     };
   public:
+    _EIGEN_GENERIC_PUBLIC_INTERFACE(VectorBlock, _Base)
 
     using Base::operator=;
     using Base::operator+=;

Eigen/src/Core/Visitor.h

@@ -164,7 +164,7 @@ struct ei_functor_traits<ei_max_coeff_visitor<Scalar> > {
 /** \returns the minimum of all coefficients of *this
   * and puts in *row and *col its location.
   *
-  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  * \sa MatrixBase::minCoeff(int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
   */
 template<typename Derived>
 typename ei_traits<Derived>::Scalar
@@ -177,6 +177,22 @@ MatrixBase<Derived>::minCoeff(int* row, int* col) const
   return minVisitor.res;
 }
 
+/** \returns the minimum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::minCoeff(int*,int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::minCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_min_coeff_visitor<Scalar> minVisitor;
+  this->visit(minVisitor);
+  *index = (RowsAtCompileTime==1) ? minVisitor.col : minVisitor.row;
+  return minVisitor.res;
+}
+
 /** \returns the maximum of all coefficients of *this
   * and puts in *row and *col its location.
   *
@@ -193,5 +209,20 @@ MatrixBase<Derived>::maxCoeff(int* row, int* col) const
   return maxVisitor.res;
 }
 
+/** \returns the maximum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::minCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::maxCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::maxCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_max_coeff_visitor<Scalar> maxVisitor;
+  this->visit(maxVisitor);
+  *index = (RowsAtCompileTime==1) ? maxVisitor.col : maxVisitor.row;
+  return maxVisitor.res;
+}
 
 #endif // EIGEN_VISITOR_H

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:

Eigen/src/Core/util/XprHelper.h

@@ -217,7 +217,7 @@ template<unsigned int Flags> struct ei_are_flags_consistent
   * overloads for complex types */
 template<typename Derived,typename Scalar,typename OtherScalar,
          bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret >
-struct ei_special_scalar_op_base
+struct ei_special_scalar_op_base : public AnyMatrixBase<Derived>
 {
   // dummy operator* so that the
   // "using ei_special_scalar_op_base::operator*" compiles
@@ -225,7 +225,7 @@ struct ei_special_scalar_op_base
 };
 
 template<typename Derived,typename Scalar,typename OtherScalar>
-struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>
+struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true>  : public AnyMatrixBase<Derived>
 {
   const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived>
   operator*(const OtherScalar& scalar) const

Eigen/src/Eigenvalues/ComplexSchur.h

@@ -31,8 +31,15 @@
   *
   * \class ComplexShur
   *
-  * \brief Performs a complex Shur decomposition of a real or complex square matrix
+  * \brief Performs a complex Schur decomposition of a real or complex square matrix
   *
+  * Given a real or complex square matrix A, this class computes the Schur decomposition:
+  * \f$ A = U T U^*\f$ where U is a unitary complex matrix, and T is a complex upper
+  * triangular matrix.
+  *
+  * The diagonal of the matrix T corresponds to the eigenvalues of the matrix A.
+  *
+  * \sa class RealSchur, class EigenSolver
   */
 template<typename _MatrixType> class ComplexSchur
 {
@@ -42,41 +49,56 @@ template<typename _MatrixType> class ComplexSchur
     typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef std::complex<RealScalar> Complex;
     typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType;
+    enum {
+      Size = MatrixType::RowsAtCompileTime
+    };
 
-    /**
+    /** \brief Default Constructor.
-      * \brief Default Constructor.
       *
       * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via ComplexSchur::compute(const MatrixType&).
+      * perform decompositions via ComplexSchur::compute().
       */
-    ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false)
+    ComplexSchur(int size = Size==Dynamic ? 0 : Size)
+      : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
     {}
 
-    ComplexSchur(const MatrixType& matrix)
+    /** Constructor computing the Schur decomposition of the matrix \a matrix.
+      * If \a skipU is true, then the matrix U is not computed. */
+    ComplexSchur(const MatrixType& matrix, bool skipU = false)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matUisUptodate(false)
     {
-      compute(matrix);
+      compute(matrix, skipU);
     }
 
-    ComplexMatrixType matrixU() const
+    /** \returns a const reference to the matrix U of the respective Schur decomposition. */
+    const ComplexMatrixType& matrixU() const
     {
       ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
+      ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
       return m_matU;
     }
 
-    ComplexMatrixType matrixT() const
+    /** \returns a const reference to the matrix T of the respective Schur decomposition.
+      * Note that this function returns a plain square matrix. If you want to reference
+      * only the upper triangular part, use:
+      * \code schur.matrixT().triangularView<Upper>() \endcode. */
+    const ComplexMatrixType& matrixT() const
     {
       ei_assert(m_isInitialized && "ComplexShur is not initialized.");
       return m_matT;
     }
 
-    void compute(const MatrixType& matrix);
+    /** Computes the Schur decomposition of the matrix \a matrix.
+      * If \a skipU is true, then the matrix U is not computed. */
+    void compute(const MatrixType& matrix, bool skipU = false);
 
   protected:
     ComplexMatrixType m_matT, m_matU;
     bool m_isInitialized;
+    bool m_matUisUptodate;
 };
 
 /** Computes the principal value of the square root of the complex \a z. */
@@ -117,17 +139,20 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
 }
 
 template<typename MatrixType>
-void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
+void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
 {
   // this code is inspired from Jampack
+
+  m_matUisUptodate = false;
   assert(matrix.cols() == matrix.rows());
   int n = matrix.cols();
 
   // Reduce to Hessenberg form
+  // TODO skip Q if skipU = true
   HessenbergDecomposition<MatrixType> hess(matrix);
 
   m_matT = hess.matrixH();
-  m_matU = hess.matrixQ();
+  if(!skipU)  m_matU = hess.matrixQ();
 
   int iu = m_matT.cols() - 1;
   int il;
@@ -206,7 +231,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
     {
       m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
       m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
-      m_matU.applyOnTheRight(i, i+1, rot);
+      if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);
 
       if(i != iu-1)
       {
@@ -232,6 +257,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
   */
 
   m_isInitialized = true;
+  m_matUisUptodate = !skipU;
 }
 
 #endif // EIGEN_COMPLEX_SCHUR_H

Eigen/src/Eigenvalues/HessenbergDecomposition.h

@@ -88,14 +88,14 @@ template<typename _MatrixType> class HessenbergDecomposition
       _compute(m_matrix, m_hCoeffs);
     }
 
-    /** \returns the householder coefficients allowing to
+    /** \returns a const reference to the householder coefficients allowing to
       * reconstruct the matrix Q from the packed data.
       *
       * \sa packedMatrix()
       */
-    CoeffVectorType householderCoefficients() const { return m_hCoeffs; }
+    const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; }
 
-    /** \returns the internal result of the decomposition.
+    /** \returns a const reference to the internal representation of the decomposition.
       *
       * The returned matrix contains the following information:
       *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H

Eigen/src/Geometry/Transform.h

@@ -395,7 +395,7 @@ public:
   Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position,
     const OrientationType& orientation, const MatrixBase<ScaleDerived> &scale);
 
-  inline const MatrixType inverse(TransformTraits traits = (TransformTraits)Mode) const;
+  inline Transform inverse(TransformTraits traits = (TransformTraits)Mode) const;
 
   /** \returns a const pointer to the column major internal matrix */
   const Scalar* data() const { return m_matrix.data(); }
@@ -874,7 +874,7 @@ Transform<Scalar,Dim,Mode>::fromPositionOrientationScale(const MatrixBase<Positi
 
 /** \nonstableyet
   *
-  * \returns the inverse transformation matrix according to some given knowledge
+  * \returns the inverse transformation according to some given knowledge
   * on \c *this.
   *
   * \param traits allows to optimize the inversion process when the transformion
@@ -892,38 +892,38 @@ Transform<Scalar,Dim,Mode>::fromPositionOrientationScale(const MatrixBase<Positi
   * \sa MatrixBase::inverse()
   */
 template<typename Scalar, int Dim, int Mode>
-const typename Transform<Scalar,Dim,Mode>::MatrixType
+Transform<Scalar,Dim,Mode>
 Transform<Scalar,Dim,Mode>::inverse(TransformTraits hint) const
 {
+  Transform res;
   if (hint == Projective)
   {
-    return m_matrix.inverse();
+    res.matrix() = m_matrix.inverse();
   }
   else
   {
-    MatrixType res;
     if (hint == Isometry)
     {
-      res.template corner<Dim,Dim>(TopLeft) = linear().transpose();
+      res.matrix().template corner<Dim,Dim>(TopLeft) = linear().transpose();
     }
     else if(hint&Affine)
     {
-      res.template corner<Dim,Dim>(TopLeft) = linear().inverse();
+      res.matrix().template corner<Dim,Dim>(TopLeft) = linear().inverse();
     }
     else
     {
       ei_assert(false && "Invalid transform traits in Transform::Inverse");
     }
     // translation and remaining parts
-    res.template corner<Dim,1>(TopRight) = - res.template corner<Dim,Dim>(TopLeft) * translation();
+    res.matrix().template corner<Dim,1>(TopRight) = - res.matrix().template corner<Dim,Dim>(TopLeft) * translation();
     if(int(Mode)!=int(AffineCompact))
     {
-      res.template block<1,Dim>(Dim,0).setZero();
+      res.matrix().template block<1,Dim>(Dim,0).setZero();
-      res.coeffRef(Dim,Dim) = 1;
+      res.matrix().coeffRef(Dim,Dim) = 1;
+    }
   }
   return res;
 }
-}
 
 /*****************************************************
 *** Specializations of take affine part            ***

Eigen/src/Householder/HouseholderSequence.h

@@ -0,0 +1,168 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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_HOUSEHOLDER_SEQUENCE_H
+#define EIGEN_HOUSEHOLDER_SEQUENCE_H
+
+/** \ingroup Householder_Module
+  * \householder_module
+  * \class HouseholderSequence
+  * \brief Represents a sequence of householder reflections with decreasing size
+  *
+  * This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$
+  * where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$,
+  * \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$
+  * and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i].
+  *
+  * Typical usages are listed below, where H is a HouseholderSequence:
+  * \code
+  * A.applyOnTheRight(H);             // A = A * H
+  * A.applyOnTheLeft(H);              // A = H * A
+  * A.applyOnTheRight(H.adjoint());   // A = A * H^*
+  * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
+  * MatrixXd Q = H;                   // conversion to a dense matrix
+  * \endcode
+  * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate.
+  *
+  * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
+  */
+
+template<typename VectorsType, typename CoeffsType>
+struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> >
+{
+  typedef typename VectorsType::Scalar Scalar;
+  enum {
+    RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+    ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+    MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    Flags = 0
+  };
+};
+
+template<typename VectorsType, typename CoeffsType> class HouseholderSequence
+  : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> >
+{
+    typedef typename VectorsType::Scalar Scalar;
+  public:
+
+    typedef HouseholderSequence<VectorsType,
+      typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+        NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
+        CoeffsType>::ret> ConjugateReturnType;
+
+    HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
+      : m_vectors(v), m_coeffs(h), m_trans(trans)
+    {}
+
+    int rows() const { return m_vectors.rows(); }
+    int cols() const { return m_vectors.rows(); }
+
+    HouseholderSequence transpose() const
+    { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); }
+
+    ConjugateReturnType conjugate() const
+    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
+
+    ConjugateReturnType adjoint() const
+    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
+
+    ConjugateReturnType inverse() const { return adjoint(); }
+
+    /** \internal */
+    template<typename DestType> void evalTo(DestType& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows());
+      int length = m_vectors.rows();
+      dst.setIdentity();
+      Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
+      for(int k = vecs-1; k >= 0; --k)
+      {
+        if(m_trans)
+          dst.corner(BottomRight, length-k, length-k)
+          .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
+        else
+          dst.corner(BottomRight, length-k, length-k)
+            .applyHouseholderOnTheLeft(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k));
+      }
+    }
+
+    /** \internal */
+    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
+      int length = m_vectors.rows();                          // size of the largest householder vector
+      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows());
+      for(int k = 0; k < vecs; ++k)
+      {
+        int actual_k = m_trans ? vecs-k-1 : k;
+        dst.corner(BottomRight, dst.rows(), length-k)
+           .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
+      }
+    }
+
+    /** \internal */
+    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
+      int length = m_vectors.rows();                          // size of the largest householder vector
+      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
+      for(int k = 0; k < vecs; ++k)
+      {
+        int actual_k = m_trans ? k : vecs-k-1;
+        dst.corner(BottomRight, length-actual_k, dst.cols())
+           .applyHouseholderOnTheLeft(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
+      }
+    }
+
+    template<typename OtherDerived>
+    typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other) const
+    {
+      typename OtherDerived::PlainMatrixType res(other);
+      applyThisOnTheLeft(res);
+      return res;
+    }
+
+    template<typename OtherDerived> friend
+    typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h)
+    {
+      typename OtherDerived::PlainMatrixType res(other);
+      h.applyThisOnTheRight(res);
+      return res;
+    }
+
+  protected:
+
+    typename VectorsType::Nested m_vectors;
+    typename CoeffsType::Nested m_coeffs;
+    bool m_trans;
+};
+
+template<typename VectorsType, typename CoeffsType>
+HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
+{
+  return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
+}
+
+#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

Eigen/src/Jacobi/Jacobi.h

@@ -123,7 +123,7 @@ bool PlanarRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
 }
 
 /** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 selfadjoint matrix
-  * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ \overline \text{this}_{pq} & \text{this}_{qq} \end{array} \right )\f$ yields
+  * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields
   * a diagonal matrix \f$ A = J^* B J \f$
   *
   * Example: \include Jacobi_makeJacobi.cpp

Eigen/src/LU/PartialLU.h

@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2006-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
@@ -232,7 +233,7 @@ struct ei_partial_lu_impl
     for(int k = 0; k < size; ++k)
     {
       int row_of_biggest_in_col;
-      lu.block(k,k,rows-k,1).cwise().abs().maxCoeff(&row_of_biggest_in_col);
+      lu.col(k).end(rows-k).cwise().abs().maxCoeff(&row_of_biggest_in_col);
       row_of_biggest_in_col += k;
 
       row_transpositions[k] = row_of_biggest_in_col;

Eigen/src/QR/ColPivotingHouseholderQR.h

@@ -62,6 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
     typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
     typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
+    typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
     /**
     * \brief Default Constructor.
@@ -99,7 +100,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
     template<typename OtherDerived, typename ResultType>
     bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixQType matrixQ(void) const;
+    HouseholderSequenceType matrixQ(void) const;
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
@@ -306,6 +307,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
     cols_transpositions.coeffRef(k) = biggest_col_in_corner;
     if(k != biggest_col_in_corner) {
       m_qr.col(k).swap(m_qr.col(biggest_col_in_corner));
+      std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_in_corner));
       ++number_of_transpositions;
     }
 
@@ -354,13 +356,8 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
 
   typename OtherDerived::PlainMatrixType c(b);
 
-  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
+  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
-  for (int k = 0; k < m_rank; ++k)
+  c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose());
-  {
-    int remainingSize = rows-k;
-    c.corner(BottomRight, remainingSize, cols)
-     .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
-  }
 
   if(!isSurjective())
   {
@@ -380,25 +377,12 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
   return true;
 }
 
-/** \returns the matrix Q */
+/** \returns the matrix Q as a sequence of householder transformations */
 template<typename MatrixType>
-typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
+typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
-  // compute the product H'_0 H'_1 ... H'_n-1,
+  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
-  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
-  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
-  int rows = m_qr.rows();
-  int cols = m_qr.cols();
-  int size = std::min(rows,cols);
-  MatrixQType res = MatrixQType::Identity(rows, rows);
-  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
-  for (int k = size-1; k >= 0; k--)
-  {
-    res.block(k, k, rows-k, rows-k)
-       .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
-  }
-  return res;
 }
 
 #endif // EIGEN_HIDE_HEAVY_CODE

Eigen/src/QR/HouseholderQR.h

@@ -62,6 +62,7 @@ template<typename MatrixType> class HouseholderQR
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
     typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
+    typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
     /**
     * \brief Default Constructor.
@@ -97,7 +98,12 @@ template<typename MatrixType> class HouseholderQR
     template<typename OtherDerived, typename ResultType>
     void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixQType matrixQ(void) const;
+    MatrixQType matrixQ() const;
+
+    HouseholderSequenceType matrixQAsHouseholderSequence() const
+    {
+      return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
+    }
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       * in a LAPACK-compatible way.
@@ -206,15 +212,7 @@ void HouseholderQR<MatrixType>::solve(
   result->resize(rows, cols);
 
   *result = b;
-  
+  result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse());
-  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
-  for (int k = 0; k < cols; ++k)
-  {
-    int remainingSize = rows-k;
-
-    result->corner(BottomRight, remainingSize, cols)
-           .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
-  }
 
   const int rank = std::min(result->rows(), result->cols());
   m_qr.corner(TopLeft, rank, rank)
@@ -227,20 +225,7 @@ template<typename MatrixType>
 typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
-  // compute the product H'_0 H'_1 ... H'_n-1,
+  return matrixQAsHouseholderSequence();
-  // where H_k is the k-th Householder transformation I - h_k v_k v_k'
-  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
-  int rows = m_qr.rows();
-  int cols = m_qr.cols();
-  int size = std::min(rows,cols);
-  MatrixQType res = MatrixQType::Identity(rows, rows);
-  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
-  for (int k = size-1; k >= 0; k--)
-  {
-    res.block(k, k, rows-k, rows-k)
-       .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
-  }
-  return res;
 }
 
 #endif // EIGEN_HIDE_HEAVY_CODE

Eigen/src/SVD/JacobiSVD.h

@@ -25,6 +25,22 @@
 #ifndef EIGEN_JACOBISVD_H
 #define EIGEN_JACOBISVD_H
 
+// forward declarations (needed by ICC)
+template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
+struct ei_svd_precondition_2x2_block_to_be_real;
+
+template<typename MatrixType, unsigned int Options,
+         bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0
+                                         && (MatrixType::RowsAtCompileTime==Dynamic
+                                             || (MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime))>
+struct ei_svd_precondition_if_more_rows_than_cols;
+
+template<typename MatrixType, unsigned int Options,
+         bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0
+                                         && (MatrixType::ColsAtCompileTime==Dynamic
+                                             || (MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime))>
+struct ei_svd_precondition_if_more_cols_than_rows;
+
 /** \ingroup SVD_Module
   * \nonstableyet
   *
@@ -118,8 +134,8 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD
     friend struct ei_svd_precondition_if_more_cols_than_rows;
 };
 
-template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
+template<typename MatrixType, unsigned int Options>
-struct ei_svd_precondition_2x2_block_to_be_real
+struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, false>
 {
   typedef JacobiSVD<MatrixType, Options> SVD;
   static void run(typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {}
@@ -195,10 +211,7 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q,
   *j_left  = rot1 * j_right->transpose();
 }
 
-template<typename MatrixType, unsigned int Options,
+template<typename MatrixType, unsigned int Options, bool PossiblyMoreRowsThanCols>
-         bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0
-                                         && (MatrixType::RowsAtCompileTime==Dynamic
-                                             || MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime)>
 struct ei_svd_precondition_if_more_rows_than_cols
 {
   typedef JacobiSVD<MatrixType, Options> SVD;
@@ -231,10 +244,7 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true>
   }
 };
 
-template<typename MatrixType, unsigned int Options,
+template<typename MatrixType, unsigned int Options, bool PossiblyMoreColsThanRows>
-         bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0
-                                         && (MatrixType::ColsAtCompileTime==Dynamic
-                                            || MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime)>
 struct ei_svd_precondition_if_more_cols_than_rows
 {
   typedef JacobiSVD<MatrixType, Options> SVD;

Eigen/src/Sparse/CholmodSupport.h

@@ -99,7 +99,7 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
   res.nrow   = mat.rows();
   res.ncol   = mat.cols();
   res.nzmax  = res.nrow * res.ncol;
-  res.d      = mat.derived().stride();
+  res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().stride();
   res.x      = mat.derived().data();
   res.z      = 0;
 
@@ -157,7 +157,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
     inline const typename Base::CholMatrixType& matrixL(void) const;
 
     template<typename Derived>
-    void solveInPlace(MatrixBase<Derived> &b) const;
+    bool solveInPlace(MatrixBase<Derived> &b) const;
 
     void compute(const MatrixType& matrix);
 
@@ -216,7 +216,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const
 
 template<typename MatrixType>
 template<typename Derived>
-void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+bool SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
 {
   const int size = m_cholmodFactor->n;
   ei_assert(size==b.rows());
@@ -228,9 +228,16 @@ void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
   // as long as our own triangular sparse solver is not fully optimal,
   // let's use CHOLMOD's one:
   cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
-  cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
+  //cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
+  cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod);
+  if(!x)
+  {
+    std::cerr << "Eigen: cholmod_solve failed\n";
+    return false;
+  }
   b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
   cholmod_free_dense(&x, &m_cholmod);
+  return true;
 }
 
 #endif // EIGEN_CHOLMODSUPPORT_H

Eigen/src/Sparse/SuperLUSupport.h

@@ -161,7 +161,7 @@ struct SluMatrix : SuperMatrix
     res.nrow      = mat.rows();
     res.ncol      = mat.cols();
 
-    res.storage.lda       = mat.stride();
+    res.storage.lda       = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.stride();
     res.storage.values    = mat.data();
     return res;
   }

cmake/EigenTesting.cmake

@@ -153,11 +153,17 @@ macro(ei_init_testing)
 endmacro(ei_init_testing)
 
 if(CMAKE_COMPILER_IS_GNUCXX)
+  option(EIGEN_COVERAGE_TESTING "Enable/disable gcov" OFF)
+  if(EIGEN_COVERAGE_TESTING)
+    set(COVERAGE_FLAGS "-fprofile-arcs -ftest-coverage")
+  else(EIGEN_COVERAGE_TESTING)
+    set(COVERAGE_FLAGS "")
+  endif(EIGEN_COVERAGE_TESTING)
   if(CMAKE_SYSTEM_NAME MATCHES Linux)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2")
-    set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -O2 -g2")
+    set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2")
-    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fno-inline-functions")
+    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${COVERAGE_FLAGS} -fno-inline-functions")
-    set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -g2")
+    set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${COVERAGE_FLAGS} -O0 -g2")
   endif(CMAKE_SYSTEM_NAME MATCHES Linux)
   set(EI_OFLAG "-O2")
 elseif(MSVC)

doc/C01_QuickStartGuide.dox

@@ -129,7 +129,7 @@ The default constructor leaves coefficients uninitialized. Any dynamic size is s
 Matrix3f A; // construct 3x3 matrix with uninitialized coefficients
 A(0,0) = 5; // OK
 MatrixXf B; // construct 0x0 matrix without allocating anything
-A(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address
+B(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address
 \endcode
 
 In the above example, B is an uninitialized matrix. What to do with such a matrix? You can call resize() on it, or you can assign another matrix to it. Like this:
@@ -261,7 +261,7 @@ v = 6 6 6
 
 \subsection TutorialCasting Casting
 
-In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitely casted to another one using the template MatrixBase::cast() function:
+In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitly casted to another one using the template MatrixBase::cast() function:
 \code
 Matrix3d md(1,2,3);
 Matrix3f mf = md.cast<float>();
@@ -328,7 +328,7 @@ In short, all arithmetic operators can be used right away as in the following ex
 mat4 -= mat1*1.5 + mat2 * (mat3/4);
 \endcode
 which includes two matrix scalar products ("mat1*1.5" and "mat3/4"), a matrix-matrix product ("mat2 * (mat3/4)"),
-a matrix addition ("+") and substraction with assignment ("-=").
+a matrix addition ("+") and subtraction with assignment ("-=").
 
 <table class="tutorial_code">
 <tr><td>
@@ -464,7 +464,7 @@ mat = 2 7 8
 
 Also note that maxCoeff and minCoeff can takes optional arguments returning the coordinates of the respective min/max coeff: \link MatrixBase::maxCoeff(int*,int*) const maxCoeff(int* i, int* j) \endlink, \link MatrixBase::minCoeff(int*,int*) const minCoeff(int* i, int* j) \endlink.
 
-<span class="note">\b Side \b note: The all() and any() functions are especially useful in combinaison with coeff-wise comparison operators (\ref CwiseAll "example").</span>
+<span class="note">\b Side \b note: The all() and any() functions are especially useful in combination with coeff-wise comparison operators (\ref CwiseAll "example").</span>
 
 
 
@@ -578,7 +578,7 @@ vec1.normalize();\endcode
 
 <a href="#" class="top">top</a>\section TutorialCoreTriangularMatrix Dealing with triangular matrices
 
-Currently, Eigen does not provide any explcit triangular matrix, with storage class. Instead, we
+Currently, Eigen does not provide any explicit triangular matrix, with storage class. Instead, we
 can reference a triangular part of a square matrix or expression to perform special treatment on it.
 This is achieved by the class TriangularView and the MatrixBase::triangularView template function.
 Note that the opposite triangular part of  the matrix is never referenced, and so it can, e.g., store
@@ -595,12 +595,12 @@ m.triangularView<Eigen::LowerTriangular>()
 m.triangularView<Eigen::UnitLowerTriangular>()\endcode
 </td></tr>
 <tr><td>
-Writting to a specific triangular part:\n (only the referenced triangular part is evaluated)
+Writing to a specific triangular part:\n (only the referenced triangular part is evaluated)
 </td><td>\code
 m1.triangularView<Eigen::LowerTriangular>() = m2 + m3 \endcode
 </td></tr>
 <tr><td>
-Convertion to a dense matrix setting the opposite triangular part to zero:
+Conversion to a dense matrix setting the opposite triangular part to zero:
 </td><td>\code
 m2 = m1.triangularView<Eigen::UnitUpperTriangular>()\endcode
 </td></tr>

test/CMakeLists.txt

@@ -94,6 +94,7 @@ ei_add_test(basicstuff)
 ei_add_test(linearstructure)
 ei_add_test(cwiseop)
 ei_add_test(redux)
+ei_add_test(visitor)
 ei_add_test(product_small)
 ei_add_test(product_large ${EI_OFLAG})
 ei_add_test(product_extra ${EI_OFLAG})
@@ -115,6 +116,7 @@ ei_add_test(product_trmv ${EI_OFLAG})
 ei_add_test(product_trmm ${EI_OFLAG})
 ei_add_test(product_trsm ${EI_OFLAG})
 ei_add_test(product_notemporary ${EI_OFLAG})
+ei_add_test(stable_norm)
 ei_add_test(bandmatrix)
 ei_add_test(cholesky " " "${GSL_LIBRARIES}")
 ei_add_test(lu ${EI_OFLAG})

test/adjoint.cpp

@@ -72,13 +72,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   if(NumTraits<Scalar>::HasFloatingPoint)
     VERIFY_IS_APPROX(v1.squaredNorm(),                v1.norm() * v1.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)),  static_cast<RealScalar>(1));
-  if(NumTraits<Scalar>::HasFloatingPoint)
-  {
-    VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(),         static_cast<RealScalar>(1));
-    VERIFY_IS_APPROX(v1.norm(),                       v1.stableNorm());
-    VERIFY_IS_APPROX(v1.blueNorm(),                   v1.stableNorm());
-    VERIFY_IS_APPROX(v1.hypotNorm(),                  v1.stableNorm());
-  }
 
   // check compatibility of dot and adjoint
   VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps));

test/cholesky.cpp

@@ -149,16 +149,16 @@ void test_cholesky()
 {
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
-//     CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
+    CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
-//     CALL_SUBTEST( cholesky(Matrix2d()) );
+    CALL_SUBTEST( cholesky(Matrix2d()) );
-//     CALL_SUBTEST( cholesky(Matrix3f()) );
+    CALL_SUBTEST( cholesky(Matrix3f()) );
-//     CALL_SUBTEST( cholesky(Matrix4d()) );
+    CALL_SUBTEST( cholesky(Matrix4d()) );
     CALL_SUBTEST( cholesky(MatrixXd(200,200)) );
     CALL_SUBTEST( cholesky(MatrixXcd(100,100)) );
   }
 
-//   CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
+  CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
-//   CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
+  CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
-//   CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
+  CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
-//   CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
+  CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
 }

test/conservative_resize.cpp

@@ -65,7 +65,7 @@ void run_matrix_tests()
     const int rows = ei_random<int>(50,75);
     const int cols = ei_random<int>(50,75);
     m = n = MatrixType::Random(50,50);
-    m.conservativeResize(rows,cols,true);
+    m.conservativeResizeLike(MatrixType::Zero(rows,cols));
     VERIFY_IS_APPROX(m.block(0,0,n.rows(),n.cols()), n);
     VERIFY( rows<=50 || m.block(50,0,rows-50,cols).sum() == Scalar(0) );
     VERIFY( cols<=50 || m.block(0,50,rows,cols-50).sum() == Scalar(0) );
@@ -102,7 +102,7 @@ void run_vector_tests()
   {
     const int size = ei_random<int>(50,100);
     m = n = MatrixType::Random(50);
-    m.conservativeResize(size,true);
+    m.conservativeResizeLike(MatrixType::Zero(size));
     VERIFY_IS_APPROX(m.segment(0,50), n);
     VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) );
   }

test/geo_orthomethods.cpp

@@ -86,10 +86,10 @@ template<typename Scalar, int Size> void orthomethods(int size=Size)
   VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
   VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1));
 
-  if (size>3)
+  if (size>=3)
   {
-    v0.template start<3>().setZero();
+    v0.template start<2>().setZero();
-    v0.end(size-3).setRandom();
+    v0.end(size-2).setRandom();
 
     VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1));
     VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1));

test/geo_transformations.cpp

@@ -296,10 +296,10 @@ template<typename Scalar, int Mode> void transformations(void)
     t0.setIdentity();
     t0.translate(v0);
     t0.linear().setRandom();
-    VERIFY_IS_APPROX(t0.inverse(Affine), t0.matrix().inverse());
+    VERIFY_IS_APPROX(t0.inverse(Affine).matrix(), t0.matrix().inverse());
     t0.setIdentity();
     t0.translate(v0).rotate(q1);
-    VERIFY_IS_APPROX(t0.inverse(Isometry), t0.matrix().inverse());
+    VERIFY_IS_APPROX(t0.inverse(Isometry).matrix(), t0.matrix().inverse());
   }
 
   // test extract rotation and aligned scaling

test/householder.cpp

@@ -84,6 +84,9 @@ template<typename MatrixType> void householder(const MatrixType& m)
   VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0)));
   VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha);
+
+  // test householder sequence
+  // TODO test HouseholderSequence
 }
 
 void test_householder()

test/main.h

@@ -40,6 +40,11 @@
 
 #define DEFAULT_REPEAT 10
 
+#ifdef __ICC
+// disable warning #279: controlling expression is constant
+#pragma warning disable 279
+#endif
+
 namespace Eigen
 {
   static std::vector<std::string> g_test_stack;

test/mixingtypes.cpp

@@ -175,9 +175,9 @@ void test_mixingtypes()
 {
   // check that our operator new is indeed called:
   CALL_SUBTEST(mixingtypes<3>());
-//   CALL_SUBTEST(mixingtypes<4>());
+  CALL_SUBTEST(mixingtypes<4>());
-//   CALL_SUBTEST(mixingtypes<Dynamic>(20));
+  CALL_SUBTEST(mixingtypes<Dynamic>(20));
-//
+
-//   CALL_SUBTEST(mixingtypes_small<4>());
+  CALL_SUBTEST(mixingtypes_small<4>());
-//   CALL_SUBTEST(mixingtypes_large(20));
+  CALL_SUBTEST(mixingtypes_large(20));
 }

test/product_extra.cpp

@@ -104,13 +104,24 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
   VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()),
                    (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval());
 
+  // test the vector-matrix product with non aligned starts
+  int i = ei_random<int>(0,m1.rows()-2);
+  int j = ei_random<int>(0,m1.cols()-2);
+  int r = ei_random<int>(1,m1.rows()-i);
+  int c = ei_random<int>(1,m1.cols()-j);
+  int i2 = ei_random<int>(0,m1.rows()-1);
+  int j2 = ei_random<int>(0,m1.cols()-1);
+
+  VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval());
+  VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval());
+
 }
 
 void test_product_extra()
 {
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
+    CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) );
     CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
-    CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,50), ei_random<int>(1,50))) );
+    CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) );
   }
 }

test/qr_colpivoting.cpp

@@ -117,8 +117,6 @@ template<typename MatrixType> void qr_verify_assert()
 void test_qr_colpivoting()
 {
   for(int i = 0; i < 1; i++) {
-    // FIXME : very weird bug here
-//     CALL_SUBTEST( qr(Matrix2f()) );
     CALL_SUBTEST( qr<MatrixXf>() );
     CALL_SUBTEST( qr<MatrixXd>() );
     CALL_SUBTEST( qr<MatrixXcd>() );

test/redux.cpp

@@ -120,7 +120,7 @@ void test_redux()
     CALL_SUBTEST( matrixRedux(MatrixXi(8, 12)) );
   }
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST( vectorRedux(VectorXf(5)) );
+    CALL_SUBTEST( vectorRedux(Vector4f()) );
     CALL_SUBTEST( vectorRedux(VectorXd(10)) );
     CALL_SUBTEST( vectorRedux(VectorXf(33)) );
   }

test/stable_norm.cpp

@@ -0,0 +1,79 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// 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/>.
+
+#include "main.h"
+
+template<typename MatrixType> void stable_norm(const MatrixType& m)
+{
+  /* this test covers the following files:
+     StableNorm.h
+  */
+
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+
+  int rows = m.rows();
+  int cols = m.cols();
+
+  Scalar big = ei_random<Scalar>() * std::numeric_limits<RealScalar>::max() * RealScalar(1e-4);
+  Scalar small = static_cast<RealScalar>(1)/big;
+
+  MatrixType  vzero = MatrixType::Zero(rows, cols),
+              vrand = MatrixType::Random(rows, cols),
+              vbig(rows, cols),
+              vsmall(rows,cols);
+
+  vbig.fill(big);
+  vsmall.fill(small);
+
+  VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
+  VERIFY_IS_APPROX(vrand.stableNorm(),      vrand.norm());
+  VERIFY_IS_APPROX(vrand.blueNorm(),        vrand.norm());
+  VERIFY_IS_APPROX(vrand.hypotNorm(),       vrand.norm());
+
+  RealScalar size = static_cast<RealScalar>(m.size());
+
+  // test overflow
+  VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vbig.norm()),   ei_sqrt(size)*big); // here the default norm must fail
+  VERIFY_IS_APPROX(static_cast<Scalar>(vbig.stableNorm()), ei_sqrt(size)*big);
+  VERIFY_IS_APPROX(static_cast<Scalar>(vbig.blueNorm()),   ei_sqrt(size)*big);
+  VERIFY_IS_APPROX(static_cast<Scalar>(vbig.hypotNorm()),  ei_sqrt(size)*big);
+
+  // test underflow
+  VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vsmall.norm()),   ei_sqrt(size)*small); // here the default norm must fail
+  VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.stableNorm()), ei_sqrt(size)*small);
+  VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.blueNorm()),   ei_sqrt(size)*small);
+  VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.hypotNorm()),  ei_sqrt(size)*small);
+}
+
+void test_stable_norm()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( stable_norm(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( stable_norm(Vector4d()) );
+    CALL_SUBTEST( stable_norm(VectorXd(ei_random<int>(10,2000))) );
+    CALL_SUBTEST( stable_norm(VectorXf(ei_random<int>(10,2000))) );
+    CALL_SUBTEST( stable_norm(VectorXcd(ei_random<int>(10,2000))) );
+  }
+}

test/visitor.cpp

@@ -0,0 +1,131 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 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 matrixVisitor(const MatrixType& p)
+{
+  typedef typename MatrixType::Scalar Scalar;
+
+  int rows = p.rows();
+  int cols = p.cols();
+
+  // construct a random matrix where all coefficients are different
+  MatrixType m;
+  m = MatrixType::Random(rows, cols);
+  for(int i = 0; i < m.size(); i++)
+    for(int i2 = 0; i2 < i; i2++)
+      while(m(i) == m(i2)) // yes, ==
+        m(i) = ei_random<Scalar>();
+  
+  Scalar minc = Scalar(1000), maxc = Scalar(-1000);
+  int minrow,mincol,maxrow,maxcol;
+  for(int j = 0; j < cols; j++)
+  for(int i = 0; i < rows; i++)
+  {
+    if(m(i,j) < minc)
+    {
+      minc = m(i,j);
+      minrow = i;
+      mincol = j;
+    }
+    if(m(i,j) > maxc)
+    {
+      maxc = m(i,j);
+      maxrow = i;
+      maxcol = j;
+    }
+  }
+  int eigen_minrow, eigen_mincol, eigen_maxrow, eigen_maxcol;
+  Scalar eigen_minc, eigen_maxc;
+  eigen_minc = m.minCoeff(&eigen_minrow,&eigen_mincol);
+  eigen_maxc = m.maxCoeff(&eigen_maxrow,&eigen_maxcol);
+  VERIFY(minrow == eigen_minrow);
+  VERIFY(maxrow == eigen_maxrow);
+  VERIFY(mincol == eigen_mincol);
+  VERIFY(maxcol == eigen_maxcol);
+  VERIFY_IS_APPROX(minc, eigen_minc);
+  VERIFY_IS_APPROX(maxc, eigen_maxc);
+  VERIFY_IS_APPROX(minc, m.minCoeff());
+  VERIFY_IS_APPROX(maxc, m.maxCoeff());
+}
+
+template<typename VectorType> void vectorVisitor(const VectorType& w)
+{
+  typedef typename VectorType::Scalar Scalar;
+
+  int size = w.size();
+
+  // construct a random vector where all coefficients are different
+  VectorType v;
+  v = VectorType::Random(size);
+  for(int i = 0; i < size; i++)
+    for(int i2 = 0; i2 < i; i2++)
+      while(v(i) == v(i2)) // yes, ==
+        v(i) = ei_random<Scalar>();
+  
+  Scalar minc = Scalar(1000), maxc = Scalar(-1000);
+  int minidx,maxidx;
+  for(int i = 0; i < size; i++)
+  {
+    if(v(i) < minc)
+    {
+      minc = v(i);
+      minidx = i;
+    }
+    if(v(i) > maxc)
+    {
+      maxc = v(i);
+      maxidx = i;
+    }
+  }
+  int eigen_minidx, eigen_maxidx;
+  Scalar eigen_minc, eigen_maxc;
+  eigen_minc = v.minCoeff(&eigen_minidx);
+  eigen_maxc = v.maxCoeff(&eigen_maxidx);
+  VERIFY(minidx == eigen_minidx);
+  VERIFY(maxidx == eigen_maxidx);
+  VERIFY_IS_APPROX(minc, eigen_minc);
+  VERIFY_IS_APPROX(maxc, eigen_maxc);
+  VERIFY_IS_APPROX(minc, v.minCoeff());
+  VERIFY_IS_APPROX(maxc, v.maxCoeff());
+}
+
+void test_visitor()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( matrixVisitor(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( matrixVisitor(Matrix2f()) );
+    CALL_SUBTEST( matrixVisitor(Matrix4d()) );
+    CALL_SUBTEST( matrixVisitor(MatrixXd(8, 12)) );
+    CALL_SUBTEST( matrixVisitor(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
+    CALL_SUBTEST( matrixVisitor(MatrixXi(8, 12)) );
+  }
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( vectorVisitor(Vector4f()) );
+    CALL_SUBTEST( vectorVisitor(VectorXd(10)) );
+    CALL_SUBTEST( vectorVisitor(RowVectorXd(10)) );
+    CALL_SUBTEST( vectorVisitor(VectorXf(33)) );
+  }
+}

unsupported/Eigen/AlignedVector3

@@ -88,7 +88,7 @@ template<typename _Scalar> class AlignedVector3
     {}
 
     inline AlignedVector3(const AlignedVector3& other)
-      : m_coeffs(other.m_coeffs)
+      : Base(), m_coeffs(other.m_coeffs)
     {}
 
     template<typename XprType, int Size=XprType::SizeAtCompileTime>
@@ -152,7 +152,6 @@ template<typename _Scalar> class AlignedVector3
     {
       ei_assert(m_coeffs.w()==Scalar(0));
       ei_assert(other.m_coeffs.w()==Scalar(0));
-      Scalar r = m_coeffs.dot(other.m_coeffs);
       return m_coeffs.dot(other.m_coeffs);
     }
 

unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h

@@ -25,6 +25,10 @@
 #ifndef EIGEN_MATRIX_EXPONENTIAL
 #define EIGEN_MATRIX_EXPONENTIAL
 
+#ifdef _MSC_VER
+  template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); }
+#endif
+
 /** \brief Compute the matrix exponential. 
  *
  * \param M      matrix whose exponential is to be computed. 
@@ -61,260 +65,243 @@ template <typename Derived>
 EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M, 
 					       typename MatrixBase<Derived>::PlainMatrixType* result);
 
-
+/** \brief Class for computing the matrix exponential.*/
-/** \internal \brief Internal helper functions for computing the
- *  matrix exponential.
- */
-namespace MatrixExponentialInternal {
-
-#ifdef _MSC_VER
-  template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); }
-#endif
-
-  /** \internal \brief Compute the (3,3)-Pad&eacute; approximant to
-   *  the exponential.
-   *  
-   *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
-   *  approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
-   *
-   *  \param M   Argument of matrix exponential
-   *  \param Id  Identity matrix of same size as M
-   *  \param tmp Temporary storage, to be provided by the caller
-   *  \param M2  Temporary storage, to be provided by the caller
-   *  \param U   Even-degree terms in numerator of Pad&eacute; approximant
-   *  \param V   Odd-degree terms in numerator of Pad&eacute; approximant
-   */
 template <typename MatrixType>
-  EIGEN_STRONG_INLINE void pade3(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, 
+class MatrixExponential {
-				 MatrixType& M2, MatrixType& U, MatrixType& V)
-  {
-    typedef typename ei_traits<MatrixType>::Scalar Scalar;
-    const Scalar b[] = {120., 60., 12., 1.};
-    M2.noalias() = M * M;
-    tmp = b[3]*M2 + b[1]*Id;
-    U.noalias() = M * tmp;
-    V = b[2]*M2 + b[0]*Id;
-  }
 
-  /** \internal \brief Compute the (5,5)-Pad&eacute; approximant to
+  public:
-   *  the exponential.
-   *  
-   *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
-   *  approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
-   *
-   *  \param M   Argument of matrix exponential
-   *  \param Id  Identity matrix of same size as M
-   *  \param tmp Temporary storage, to be provided by the caller
-   *  \param M2  Temporary storage, to be provided by the caller
-   *  \param U   Even-degree terms in numerator of Pad&eacute; approximant
-   *  \param V   Odd-degree terms in numerator of Pad&eacute; approximant
-   */
-  template <typename MatrixType>
-  EIGEN_STRONG_INLINE void pade5(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, 
-				 MatrixType& M2, MatrixType& U, MatrixType& V)
-  {
-    typedef typename ei_traits<MatrixType>::Scalar Scalar;
-    const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
-    M2.noalias() = M * M;
-    MatrixType M4 = M2 * M2;
-    tmp = b[5]*M4 + b[3]*M2 + b[1]*Id;
-    U.noalias() = M * tmp;
-    V = b[4]*M4 + b[2]*M2 + b[0]*Id;
-  }
   
-  /** \internal \brief Compute the (7,7)-Pad&eacute; approximant to
+    /** \brief Compute the matrix exponential. 
-   *  the exponential.
-   *  
-   *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
-   *  approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
-   *
-   *  \param M   Argument of matrix exponential
-   *  \param Id  Identity matrix of same size as M
-   *  \param tmp Temporary storage, to be provided by the caller
-   *  \param M2  Temporary storage, to be provided by the caller
-   *  \param U   Even-degree terms in numerator of Pad&eacute; approximant
-   *  \param V   Odd-degree terms in numerator of Pad&eacute; approximant
-   */
-  template <typename MatrixType>
-  EIGEN_STRONG_INLINE void pade7(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, 
-				 MatrixType& M2, MatrixType& U, MatrixType& V)
-  {
-    typedef typename ei_traits<MatrixType>::Scalar Scalar;
-    const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
-    M2.noalias() = M * M;
-    MatrixType M4 = M2 * M2;
-    MatrixType M6 = M4 * M2;
-    tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
-    U.noalias() = M * tmp;
-    V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
-  }
-  
-  /** \internal \brief Compute the (9,9)-Pad&eacute; approximant to
-   *  the exponential.
-   *  
-   *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
-   *  approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
-   *
-   *  \param M   Argument of matrix exponential
-   *  \param Id  Identity matrix of same size as M
-   *  \param tmp Temporary storage, to be provided by the caller
-   *  \param M2  Temporary storage, to be provided by the caller
-   *  \param U   Even-degree terms in numerator of Pad&eacute; approximant
-   *  \param V   Odd-degree terms in numerator of Pad&eacute; approximant
-   */
-  template <typename MatrixType>
-  EIGEN_STRONG_INLINE void pade9(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, 
-				 MatrixType& M2, MatrixType& U, MatrixType& V)
-  {
-    typedef typename ei_traits<MatrixType>::Scalar Scalar;
-    const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
-  		      2162160., 110880., 3960., 90., 1.};
-    M2.noalias() = M * M;
-    MatrixType M4 = M2 * M2;
-    MatrixType M6 = M4 * M2;
-    MatrixType M8 = M6 * M2;
-    tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
-    U.noalias() = M * tmp;
-    V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
-  }
-  
-  /** \internal \brief Compute the (13,13)-Pad&eacute; approximant to
-   *  the exponential.
-   *  
-   *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
-   *  approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$.
-   *
-   *  \param M   Argument of matrix exponential
-   *  \param Id  Identity matrix of same size as M
-   *  \param tmp Temporary storage, to be provided by the caller
-   *  \param M2  Temporary storage, to be provided by the caller
-   *  \param U   Even-degree terms in numerator of Pad&eacute; approximant
-   *  \param V   Odd-degree terms in numerator of Pad&eacute; approximant
-   */
-  template <typename MatrixType>
-  EIGEN_STRONG_INLINE void pade13(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, 
-				  MatrixType& M2, MatrixType& U, MatrixType& V)
-  {
-    typedef typename ei_traits<MatrixType>::Scalar Scalar;
-    const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., 
-  		      1187353796428800., 129060195264000., 10559470521600., 670442572800., 
-  		      33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
-    M2.noalias() = M * M;
-    MatrixType M4 = M2 * M2;
-    MatrixType M6 = M4 * M2;
-    V = b[13]*M6 + b[11]*M4 + b[9]*M2;
-    tmp.noalias() = M6 * V;
-    tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id;
-    U.noalias() = M * tmp;
-    tmp = b[12]*M6 + b[10]*M4 + b[8]*M2;
-    V.noalias() = M6 * tmp;
-    V += b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id;
-  }
-  
-  /** \internal \brief Helper class for computing Pad&eacute;
-   *  approximants to the exponential.
-   */
-  template <typename MatrixType, typename RealScalar = typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real>
-  struct computeUV_selector
-  {
-    /** \internal \brief Compute Pad&eacute; approximant to the exponential. 
-     *  
-     *  Computes \p U, \p V and \p squarings such that \f$ (V+U)(V-U)^{-1} \f$ 
-     *  is a Pad&eacute; of \f$ \exp(2^{-\mbox{squarings}}M) \f$
-     *  around \f$ M = 0 \f$. The degree of the Pad&eacute;
-     *  approximant and the value of squarings are chosen such that
-     *  the approximation error is no more than the round-off error.
-     *
-     *  \param M         Argument of matrix exponential
-     *  \param Id        Identity matrix of same size as M
-     *  \param tmp1      Temporary storage, to be provided by the caller
-     *  \param tmp2      Temporary storage, to be provided by the caller
-     *  \param U         Even-degree terms in numerator of Pad&eacute; approximant
-     *  \param V         Odd-degree terms in numerator of Pad&eacute; approximant
-     *  \param l1norm    L<sub>1</sub> norm of M
-     *  \param squarings Pointer to integer containing number of times
-     *                   that the result needs to be squared to find the
-     *                   matrix exponential 
-     */
-    static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, 
-		    MatrixType& U, MatrixType& V, float l1norm, int* squarings);
-  };
-  
-  template <typename MatrixType>
-  struct computeUV_selector<MatrixType, float>
-  {
-    static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, 
-		    MatrixType& U, MatrixType& V, float l1norm, int* squarings)
-    {
-      *squarings = 0;
-      if (l1norm < 4.258730016922831e-001) {
-        pade3(M, Id, tmp1, tmp2, U, V);
-      } else if (l1norm < 1.880152677804762e+000) {
-        pade5(M, Id, tmp1, tmp2, U, V);
-      } else {
-        const float maxnorm = 3.925724783138660f;
-        *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
-        MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings);
-        pade7(A, Id, tmp1, tmp2, U, V);
-      }
-    }
-  };
-  
-  template <typename MatrixType>
-  struct computeUV_selector<MatrixType, double>
-  {
-    static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, 
-		    MatrixType& U, MatrixType& V, float l1norm, int* squarings)
-    {
-      *squarings = 0;
-      if (l1norm < 1.495585217958292e-002) {
-        pade3(M, Id, tmp1, tmp2, U, V);
-      } else if (l1norm < 2.539398330063230e-001) {
-        pade5(M, Id, tmp1, tmp2, U, V);
-      } else if (l1norm < 9.504178996162932e-001) {
-        pade7(M, Id, tmp1, tmp2, U, V);
-      } else if (l1norm < 2.097847961257068e+000) {
-        pade9(M, Id, tmp1, tmp2, U, V);
-      } else {
-        const double maxnorm = 5.371920351148152;
-        *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm)));
-        MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings);
-        pade13(A, Id, tmp1, tmp2, U, V);
-      }
-    }
-  };
-  
-  /** \internal \brief Compute the matrix exponential. 
      *
      * \param M      matrix whose exponential is to be computed. 
      * \param result pointer to the matrix in which to store the result.
      */
+    MatrixExponential(const MatrixType &M, MatrixType *result);  
+
+  private:
+
+    // Prevent copying
+    MatrixExponential(const MatrixExponential&);
+    MatrixExponential& operator=(const MatrixExponential&);
+
+    /** \brief Compute the (3,3)-Pad&eacute; approximant to the exponential.
+     *  
+     *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
+     *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
+     *
+     *  \param A   Argument of matrix exponential
+     */
+    void pade3(const MatrixType &A);
+
+    /** \brief Compute the (5,5)-Pad&eacute; approximant to the exponential.
+     *  
+     *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
+     *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
+     *
+     *  \param A   Argument of matrix exponential
+     */
+    void pade5(const MatrixType &A);
+
+    /** \brief Compute the (7,7)-Pad&eacute; approximant to the exponential.
+     *  
+     *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
+     *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
+     *
+     *  \param A   Argument of matrix exponential
+     */
+    void pade7(const MatrixType &A);
+
+    /** \brief Compute the (9,9)-Pad&eacute; approximant to the exponential.
+     *  
+     *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
+     *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
+     *
+     *  \param A   Argument of matrix exponential
+     */
+    void pade9(const MatrixType &A);
+
+    /** \brief Compute the (13,13)-Pad&eacute; approximant to the exponential.
+     *  
+     *  After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Pad&eacute;
+     *  approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$.
+     *
+     *  \param A   Argument of matrix exponential
+     */
+    void pade13(const MatrixType &A);
+
+    /** \brief Compute Pad&eacute; approximant to the exponential. 
+     *  
+     * Computes \c m_U, \c m_V and \c m_squarings such that 
+     * \f$ (V+U)(V-U)^{-1} \f$ is a Pad&eacute; of 
+     * \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$. The
+     * degree of the Pad&eacute; approximant and the value of
+     * squarings are chosen such that the approximation error is no
+     * more than the round-off error.
+     *
+     * The argument of this function should correspond with the (real
+     * part of) the entries of \c m_M.  It is used to select the
+     * correct implementation using overloading.
+     */
+    void computeUV(double);
+
+    /** \brief Compute Pad&eacute; approximant to the exponential. 
+     *
+     *  \sa computeUV(double);
+     */
+    void computeUV(float);
+
+    typedef typename ei_traits<MatrixType>::Scalar Scalar;
+    typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar;
+
+    /** \brief Pointer to matrix whose exponential is to be computed. */
+    const MatrixType* m_M; 
+
+    /** \brief Even-degree terms in numerator of Pad&eacute; approximant. */
+    MatrixType m_U;
+
+    /** \brief Odd-degree terms in numerator of Pad&eacute; approximant. */
+    MatrixType m_V;
+
+    /** \brief Used for temporary storage. */
+    MatrixType m_tmp1;
+
+    /** \brief Used for temporary storage. */
+    MatrixType m_tmp2;
+
+    /** \brief Identity matrix of the same size as \c m_M. */
+    MatrixType m_Id;
+
+    /** \brief Number of squarings required in the last step. */
+    int m_squarings;
+
+    /** \brief L1 norm of m_M. */
+    float m_l1norm;
+};
+
 template <typename MatrixType>
-  void compute(const MatrixType &M, MatrixType* result)
+MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M, MatrixType *result) :
+  m_M(&M), 
+  m_U(M.rows(),M.cols()), 
+  m_V(M.rows(),M.cols()), 
+  m_tmp1(M.rows(),M.cols()), 
+  m_tmp2(M.rows(),M.cols()), 
+  m_Id(MatrixType::Identity(M.rows(), M.cols())), 
+  m_squarings(0), 
+  m_l1norm(static_cast<float>(M.cwise().abs().colwise().sum().maxCoeff()))
 {
-    MatrixType num(M.rows(), M.cols());
+  computeUV(RealScalar());
-    MatrixType den(M.rows(), M.cols());
+  m_tmp1 = m_U + m_V;	// numerator of Pade approximant
-    MatrixType U(M.rows(), M.cols());
+  m_tmp2 = -m_U + m_V;	// denominator of Pade approximant
-    MatrixType V(M.rows(), M.cols());
+  m_tmp2.partialLu().solve(m_tmp1, result);
-    MatrixType Id = MatrixType::Identity(M.rows(), M.cols());
+  for (int i=0; i<m_squarings; i++)
-    float l1norm = static_cast<float>(M.cwise().abs().colwise().sum().maxCoeff());
-    int squarings;
-    computeUV_selector<MatrixType>::run(M, Id, num, den, U, V, l1norm, &squarings);
-    num = U + V;			// numerator of Pade approximant
-    den = -U + V;			// denominator of Pade approximant
-    den.partialLu().solve(num, result);
-    for (int i=0; i<squarings; i++)
     *result *= *result;		// undo scaling by repeated squaring
 }
 
-} // end of namespace MatrixExponentialInternal
+template <typename MatrixType>
+EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A)
+{
+  const Scalar b[] = {120., 60., 12., 1.};
+  m_tmp1.noalias() = A * A;
+  m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id;
+  m_U.noalias() = A * m_tmp2;
+  m_V = b[2]*m_tmp1 + b[0]*m_Id;
+}
+
+template <typename MatrixType>
+EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A)
+{
+  const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.};
+  MatrixType A2 = A * A;
+  m_tmp1.noalias() = A2 * A2;
+  m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id;
+  m_U.noalias() = A * m_tmp2;
+  m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id;
+}
+
+template <typename MatrixType>
+EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A)
+{
+  const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.};
+  MatrixType A2 = A * A;
+  MatrixType A4 = A2 * A2;
+  m_tmp1.noalias() = A4 * A2;
+  m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
+  m_U.noalias() = A * m_tmp2;
+  m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
+}
+
+template <typename MatrixType>
+EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A)
+{
+  const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240.,
+  		      2162160., 110880., 3960., 90., 1.};
+  MatrixType A2 = A * A;
+  MatrixType A4 = A2 * A2;
+  MatrixType A6 = A4 * A2;
+  m_tmp1.noalias() = A6 * A2;
+  m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
+  m_U.noalias() = A * m_tmp2;
+  m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
+}
+
+template <typename MatrixType>
+EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A)
+{
+  const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., 
+  		      1187353796428800., 129060195264000., 10559470521600., 670442572800., 
+  		      33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.};
+  MatrixType A2 = A * A;
+  MatrixType A4 = A2 * A2;
+  m_tmp1.noalias() = A4 * A2;
+  m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage
+  m_tmp2.noalias() = m_tmp1 * m_V;
+  m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id;
+  m_U.noalias() = A * m_tmp2;
+  m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2;
+  m_V.noalias() = m_tmp1 * m_tmp2;
+  m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id;
+}
+
+template <typename MatrixType>
+void MatrixExponential<MatrixType>::computeUV(float)
+{
+  if (m_l1norm < 4.258730016922831e-001) {
+    pade3(*m_M);
+  } else if (m_l1norm < 1.880152677804762e+000) {
+    pade5(*m_M);
+  } else {
+    const float maxnorm = 3.925724783138660f;
+    m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    MatrixType A = *m_M / std::pow(Scalar(2), m_squarings);
+    pade7(A);
+  }
+}
+
+template <typename MatrixType>
+void MatrixExponential<MatrixType>::computeUV(double)
+{
+  if (m_l1norm < 1.495585217958292e-002) {
+    pade3(*m_M);
+  } else if (m_l1norm < 2.539398330063230e-001) {
+    pade5(*m_M);
+  } else if (m_l1norm < 9.504178996162932e-001) {
+    pade7(*m_M);
+  } else if (m_l1norm < 2.097847961257068e+000) {
+    pade9(*m_M);
+  } else {
+    const double maxnorm = 5.371920351148152;
+    m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    MatrixType A = *m_M / std::pow(Scalar(2), m_squarings);
+    pade13(A);
+  }
+}
 
 template <typename Derived>
 EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M, 
 					       typename MatrixBase<Derived>::PlainMatrixType* result)
 {
   ei_assert(M.rows() == M.cols());
-  MatrixExponentialInternal::compute(M.eval(), result);
+  MatrixExponential<typename MatrixBase<Derived>::PlainMatrixType>(M, result);
 }
 
 #endif // EIGEN_MATRIX_EXPONENTIAL