diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h
index 71549c43e..69e906ee7 100644
--- a/Eigen/src/Cholesky/Cholesky.h
+++ b/Eigen/src/Cholesky/Cholesky.h
@@ -58,7 +58,7 @@ template<typename MatrixType> class Cholesky
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
     template<typename Derived>
-    typename Derived::Eval solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
+    typename MatrixBase<Derived>::PlainMatrixType_ColMajor solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
 
     template<typename RhsDerived, typename ResDerived>
     bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
@@ -119,11 +119,11 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
 /** \deprecated */
 template<typename MatrixType>
 template<typename Derived>
-typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b) const
+typename MatrixBase<Derived>::PlainMatrixType_ColMajor Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows());
-  typename ei_eval_to_column_major<Derived>::type x(b);
+  typename MatrixBase<Derived>::PlainMatrixType_ColMajor x(b);
   solveInPlace(x);
   return x;
 }
@@ -156,10 +156,10 @@ bool Cholesky<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
   * \deprecated has been renamed llt()
   */
 template<typename Derived>
-inline const Cholesky<typename MatrixBase<Derived>::EvalType>
+inline const Cholesky<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::cholesky() const
 {
-  return Cholesky<typename ei_eval<Derived>::type>(derived());
+  return Cholesky<PlainMatrixType>(derived());
 }
 
 #endif // EIGEN_CHOLESKY_H
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index bf9951709..b40ba18c0 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -175,7 +175,7 @@ bool CholeskyWithoutSquareRoot<MatrixType>::solveInPlace(MatrixBase<Derived> &bA
   * \deprecated has been renamed ldlt()
   */
 template<typename Derived>
-inline const CholeskyWithoutSquareRoot<typename MatrixBase<Derived>::EvalType>
+inline const CholeskyWithoutSquareRoot<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::choleskyNoSqrt() const
 {
   return derived();
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index aa967f6b9..a275e093f 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -189,7 +189,7 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
   * \returns the Cholesky decomposition without square root of \c *this
   */
 template<typename Derived>
-inline const LDLT<typename MatrixBase<Derived>::EvalType>
+inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::ldlt() const
 {
   return derived();
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 2fc658bb7..6e412e20c 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -177,10 +177,10 @@ bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
   * \returns the LLT decomposition of \c *this
   */
 template<typename Derived>
-inline const LLT<typename MatrixBase<Derived>::EvalType>
+inline const LLT<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::llt() const
 {
-  return LLT<typename ei_eval<Derived>::type>(derived());
+  return LLT<PlainMatrixType>(derived());
 }
 
 #endif // EIGEN_LLT_H
diff --git a/Eigen/src/Core/DiagonalProduct.h b/Eigen/src/Core/DiagonalProduct.h
index ca0b56872..5e23fb066 100644
--- a/Eigen/src/Core/DiagonalProduct.h
+++ b/Eigen/src/Core/DiagonalProduct.h
@@ -32,7 +32,7 @@
  */
 template<typename T, int N> struct ei_nested_diagonal : ei_nested<T,N> {};
 template<typename T, int N> struct ei_nested_diagonal<DiagonalMatrix<T>,N >
- : ei_nested<DiagonalMatrix<T>, N, DiagonalMatrix<NestByValue<typename ei_eval<T>::type> > >
+ : ei_nested<DiagonalMatrix<T>, N, DiagonalMatrix<NestByValue<typename ei_plain_matrix_type<T>::type> > >
 {};
 
 // specialization of ProductReturnType
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index c4703adc3..86bebe246 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -318,7 +318,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<
   * \sa norm(), normalize()
   */
 template<typename Derived>
-inline const typename MatrixBase<Derived>::EvalType
+inline const typename MatrixBase<Derived>::PlainMatrixType
 MatrixBase<Derived>::normalized() const
 {
   typedef typename ei_nested<Derived>::type Nested;
diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h
index ca00cae3d..2b00d5bc5 100644
--- a/Eigen/src/Core/IO.h
+++ b/Eigen/src/Core/IO.h
@@ -122,9 +122,10 @@ MatrixBase<Derived>::format(const IOFormat& fmt) const
 /** \internal
   * print the matrix \a _m to the output stream \a s using the output format \a fmt */
 template<typename Derived>
-std::ostream & ei_print_matrix(std::ostream & s, const MatrixBase<Derived> & _m, const IOFormat& fmt)
+std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOFormat& fmt)
 {
   const typename Derived::Nested m = _m;
+
   int width = 0;
   if (fmt.flags & AlignCols)
   {
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index fd4b5cb4a..3cde1e28b 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -426,14 +426,6 @@ class Matrix
     /** Destructor */
     inline ~Matrix() {}
 
-    /** Override MatrixBase::eval() since matrices don't need to be evaluated, it is enough to just read them.
-      * This prevents a useless copy when doing e.g. "m1 = m2.eval()"
-      */
-    inline const Matrix& eval() const
-    {
-      return *this;
-    }
-
     /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
       * data pointers.
       */
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f916dbf2a..bb3cc0532 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -179,9 +179,20 @@ template<typename Derived> class MatrixBase
     int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
-    /** \internal the type to which the expression gets evaluated (needed by MSVC) */
-    typedef typename ei_eval<Derived>::type EvalType;
-    /** \internal Represents a constant matrix */
+    /** \internal the plain matrix type corresponding to this expression. Note that is not necessarily
+      * exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const
+      * reference to a matrix, not a matrix! It guaranteed however, that the return type of eval() is either
+      * PlainMatrixType or const PlainMatrixType&.
+      */
+    typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
+    /** \internal the column-major plain matrix type corresponding to this expression. Note that is not necessarily
+      * exactly the return type of eval(): in the case of plain matrices, the return type of eval() is a const
+      * reference to a matrix, not a matrix!
+      * The only difference from PlainMatrixType is that PlainMatrixType_ColMajor is guaranteed to be column-major.
+      */
+    typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType_ColMajor;
+
+    /** \internal Represents a matrix with all coefficients equal to one another*/
     typedef CwiseNullaryOp<ei_scalar_constant_op<Scalar>,Derived> ConstantReturnType;
     /** \internal Represents a scalar multiple of a matrix */
     typedef CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, Derived> ScalarMultipleReturnType;
@@ -331,7 +342,7 @@ template<typename Derived> class MatrixBase
     Derived& operator*=(const MatrixBase<OtherDerived>& other);
 
     template<typename OtherDerived>
-    typename ei_eval_to_column_major<OtherDerived>::type
+    typename ei_plain_matrix_type_column_major<OtherDerived>::type
 		solveTriangular(const MatrixBase<OtherDerived>& other) const;
 
     template<typename OtherDerived>
@@ -343,7 +354,7 @@ template<typename Derived> class MatrixBase
     RealScalar squaredNorm() const;
     RealScalar norm2() const;
     RealScalar norm()  const;
-    const EvalType normalized() const;
+    const PlainMatrixType normalized() const;
     void normalize();
 
     Eigen::Transpose<Derived> transpose();
@@ -481,6 +492,8 @@ template<typename Derived> class MatrixBase
 
     /** \returns the matrix or vector obtained by evaluating this expression.
       *
+      * Notice that in the case of a plain matrix or vector (not an expression) this function just returns
+      * a const reference, in order to avoid a useless copy.
       */
     EIGEN_ALWAYS_INLINE const typename ei_eval<Derived>::type eval() const
     {
@@ -573,35 +586,35 @@ template<typename Derived> class MatrixBase
 
 /////////// LU module ///////////
 
-    const LU<EvalType> lu() const;
-    const EvalType inverse() const;
-    void computeInverse(EvalType *result) const;
+    const LU<PlainMatrixType> lu() const;
+    const PlainMatrixType inverse() const;
+    void computeInverse(PlainMatrixType *result) const;
     Scalar determinant() const;
 
 /////////// Cholesky module ///////////
 
-    const LLT<EvalType>  llt() const;
-    const LDLT<EvalType> ldlt() const;
+    const LLT<PlainMatrixType>  llt() const;
+    const LDLT<PlainMatrixType> ldlt() const;
     // deprecated:
-    const Cholesky<EvalType> cholesky() const;
-    const CholeskyWithoutSquareRoot<EvalType> choleskyNoSqrt() const;
+    const Cholesky<PlainMatrixType> cholesky() const;
+    const CholeskyWithoutSquareRoot<PlainMatrixType> choleskyNoSqrt() const;
 
 /////////// QR module ///////////
 
-    const QR<EvalType> qr() const;
+    const QR<PlainMatrixType> qr() const;
 
     EigenvaluesReturnType eigenvalues() const;
     RealScalar operatorNorm() const;
 
 /////////// SVD module ///////////
 
-    SVD<EvalType> svd() const;
+    SVD<PlainMatrixType> svd() const;
 
 /////////// Geometry module ///////////
 
     template<typename OtherDerived>
-    EvalType cross(const MatrixBase<OtherDerived>& other) const;
-    EvalType unitOrthogonal(void) const;
+    PlainMatrixType cross(const MatrixBase<OtherDerived>& other) const;
+    PlainMatrixType unitOrthogonal(void) const;
     Matrix<Scalar,3,1> eulerAngles(int a0, int a1, int a2) const;
 
     #ifdef EIGEN_MATRIXBASE_PLUGIN
diff --git a/Eigen/src/Core/Part.h b/Eigen/src/Core/Part.h
index 3cb55fe1d..e51eaeb29 100644
--- a/Eigen/src/Core/Part.h
+++ b/Eigen/src/Core/Part.h
@@ -157,7 +157,7 @@ inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator=(const Other& ot
 {
   if(Other::Flags & EvalBeforeAssigningBit)
   {
-    typename ei_eval<Other>::type other_evaluated(other.rows(), other.cols());
+    typename MatrixBase<Other>::PlainMatrixType other_evaluated(other.rows(), other.cols());
     other_evaluated.template part<Mode>().lazyAssign(other);
     lazyAssign(other_evaluated);
   }
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index a844470a7..4e5bea050 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -69,7 +69,7 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
   typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
 
   typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime,
-                             typename ei_eval_to_column_major<Rhs>::type
+                             typename ei_plain_matrix_type_column_major<Rhs>::type
                    >::type RhsNested;
 
   typedef Product<LhsNested, RhsNested, CacheFriendlyProduct> Type;
@@ -735,7 +735,7 @@ template<typename T> struct ei_product_copy_rhs
   typedef typename ei_meta_if<
          (ei_traits<T>::Flags & RowMajorBit)
       || (!(ei_traits<T>::Flags & DirectAccessBit)),
-      typename ei_eval_to_column_major<T>::type,
+      typename ei_plain_matrix_type_column_major<T>::type,
       const T&
     >::ret type;
 };
@@ -744,7 +744,7 @@ template<typename T> struct ei_product_copy_lhs
 {
   typedef typename ei_meta_if<
       (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
-      typename ei_eval<T>::type,
+      typename ei_plain_matrix_type<T>::type,
       const T&
     >::ret type;
 };
diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h
index 20c0408bd..b58dab01d 100644
--- a/Eigen/src/Core/SolveTriangular.h
+++ b/Eigen/src/Core/SolveTriangular.h
@@ -236,7 +236,7 @@ void MatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other
   enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit };
 
   typedef typename ei_meta_if<copy,
-    typename ei_eval_to_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy;
+    typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy;
   OtherCopy otherCopy(other.derived());
 
   ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy);
@@ -278,10 +278,10 @@ void MatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other
   */
 template<typename Derived>
 template<typename OtherDerived>
-typename ei_eval_to_column_major<OtherDerived>::type
+typename ei_plain_matrix_type_column_major<OtherDerived>::type
 MatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
 {
-  typename ei_eval_to_column_major<OtherDerived>::type res(other);
+  typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other);
   solveTriangularInPlace(res);
   return res;
 }
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 312d863e7..dc18a425c 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -171,7 +171,6 @@ typedef typename Eigen::ei_traits<Derived>::Scalar Scalar; \
 typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
 typedef typename Base::PacketScalar PacketScalar; \
 typedef typename Eigen::ei_nested<Derived>::type Nested; \
-typedef typename Eigen::ei_eval<Derived>::type Eval; \
 enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
        ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \
        MaxRowsAtCompileTime = Eigen::ei_traits<Derived>::MaxRowsAtCompileTime, \
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index 67d1f8c1b..ae8703958 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -112,6 +112,10 @@ template<int _Rows, int _Cols> struct ei_size_at_compile_time
   enum { ret = (_Rows==Dynamic || _Cols==Dynamic) ? Dynamic : _Rows * _Cols };
 };
 
+/* ei_eval : the return type of eval(). For matrices, this is just a const reference
+ * in order to avoid a useless copy
+ */
+
 template<typename T, int Sparseness = ei_traits<T>::Flags&SparseBit> class ei_eval;
 
 template<typename T> struct ei_eval<T,IsDense>
@@ -125,8 +129,30 @@ template<typename T> struct ei_eval<T,IsDense>
           > type;
 };
 
+// for matrices, no need to evaluate, just use a const reference to avoid a useless copy
+template<typename _Scalar, int _Rows, int _Cols, int _StorageOrder, int _MaxRows, int _MaxCols>
+struct ei_eval<Matrix<_Scalar, _Rows, _Cols, _StorageOrder, _MaxRows, _MaxCols>, IsDense>
+{
+  typedef const Matrix<_Scalar, _Rows, _Cols, _StorageOrder, _MaxRows, _MaxCols>& type;
+};
 
-template<typename T> struct ei_eval_to_column_major
+/* ei_plain_matrix_type : the difference from ei_eval is that ei_plain_matrix_type is always a plain matrix type,
+ * whereas ei_eval is a const reference in the case of a matrix
+ */
+template<typename T> struct ei_plain_matrix_type
+{
+  typedef Matrix<typename ei_traits<T>::Scalar,
+                ei_traits<T>::RowsAtCompileTime,
+                ei_traits<T>::ColsAtCompileTime,
+                ei_traits<T>::Flags&RowMajorBit ? RowMajor : ColMajor,
+                ei_traits<T>::MaxRowsAtCompileTime,
+                ei_traits<T>::MaxColsAtCompileTime
+          > type;
+};
+
+/* ei_plain_matrix_type_column_major : same as ei_plain_matrix_type but guaranteed to be column-major
+ */
+template<typename T> struct ei_plain_matrix_type_column_major
 {
   typedef Matrix<typename ei_traits<T>::Scalar,
                 ei_traits<T>::RowsAtCompileTime,
@@ -158,7 +184,7 @@ template<typename T> struct ei_must_nest_by_value<NestByValue<T> > { enum { ret
   * const Matrix3d&, because the internal logic of ei_nested determined that since a was already a matrix, there was no point
   * in copying it into another matrix.
   */
-template<typename T, int n=1, typename EvalType = typename ei_eval<T>::type> struct ei_nested
+template<typename T, int n=1, typename PlainMatrixType = typename ei_eval<T>::type> struct ei_nested
 {
   enum {
     CostEval   = (n+1) * int(NumTraits<typename ei_traits<T>::Scalar>::ReadCost),
@@ -170,7 +196,7 @@ template<typename T, int n=1, typename EvalType = typename ei_eval<T>::type> str
     typename ei_meta_if<
       (int(ei_traits<T>::Flags) & EvalBeforeNestingBit)
       || ( int(CostEval) <= int(CostNoEval) ),
-      EvalType,
+      PlainMatrixType,
       const T&
     >::ret
   >::ret type;
diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h
index d60664cc0..047152d0b 100644
--- a/Eigen/src/Geometry/OrthoMethods.h
+++ b/Eigen/src/Geometry/OrthoMethods.h
@@ -34,7 +34,7 @@
   */
 template<typename Derived>
 template<typename OtherDerived>
-inline typename MatrixBase<Derived>::EvalType
+inline typename MatrixBase<Derived>::PlainMatrixType
 MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
 {
   EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3)
@@ -43,7 +43,7 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
   // optimize such a small temporary very well (even within a complex expression)
   const typename ei_nested<Derived,2>::type lhs(derived());
   const typename ei_nested<OtherDerived,2>::type rhs(other.derived());
-  return typename ei_eval<Derived>::type(
+  return typename ei_plain_matrix_type<Derived>::type(
     lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1),
     lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2),
     lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)
@@ -53,7 +53,7 @@ MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
 template<typename Derived, int Size = Derived::SizeAtCompileTime>
 struct ei_unitOrthogonal_selector
 {
-  typedef typename ei_eval<Derived>::type VectorType;
+  typedef typename ei_plain_matrix_type<Derived>::type VectorType;
   typedef typename ei_traits<Derived>::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   inline static VectorType run(const Derived& src)
@@ -96,7 +96,7 @@ struct ei_unitOrthogonal_selector
 template<typename Derived>
 struct ei_unitOrthogonal_selector<Derived,2>
 {
-  typedef typename ei_eval<Derived>::type VectorType;
+  typedef typename ei_plain_matrix_type<Derived>::type VectorType;
   inline static VectorType run(const Derived& src)
   { return VectorType(-ei_conj(src.y()), ei_conj(src.x())).normalized(); }
 };
@@ -109,7 +109,7 @@ struct ei_unitOrthogonal_selector<Derived,2>
   * \sa cross()
   */
 template<typename Derived>
-typename MatrixBase<Derived>::EvalType
+typename MatrixBase<Derived>::PlainMatrixType
 MatrixBase<Derived>::unitOrthogonal() const
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index d03d176f9..8be2df67a 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -548,7 +548,7 @@ inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const ScalingType
 {
   m_matrix.setZero();
   linear().diagonal() = s.coeffs();
-  m_matrix(Dim,Dim) = Scalar(1);
+  m_matrix.coeffRef(Dim,Dim) = Scalar(1);
   return *this;
 }
 
@@ -567,7 +567,7 @@ inline Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const RotationBas
   linear() = ei_toRotationMatrix<Scalar,Dim>(r);
   translation().setZero();
   m_matrix.template block<1,Dim>(Dim,0).setZero();
-  m_matrix(Dim,Dim) = Scalar(1);
+  m_matrix.coeffRef(Dim,Dim) = Scalar(1);
   return *this;
 }
 
@@ -610,7 +610,7 @@ Transform<Scalar,Dim>::extractRotation(TransformTraits traits) const
     LinearMatrixType matQ = qr.matrixQ();
     LinearMatrixType matR = qr.matrixR();
     for (int i=0 ; i<Dim; ++i)
-      if (matR(i,i)<0)
+      if (matR.coeff(i,i)<0)
         matQ.col(i) = -matQ.col(i);
     return matQ;
   }
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index cd59290ab..14d702cb1 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -92,7 +92,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp
     * R' = -S' * (R*P_inverse)
     */
   typedef Block<MatrixType,2,2> XprBlock22;
-  typedef typename XprBlock22::Eval Block22;
+  typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
   Block22 P_inverse;
   if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
   {
@@ -216,12 +216,11 @@ struct ei_compute_inverse<MatrixType, 4>
   * \sa inverse()
   */
 template<typename Derived>
-inline void MatrixBase<Derived>::computeInverse(EvalType *result) const
+inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
 {
-  typedef typename ei_eval<Derived>::type MatrixType;
   ei_assert(rows() == cols());
   EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,numeric_type_must_be_floating_point)
-  ei_compute_inverse<MatrixType>::run(eval(), result);
+  ei_compute_inverse<PlainMatrixType>::run(eval(), result);
 }
 
 /** \lu_module
@@ -239,9 +238,9 @@ inline void MatrixBase<Derived>::computeInverse(EvalType *result) const
   * \sa computeInverse()
   */
 template<typename Derived>
-inline const typename MatrixBase<Derived>::EvalType MatrixBase<Derived>::inverse() const
+inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const
 {
-  EvalType result(rows(), cols());
+  PlainMatrixType result(rows(), cols());
   computeInverse(&result);
   return result;
 }
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index bd2f30e84..41e33b48f 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -76,7 +76,7 @@ template<typename MatrixType> class LU
                   MatrixType::ColsAtCompileTime, // the number of rows in the "kernel matrix" is the number of cols of the original matrix
                                                  // so that the product "matrix * kernel = zero" makes sense
                   Dynamic,                       // we don't know at compile-time the dimension of the kernel
-                  MatrixType::Flags&RowMajorBit, // small optimization as we construct the kernel row by row
+                  MatrixType::Flags&RowMajorBit,
                   MatrixType::MaxColsAtCompileTime, // see explanation for 2nd template parameter
                   MatrixType::MaxColsAtCompileTime // the kernel is a subspace of the domain space, whose dimension is the number
                                                    // of columns of the original matrix
@@ -164,7 +164,8 @@ template<typename MatrixType> class LU
       *
       * \sa kernel(), computeImage(), image()
       */
-    void computeKernel(KernelResultType *result) const;
+    template<typename KernelMatrixType>
+    void computeKernel(KernelMatrixType *result) const;
 
     /** Computes a basis of the image of the matrix, also called the column-space or range of he matrix.
       *
@@ -179,7 +180,8 @@ template<typename MatrixType> class LU
       *
       * \sa image(), computeKernel(), kernel()
       */
-    void computeImage(ImageResultType *result) const;
+    template<typename ImageMatrixType>
+    void computeImage(ImageMatrixType *result) const;
 
     /** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
       * will form a basis of the kernel.
@@ -411,7 +413,8 @@ typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
 }
 
 template<typename MatrixType>
-void LU<MatrixType>::computeKernel(KernelResultType *result) const
+template<typename KernelMatrixType>
+void LU<MatrixType>::computeKernel(KernelMatrixType *result) const
 {
   ei_assert(!isInvertible());
   const int dimker = dimensionOfKernel(), cols = m_lu.cols();
@@ -434,7 +437,7 @@ void LU<MatrixType>::computeKernel(KernelResultType *result) const
     */
 
   Matrix<Scalar, Dynamic, Dynamic, MatrixType::Flags&RowMajorBit,
-         MatrixType::MaxColsAtCompileTime, MaxSmallDimAtCompileTime>
+         MatrixType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime>
     y(-m_lu.corner(TopRight, m_rank, dimker));
 
   m_lu.corner(TopLeft, m_rank, m_rank)
@@ -457,7 +460,8 @@ LU<MatrixType>::kernel() const
 }
 
 template<typename MatrixType>
-void LU<MatrixType>::computeImage(ImageResultType *result) const
+template<typename ImageMatrixType>
+void LU<MatrixType>::computeImage(ImageMatrixType *result) const
 {
   ei_assert(m_rank > 0);
   result->resize(m_originalMatrix.rows(), m_rank);
@@ -493,7 +497,7 @@ bool LU<MatrixType>::solve(
   ei_assert(b.rows() == rows);
   const int smalldim = std::min(rows, m_lu.cols());
 
-  typename OtherDerived::Eval c(b.rows(), b.cols());
+  typename OtherDerived::PlainMatrixType c(b.rows(), b.cols());
 
   // Step 1
   for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i);
@@ -540,10 +544,10 @@ bool LU<MatrixType>::solve(
   * \sa class LU
   */
 template<typename Derived>
-inline const LU<typename MatrixBase<Derived>::EvalType>
+inline const LU<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::lu() const
 {
-  return eval();
+  return LU<PlainMatrixType>(eval());
 }
 
 #endif // EIGEN_LU_H
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 94b817a02..c19205816 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -169,10 +169,10 @@ MatrixType QR<MatrixType>::matrixQ(void) const
   * \sa class QR
   */
 template<typename Derived>
-const QR<typename MatrixBase<Derived>::EvalType>
+const QR<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::qr() const
 {
-  return eval();
+  return QR<PlainMatrixType>(eval());
 }
 
 
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index 05060063c..c6bda1115 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -267,7 +267,7 @@ inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_
 MatrixBase<Derived>::eigenvalues() const
 {
   ei_assert(Flags&SelfAdjointBit);
-  return SelfAdjointEigenSolver<typename Derived::Eval>(eval(),false).eigenvalues();
+  return SelfAdjointEigenSolver<typename Derived::PlainMatrixType>(eval(),false).eigenvalues();
 }
 
 template<typename Derived, bool IsSelfAdjoint>
@@ -287,7 +287,7 @@ template<typename Derived> struct ei_operatorNorm_selector<Derived, false>
   static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
   operatorNorm(const MatrixBase<Derived>& m)
   {
-    typename Derived::Eval m_eval(m);
+    typename Derived::PlainMatrixType m_eval(m);
     // FIXME if it is really guaranteed that the eigenvalues are already sorted,
     // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
     return ei_sqrt(
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 988316649..b8432c943 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -538,10 +538,10 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul
   * \returns the SVD decomposition of \c *this
   */
 template<typename Derived>
-inline SVD<typename MatrixBase<Derived>::EvalType>
+inline SVD<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::svd() const
 {
-  return SVD<typename ei_eval<Derived>::type>(derived());
+  return SVD<PlainMatrixType>(derived());
 }
 
 #endif // EIGEN_SVD_H
diff --git a/doc/snippets/LU_computeImage.cpp b/doc/snippets/LU_computeImage.cpp
new file mode 100644
index 000000000..5c812cc4c
--- /dev/null
+++ b/doc/snippets/LU_computeImage.cpp
@@ -0,0 +1,13 @@
+MatrixXd m(3,3);
+m << 1,1,0,
+     1,3,2,
+     0,1,1;
+cout << "Here is the matrix m:" << endl << m << endl;
+LU<Matrix3d> lu(m);
+// allocate the matrix img with the correct size to avoid reallocation
+MatrixXd img(m.rows(), lu.rank());
+lu.computeImage(&img);
+cout << "Notice that the middle column is the sum of the two others, so the "
+     << "columns are linearly dependent." << endl;
+cout << "Here is a matrix whose columns have the same span but are linearly independent:"
+     << endl << img << endl;
diff --git a/doc/snippets/LU_image.cpp b/doc/snippets/LU_image.cpp
new file mode 100644
index 000000000..0d1088a2f
--- /dev/null
+++ b/doc/snippets/LU_image.cpp
@@ -0,0 +1,9 @@
+Matrix3d m;
+m << 1,1,0,
+     1,3,2,
+     0,1,1;
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Notice that the middle column is the sum of the two others, so the "
+     << "columns are linearly dependent." << endl;
+cout << "Here is a matrix whose columns have the same span but are linearly independent:"
+     << endl << m.lu().image() << endl;
diff --git a/test/main.h b/test/main.h
index 433c8688b..ec6724ffd 100644
--- a/test/main.h
+++ b/test/main.h
@@ -162,7 +162,7 @@ namespace Eigen {
 
 template<typename T> inline typename NumTraits<T>::Real test_precision();
 template<> inline int test_precision<int>() { return 0; }
-template<> inline float test_precision<float>() { return 1e-4f; }
+template<> inline float test_precision<float>() { return 1e-3f; }
 template<> inline double test_precision<double>() { return 1e-6; }
 template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); }
 template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); }