* split Product to a DiagonalProduct template specialization

to optimize matrix-diag and diag-matrix products without
  making Product over complicated.
* compilation fixes in Tridiagonalization and HessenbergDecomposition
  in the case of 2x2 matrices.
* added an Orientation2D small class with similar interface than Quaternion
  (used by Transform to handle 2D and 3D orientations seamlessly)
* added a couple of features in Transform.

Eigen/Core

@@ -41,6 +41,7 @@ namespace Eigen {
 #include "src/Core/CwiseNullaryOp.h"
 #include "src/Core/InverseProduct.h"
 #include "src/Core/Product.h"
+#include "src/Core/DiagonalProduct.h"
 #include "src/Core/Block.h"
 #include "src/Core/Minor.h"
 #include "src/Core/Transpose.h"

Eigen/src/Core/DiagonalProduct.h

@@ -0,0 +1,108 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_DIAGONALPRODUCT_H
+#define EIGEN_DIAGONALPRODUCT_H
+
+template<typename Lhs, typename Rhs>
+struct ei_traits<Product<Lhs, Rhs, DiagonalProduct> >
+{
+  typedef typename Lhs::Scalar Scalar;
+  typedef typename ei_nested<Lhs>::type LhsNested;
+  typedef typename ei_nested<Rhs>::type RhsNested;
+  typedef typename ei_unref<LhsNested>::type _LhsNested;
+  typedef typename ei_unref<RhsNested>::type _RhsNested;
+  enum {
+    LhsFlags = _LhsNested::Flags,
+    RhsFlags = _RhsNested::Flags,
+    RowsAtCompileTime = Lhs::RowsAtCompileTime,
+    ColsAtCompileTime = Rhs::ColsAtCompileTime,
+    MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = Rhs::MaxColsAtCompileTime,
+    _RhsVectorizable =  (RhsFlags & RowMajorBit) && (RhsFlags & VectorizableBit)
+                     && (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
+    _LhsVectorizable =  (!(LhsFlags & RowMajorBit)) && (LhsFlags & VectorizableBit)
+                     && (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
+    _LostBits = ~(((RhsFlags & RowMajorBit) && (!_LhsVectorizable) ? 0 : RowMajorBit)
+                | ((RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic) ? 0 : LargeBit)),
+    Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & _LostBits)
+          | (_LhsVectorizable || _RhsVectorizable ? VectorizableBit : 0),
+    CoeffReadCost = NumTraits<Scalar>::MulCost + _LhsNested::CoeffReadCost + _RhsNested::CoeffReadCost
+  };
+};
+
+template<typename Lhs, typename Rhs> class Product<Lhs, Rhs, DiagonalProduct> : ei_no_assignment_operator,
+  public MatrixBase<Product<Lhs, Rhs, DiagonalProduct> >
+{
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
+    typedef typename ei_traits<Product>::LhsNested LhsNested;
+    typedef typename ei_traits<Product>::RhsNested RhsNested;
+    typedef typename ei_traits<Product>::_LhsNested _LhsNested;
+    typedef typename ei_traits<Product>::_RhsNested _RhsNested;
+
+    enum {
+      PacketSize = ei_packet_traits<Scalar>::size
+    };
+
+    inline Product(const Lhs& lhs, const Rhs& rhs)
+      : m_lhs(lhs), m_rhs(rhs)
+    {
+      ei_assert(lhs.cols() == rhs.rows());
+    }
+
+  private:
+
+    inline int _rows() const { return m_lhs.rows(); }
+    inline int _cols() const { return m_rhs.cols(); }
+
+    const Scalar _coeff(int row, int col) const
+    {
+      int unique = ((Rhs::Flags&Diagonal)==Diagonal) ? col : row;
+      return m_lhs.coeff(row, unique) * m_rhs.coeff(unique, col);
+    }
+
+    template<int LoadMode>
+    const PacketScalar _packetCoeff(int row, int col) const
+    {
+      if ((Rhs::Flags&Diagonal)==Diagonal)
+      {
+        ei_assert((_LhsNested::Flags&RowMajorBit)==0);
+        return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
+      }
+      else
+      {
+        ei_assert(_RhsNested::Flags&RowMajorBit);
+        return ei_pmul(ei_pset1(m_lhs.coeff(row, row)), m_rhs.template packetCoeff<LoadMode>(row, col));
+      }
+    }
+
+  protected:
+    const LhsNested m_lhs;
+    const RhsNested m_rhs;
+};
+
+#endif // EIGEN_DIAGONALPRODUCT_H

Eigen/src/Core/Product.h

@@ -144,11 +144,12 @@ struct ei_packet_product_impl<false, Index, Dynamic, Lhs, Rhs, PacketScalar>
   */
 template<typename Lhs, typename Rhs> struct ei_product_eval_mode
 {
-  enum{ value =  Lhs::MaxRowsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-              && Rhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-              && Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
-              && (Rhs::Flags&Diagonal)!=Diagonal
-              ? CacheFriendlyProduct : NormalProduct };
+  enum{ value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal)
+              ? DiagonalProduct
+              :    Lhs::MaxRowsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+                && Rhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+                && Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+                ? CacheFriendlyProduct : NormalProduct };
 };
 
 template<typename T> class ei_product_eval_to_column_major
@@ -265,44 +266,25 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
 
     const Scalar _coeff(int row, int col) const
     {
-      if ((Rhs::Flags&Diagonal)==Diagonal)
-      {
-        return  m_lhs.coeff(row, col) * m_rhs.coeff(col, col);
-      }
-      else if ((Lhs::Flags&Diagonal)==Diagonal)
-      {
-        return  m_lhs.coeff(row, row) * m_rhs.coeff(row, col);
-      }
-      else
-      {
-        Scalar res;
-        const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
-        ei_product_impl<unroll ? Lhs::ColsAtCompileTime-1 : Dynamic,
-                            unroll ? Lhs::ColsAtCompileTime : Dynamic,
-                            _LhsNested, _RhsNested>
-          ::run(row, col, m_lhs, m_rhs, res);
-        return res;
-      }
+      Scalar res;
+      const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
+      ei_product_impl<unroll ? Lhs::ColsAtCompileTime-1 : Dynamic,
+                          unroll ? Lhs::ColsAtCompileTime : Dynamic,
+                          _LhsNested, _RhsNested>
+        ::run(row, col, m_lhs, m_rhs, res);
+      return res;
     }
 
     template<int LoadMode>
     const PacketScalar _packetCoeff(int row, int col) const
     {
-      if ((Rhs::Flags&Diagonal)==Diagonal)
-      {
-        ei_assert((_LhsNested::Flags&RowMajorBit)==0);
-        return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
-      }
-      else
-      {
-        const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
-        PacketScalar res;
-        ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1,
-                            unroll ? Lhs::ColsAtCompileTime : Dynamic,
-                            _LhsNested, _RhsNested, PacketScalar>
-          ::run(row, col, m_lhs, m_rhs, res);
-        return res;
-      }
+      const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
+      PacketScalar res;
+      ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1,
+                          unroll ? Lhs::ColsAtCompileTime : Dynamic,
+                          _LhsNested, _RhsNested, PacketScalar>
+        ::run(row, col, m_lhs, m_rhs, res);
+      return res;
     }
 
     template<typename Lhs_, typename Rhs_, int EvalMode_, typename DestDerived_, bool DirectAccess_>

Eigen/src/Core/Transpose.h

@@ -52,8 +52,6 @@ struct ei_traits<Transpose<MatrixType> >
           & ~( Like1DArrayBit | LowerTriangularBit | UpperTriangularBit)
           | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
           | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
-    // Note that the above test cannot be simplified using ^ because a diagonal matrix
-    // has both LowerTriangularBit and UpperTriangularBit and both must be preserved.
     CoeffReadCost = _MatrixTypeNested::CoeffReadCost
   };
 };

Eigen/src/Core/util/Constants.h

@@ -80,7 +80,7 @@ const unsigned int Like1DArrayBit = 0x20;
 /** \ingroup flags
   *
   * means all diagonal coefficients are equal to 0 */
-const unsigned int ZeroDiagBit = 0x40;      
+const unsigned int ZeroDiagBit = 0x40;
 
 /** \ingroup flags
   *
@@ -140,7 +140,7 @@ enum { Aligned=0, UnAligned=1 };
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 enum DirectionType { Vertical, Horizontal };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, LazyProduct};
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, LazyProduct};
 
 
 #endif // EIGEN_CONSTANTS_H

Eigen/src/Geometry/Transform.h

@@ -25,6 +25,94 @@
 #ifndef EIGEN_TRANSFORM_H
 #define EIGEN_TRANSFORM_H
 
+/** \class Orientation2D
+  *
+  * \brief Represents an orientation/rotation in a 2 dimensional space.
+  *
+  * \param _Scalar the scalar type, i.e., the type of the coefficients
+  *
+  * This class is equivalent to a single scalar representating the rotation angle
+  * in radian with some additional features such as the conversion from/to
+  * rotation matrix. Moreover this class aims to provide a similar interface
+  * to Quaternion in order to facilitate the writting of generic algorithm
+  * dealing with rotations.
+  *
+  * \sa class Quaternion, class Transform
+  */
+template<typename _Scalar>
+class Orientation2D
+{
+public:
+  enum { Dim = 2 };
+  /** the scalar type of the coefficients */
+  typedef _Scalar Scalar;
+  typedef Matrix<Scalar,2,2> Matrix2;
+
+protected:
+
+  Scalar m_angle;
+
+public:
+
+  inline Orientation2D(Scalar a) : m_angle(a) {}
+  inline operator Scalar& () { return m_angle; }
+  inline operator Scalar () const { return m_angle; }
+
+  template<typename Derived>
+  Orientation2D& fromRotationMatrix(const MatrixBase<Derived>& m);
+  Matrix2 toRotationMatrix(void) const;
+
+  Orientation2D slerp(Scalar t, const Orientation2D& other) const;
+};
+
+/** returns the default type used to represent an orientation.
+  */
+template<typename Scalar, int Dim>
+struct ei_get_orientation_type;
+
+template<typename Scalar>
+struct ei_get_orientation_type<Scalar,2>
+{ typedef Orientation2D<Scalar> type; };
+
+template<typename Scalar>
+struct ei_get_orientation_type<Scalar,3>
+{ typedef Quaternion<Scalar> type; };
+
+/** Set \c *this from a 2x2 rotation matrix \a mat.
+  * In other words, this function extract the rotation angle
+  * from the rotation matrix.
+  */
+template<typename Scalar>
+template<typename Derived>
+Orientation2D<Scalar>& Orientation2D<Scalar>::fromRotationMatrix(const MatrixBase<Derived>& mat)
+{
+  EIGEN_STATIC_ASSERT(Derived::RowsAtCompileTime==2 && Derived::ColsAtCompileTime==2,you_did_a_programming_error);
+  m_angle = ei_atan2(mat.coeff(1,0), mat.coeff(0,0));
+  return *this;
+}
+
+/** Constructs and \returns an equivalent 2x2 rotation matrix.
+  */
+template<typename Scalar>
+typename Orientation2D<Scalar>::Matrix2
+Orientation2D<Scalar>::toRotationMatrix(void) const
+{
+  Scalar sinA = ei_sin(m_angle);
+  Scalar cosA = ei_cos(m_angle);
+  return Matrix2(cosA, -sinA, sinA, cosA);
+}
+
+/** \returns the spherical interpolation between \c *this and \a other using
+  * parameter \a t. It is equivalent to a linear interpolation.
+  */
+template<typename Scalar>
+Orientation2D<Scalar>
+Orientation2D<Scalar>::slerp(Scalar t, const Orientation2D& other) const
+{
+  return m_angle * (1-t) + t * other;
+}
+
+
 /** \class Transform
   *
   * \brief Represents an homogeneous transformation in a N dimensional space
@@ -32,6 +120,11 @@
   * \param _Scalar the scalar type, i.e., the type of the coefficients
   * \param _Dim the dimension of the space
   *
+  * The homography is internally represented and stored as a (Dim+1)^2 matrix which
+  * is available through the matrix() method.
+  *
+  * Conversion methods from/to Qt's QMatrix are available if the preprocessor token
+  * EIGEN_QT_SUPPORT is defined.
   *
   */
 template<typename _Scalar, int _Dim>
@@ -47,6 +140,7 @@ public:
   typedef Block<MatrixType,Dim,Dim> AffineMatrixRef;
   typedef Matrix<Scalar,Dim,1> VectorType;
   typedef Block<MatrixType,Dim,1> VectorRef;
+  typedef typename ei_get_orientation_type<Scalar,Dim>::type Orientation;
 
 protected:
 
@@ -59,13 +153,42 @@ protected:
 
 public:
 
+  /** Default constructor without initialization of the coefficients. */
+  Transform() { }
+
+  inline Transform(const Transform& other)
+  { m_matrix = other.m_matrix; }
+
+  inline Transform& operator=(const Transform& other)
+  { m_matrix = other.m_matrix; }
+
+  template<typename OtherDerived>
+  inline explicit Transform(const MatrixBase<OtherDerived>& other)
+  { m_matrix = other; }
+
+  template<typename OtherDerived>
+  inline Transform& operator=(const MatrixBase<OtherDerived>& other)
+  { m_matrix = other; }
+
+  #ifdef EIGEN_QT_SUPPORT
+  inline Transform(const QMatrix& other);
+  inline Transform& operator=(const QMatrix& other);
+  inline QMatrix toQMatrix(void) const;
+  #endif
+
+  /** \returns a read-only expression of the transformation matrix */
   inline const MatrixType matrix() const { return m_matrix; }
+  /** \returns a writable expression of the transformation matrix */
   inline MatrixType matrix() { return m_matrix; }
 
+  /** \returns a read-only expression of the affine (linear) part of the transformation */
   inline const AffineMatrixRef affine() const { return m_matrix.template block<Dim,Dim>(0,0); }
+  /** \returns a writable expression of the affine (linear) part of the transformation */
   inline AffineMatrixRef affine() { return m_matrix.template block<Dim,Dim>(0,0); }
 
+  /** \returns a read-only expression of the translation vector of the transformation */
   inline const VectorRef translation() const { return m_matrix.template block<Dim,1>(0,Dim); }
+  /** \returns a writable expression of the translation vector of the transformation */
   inline VectorRef translation() { return m_matrix.template block<Dim,1>(0,Dim); }
 
   template<typename OtherDerived>
@@ -78,6 +201,11 @@ public:
   const typename ProductReturnType<OtherDerived>::Type
   operator * (const MatrixBase<OtherDerived> &other) const;
 
+  /** Contatenates two transformations */
+  Product<MatrixType,MatrixType>
+  operator * (const Transform& other) const
+  { return m_matrix * other.matrix(); }
+
   void setIdentity() { m_matrix.setIdentity(); }
 
   template<typename OtherDerived>
@@ -92,13 +220,58 @@ public:
   template<typename OtherDerived>
   Transform& pretranslate(const MatrixBase<OtherDerived> &other);
 
+  template<typename OtherDerived>
+  Transform& shear(Scalar sx, Scalar sy);
+
+  template<typename OtherDerived>
+  Transform& preshear(Scalar sx, Scalar sy);
+
   AffineMatrixType extractRotation() const;
   AffineMatrixType extractRotationNoShear() const;
 
+  template<typename PositionDerived, typename ScaleDerived>
+  Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position,
+    const Orientation& orientation, const MatrixBase<ScaleDerived> &scale);
+
+  const Inverse<MatrixType, false> inverse() const
+  { return m_matrix.inverse(); }
+
 protected:
 
 };
 
+#ifdef EIGEN_QT_SUPPORT
+/** Initialises \c *this from a QMatrix assuming the dimension is 2.
+  */
+template<typename Scalar, int Dim>
+Transform<Scalar,Dim>::Transform(const QMatrix& other)
+{
+  *this = other;
+}
+
+/** Set \c *this from a QMatrix assuming the dimension is 2.
+  */
+template<typename Scalar, int Dim>
+Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QMatrix& other)
+{
+  EIGEN_STATIC_ASSERT(Dim==2, you_did_a_programming_error);
+  m_matrix << other.m11(), other.m21(), other.dx(),
+              other.m12(), other.m22(), other.dy(),
+              0, 0, 1;
+}
+
+/** \returns a QMatrix from \c *this assuming the dimension is 2.
+  */
+template<typename Scalar, int Dim>
+QMatrix Transform<Scalar,Dim>::toQMatrix(void) const
+{
+  EIGEN_STATIC_ASSERT(Dim==2, you_did_a_programming_error);
+  return QMatrix( other.coeffRef(0,0), other.coeffRef(1,0),
+                  other.coeffRef(0,1), other.coeffRef(1,1),
+                  other.coeffRef(0,2), other.coeffRef(1,2),
+}
+#endif
+
 template<typename Scalar, int Dim>
 template<typename OtherDerived>
 const typename Transform<Scalar,Dim>::template ProductReturnType<OtherDerived>::Type
@@ -133,7 +306,7 @@ Transform<Scalar,Dim>::prescale(const MatrixBase<OtherDerived> &other)
 {
   EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime)
     && int(OtherDerived::SizeAtCompileTime)==int(Dim), you_did_a_programming_error);
-  m_matrix.template block<3,4>(0,0) = (other.asDiagonal().eval() * m_matrix.template block<3,4>(0,0)).lazy();
+  m_matrix.template block<3,4>(0,0) = (other.asDiagonal() * m_matrix.template block<3,4>(0,0)).lazy();
   return *this;
 }
 
@@ -167,6 +340,39 @@ Transform<Scalar,Dim>::pretranslate(const MatrixBase<OtherDerived> &other)
   return *this;
 }
 
+/** Applies on the right the shear transformation represented
+  * by the vector \a other to \c *this and returns a reference to \c *this.
+  * \warning 2D only.
+  * \sa preshear()
+  */
+template<typename Scalar, int Dim>
+template<typename OtherDerived>
+Transform<Scalar,Dim>&
+Transform<Scalar,Dim>::shear(Scalar sx, Scalar sy)
+{
+  EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime)
+    && int(OtherDerived::SizeAtCompileTime)==int(Dim) && int(Dim)==2, you_did_a_programming_error);
+  VectorType tmp = affine().col(0)*sy + affine().col(1);
+  affine() << affine().col(0) + affine().col(1)*sx, tmp;
+  return *this;
+}
+
+/** Applies on the left the shear transformation represented
+  * by the vector \a other to \c *this and returns a reference to \c *this.
+  * \warning 2D only.
+  * \sa shear()
+  */
+template<typename Scalar, int Dim>
+template<typename OtherDerived>
+Transform<Scalar,Dim>&
+Transform<Scalar,Dim>::preshear(Scalar sx, Scalar sy)
+{
+  EIGEN_STATIC_ASSERT(int(OtherDerived::IsVectorAtCompileTime)
+    && int(OtherDerived::SizeAtCompileTime)==int(Dim), you_did_a_programming_error);
+  m_matrix.template block<3,4>(0,0) = AffineMatrixType(1, sx, sy, 1) * m_matrix.template block<3,4>(0,0);
+  return *this;
+}
+
 /** \returns the rotation part of the transformation using a QR decomposition.
   * \sa extractRotationNoShear()
   */
@@ -189,6 +395,22 @@ Transform<Scalar,Dim>::extractRotationNoShear() const
             .verticalRedux(ei_scalar_sum_op<Scalar>()).cwiseSqrt();
 }
 
+/** Convenient method to set \c *this from a position, orientation and scale
+  * of a 3D object.
+  */
+template<typename Scalar, int Dim>
+template<typename PositionDerived, typename ScaleDerived>
+Transform<Scalar,Dim>&
+Transform<Scalar,Dim>::fromPositionOrientationScale(const MatrixBase<PositionDerived> &position,
+  const Orientation& orientation, const MatrixBase<ScaleDerived> &scale)
+{
+  affine() = orientation.toRotationMatrix();
+  translation() = position;
+  m_matrix(Dim,Dim) = 1.;
+  m_matrix.template block<1,Dim>(Dim,0).setZero();
+  affine() *= scale.asDiagonal();
+}
+
 //----------
 
 template<typename Scalar, int Dim>
@@ -216,12 +438,17 @@ template<typename Other>
 struct Transform<Scalar,Dim>::ei_transform_product_impl<Other,Dim,1>
 {
   typedef typename Transform<Scalar,Dim>::AffineMatrixRef MatrixType;
-  typedef const CwiseBinaryOp<
-    ei_scalar_sum_op<Scalar>,
-    NestByValue<Product<NestByValue<MatrixType>,Other> >,
-    NestByValue<typename Transform<Scalar,Dim>::VectorRef> > ResultType;
+  typedef const CwiseUnaryOp<
+      ei_scalar_multiple_op<Scalar>,
+      NestByValue<CwiseBinaryOp<
+        ei_scalar_sum_op<Scalar>,
+        NestByValue<Product<NestByValue<MatrixType>,Other> >,
+        NestByValue<typename Transform<Scalar,Dim>::VectorRef> > >
+      > ResultType;
+  // FIXME shall we offer an optimized version when the last row is know to be 0,0...,0,1 ?
   static ResultType run(const Transform<Scalar,Dim>& tr, const Other& other)
-  { return (tr.affine().nestByValue() * other).nestByValue() + tr.translation().nestByValue(); }
+  { return ((tr.affine().nestByValue() * other).nestByValue() + tr.translation().nestByValue()).nestByValue()
+          * (Scalar(1) / ( (tr.matrix().template block<1,Dim>(Dim,0) * other).coeff(0) + tr.matrix().coeff(Dim,Dim))); }
 };
 
 #endif // EIGEN_TRANSFORM_H

Eigen/src/QR/HessenbergDecomposition.h

@@ -239,7 +239,8 @@ HessenbergDecomposition<MatrixType>::matrixH(void) const
   // and fill it (to avoid temporaries)
   int n = m_matrix.rows();
   MatrixType matH = m_matrix;
-  matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
+  if (n>2)
+    matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
   return matH;
 }
 

Eigen/src/QR/Tridiagonalization.h

@@ -160,8 +160,11 @@ Tridiagonalization<MatrixType>::matrixT(void) const
   int n = m_matrix.rows();
   MatrixType matT = m_matrix;
   matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate();
-  matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero();
-  matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
+  if (n>2)
+  {
+    matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero();
+    matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
+  }
   return matT;
 }
 

test/geometry.cpp

@@ -116,6 +116,9 @@ template<typename Scalar> void geometry(void)
   t1.translate(-v0);
 
   VERIFY((t0.matrix() * t1.matrix()).isIdentity());
+
+  t1.fromPositionOrientationScale(v0, q1, v1);
+  VERIFY_IS_APPROX(t1.matrix(), t0.matrix());
 }
 
 void test_geometry()