sync with default branch

31 files changed, 337 insertions, 68 deletions

diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h
index e846b38e5..2f94d12dc 100644
--- a/Eigen/src/Core/BandMatrix.h
+++ b/Eigen/src/Core/BandMatrix.h
@@ -55,7 +55,7 @@ struct ei_traits<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
     ColsAtCompileTime = Cols,
     MaxRowsAtCompileTime = Rows,
     MaxColsAtCompileTime = Cols,
-    Flags = 0
+    Flags = LvalueBit
   };
 };
 
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index f888d5ea0..498757ad9 100644
--- a/Eigen/src/Core/Block.h
+++ b/Eigen/src/Core/Block.h
@@ -93,7 +93,7 @@ struct ei_traits<Block<XprType, BlockRows, BlockCols, HasDirectAccess> > : ei_tr
                        && (InnerStrideAtCompileTime == 1)
                         ? PacketAccessBit : 0,
     FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0,
-    Flags0 = ei_traits<XprType>::Flags & (HereditaryBits | MaskPacketAccessBit | DirectAccessBit),
+    Flags0 = ei_traits<XprType>::Flags & (HereditaryBits | MaskPacketAccessBit | LvalueBit | DirectAccessBit),
     Flags1 = Flags0 | FlagsLinearAccessBit,
     Flags = (Flags1 & ~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0)
   };
@@ -679,9 +679,62 @@ DenseBase<Derived>::bottomRows() const
 
 
 
+/** \returns a block consisting of a range of rows of *this.
+  *
+  * \param startRow the index of the first row in the block
+  * \param numRows the number of rows in the block
+  *
+  * Example: \include MatrixBase_middleRows_int.cpp
+  * Output: \verbinclude MatrixBase_middleRows_int.out
+  *
+  * \sa class Block, block(Index,Index,Index,Index)
+  */
+template<typename Derived>
+inline typename DenseBase<Derived>::RowsBlockXpr DenseBase<Derived>
+  ::middleRows(Index startRow, Index numRows)
+{
+  return RowsBlockXpr(derived(), startRow, 0, numRows, cols());
+}
 
+/** This is the const version of middleRows(Index,Index).*/
+template<typename Derived>
+inline const typename DenseBase<Derived>::RowsBlockXpr
+DenseBase<Derived>::middleRows(Index startRow, Index numRows) const
+{
+  return RowsBlockXpr(derived(), startRow, 0, numRows, cols());
+}
 
-/** \returns a block consisting of the top columns of *this.
+/** \returns a block consisting of a range of rows of *this.
+  *
+  * \param N the number of rows in the block
+  * \param startRow the index of the first row in the block
+  *
+  * Example: \include MatrixBase_template_int_middleRows.cpp
+  * Output: \verbinclude MatrixBase_template_int_middleRows.out
+  *
+  * \sa class Block, block(Index,Index,Index,Index)
+  */
+template<typename Derived>
+template<int N>
+inline typename DenseBase<Derived>::template NRowsBlockXpr<N>::Type
+DenseBase<Derived>::middleRows(Index startRow)
+{
+  return typename DenseBase<Derived>::template NRowsBlockXpr<N>::Type(derived(), startRow, 0, N, cols());
+}
+
+/** This is the const version of middleRows<int>().*/
+template<typename Derived>
+template<int N>
+inline const typename DenseBase<Derived>::template NRowsBlockXpr<N>::Type
+DenseBase<Derived>::middleRows(Index startRow) const
+{
+  return typename DenseBase<Derived>::template NRowsBlockXpr<N>::Type(derived(), startRow, 0, N, cols());
+}
+
+
+
+
+/** \returns a block consisting of the left columns of *this.
   *
   * \param n the number of columns in the block
   *
@@ -788,6 +841,61 @@ DenseBase<Derived>::rightCols() const
 
 
 
+/** \returns a block consisting of a range of columns of *this.
+  *
+  * \param startCol the index of the first column in the block
+  * \param numCols the number of columns in the block
+  *
+  * Example: \include MatrixBase_middleCols_int.cpp
+  * Output: \verbinclude MatrixBase_middleCols_int.out
+  *
+  * \sa class Block, block(Index,Index,Index,Index)
+  */
+template<typename Derived>
+inline typename DenseBase<Derived>::ColsBlockXpr DenseBase<Derived>
+  ::middleCols(Index startCol, Index numCols)
+{
+  return ColsBlockXpr(derived(), 0, startCol, rows(), numCols);
+}
+
+/** This is the const version of middleCols(Index,Index).*/
+template<typename Derived>
+inline const typename DenseBase<Derived>::ColsBlockXpr
+DenseBase<Derived>::middleCols(Index startCol, Index numCols) const
+{
+  return ColsBlockXpr(derived(), 0, startCol, rows(), numCols);
+}
+
+/** \returns a block consisting of a range of columns of *this.
+  *
+  * \param N the number of columns in the block
+  * \param startCol the index of the first column in the block
+  *
+  * Example: \include MatrixBase_template_int_middleCols.cpp
+  * Output: \verbinclude MatrixBase_template_int_middleCols.out
+  *
+  * \sa class Block, block(Index,Index,Index,Index)
+  */
+template<typename Derived>
+template<int N>
+inline typename DenseBase<Derived>::template NColsBlockXpr<N>::Type
+DenseBase<Derived>::middleCols(Index startCol)
+{
+  return typename NColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), N);
+}
+
+/** This is the const version of middleCols<int>().*/
+template<typename Derived>
+template<int N>
+inline const typename DenseBase<Derived>::template NColsBlockXpr<N>::Type
+DenseBase<Derived>::middleCols(Index startCol) const
+{
+  return typename NColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), N);
+}
+
+
+
+
 
 /** \returns a fixed-size expression of a block in *this.
   *
diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h
index 25041d1b7..97331fa26 100644
--- a/Eigen/src/Core/CwiseNullaryOp.h
+++ b/Eigen/src/Core/CwiseNullaryOp.h
@@ -240,17 +240,30 @@ DenseBase<Derived>::Constant(const Scalar& value)
   * Example: \include DenseBase_LinSpaced_seq.cpp
   * Output: \verbinclude DenseBase_LinSpaced_seq.out
   *
-  * \sa setLinSpaced(const Scalar&,const Scalar&,Index), LinSpaced(Scalar,Scalar,Index), CwiseNullaryOp
+  * \sa setLinSpaced(Index,const Scalar&,const Scalar&), LinSpaced(Index,Scalar,Scalar), CwiseNullaryOp
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::SequentialLinSpacedReturnType
-DenseBase<Derived>::LinSpaced(Sequential_t, const Scalar& low, const Scalar& high, Index size)
+DenseBase<Derived>::LinSpaced(Sequential_t, Index size, const Scalar& low, const Scalar& high)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
   return DenseBase<Derived>::NullaryExpr(size, ei_linspaced_op<Scalar,false>(low,high,size));
 }
 
 /**
+  * \copydoc DenseBase<Derived>::LinSpaced(Sequential_t, Index, const Scalar&, const Scalar&)
+  * Special version for fixed size types which does not require the size parameter.
+  */
+template<typename Derived>
+EIGEN_STRONG_INLINE const typename DenseBase<Derived>::SequentialLinSpacedReturnType
+DenseBase<Derived>::LinSpaced(Sequential_t, const Scalar& low, const Scalar& high)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived)
+  return DenseBase<Derived>::NullaryExpr(Derived::SizeAtCompileTime, ei_linspaced_op<Scalar,false>(low,high,Derived::SizeAtCompileTime));
+}
+
+/**
   * \brief Sets a linearly space vector.
   *
   * The function generates 'size' equally spaced values in the closed interval [low,high].
@@ -260,16 +273,29 @@ DenseBase<Derived>::LinSpaced(Sequential_t, const Scalar& low, const Scalar& hig
   * Example: \include DenseBase_LinSpaced.cpp
   * Output: \verbinclude DenseBase_LinSpaced.out
   *
-  * \sa setLinSpaced(const Scalar&,const Scalar&,Index), LinSpaced(Sequential_t,const Scalar&,const Scalar&,Index), CwiseNullaryOp
+  * \sa setLinSpaced(Index,const Scalar&,const Scalar&), LinSpaced(Sequential_t,Index,const Scalar&,const Scalar&,Index), CwiseNullaryOp
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::RandomAccessLinSpacedReturnType
-DenseBase<Derived>::LinSpaced(const Scalar& low, const Scalar& high, Index size)
+DenseBase<Derived>::LinSpaced(Index size, const Scalar& low, const Scalar& high)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
   return DenseBase<Derived>::NullaryExpr(size, ei_linspaced_op<Scalar,true>(low,high,size));
 }
 
+/**
+  * \copydoc DenseBase<Derived>::LinSpaced(Index, const Scalar&, const Scalar&)
+  * Special version for fixed size types which does not require the size parameter.
+  */
+template<typename Derived>
+EIGEN_STRONG_INLINE const typename DenseBase<Derived>::RandomAccessLinSpacedReturnType
+DenseBase<Derived>::LinSpaced(const Scalar& low, const Scalar& high)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived)
+  return DenseBase<Derived>::NullaryExpr(Derived::SizeAtCompileTime, ei_linspaced_op<Scalar,true>(low,high,Derived::SizeAtCompileTime));
+}
+
 /** \returns true if all coefficients in this matrix are approximately equal to \a value, to within precision \a prec */
 template<typename Derived>
 bool DenseBase<Derived>::isApproxToConstant
@@ -360,7 +386,7 @@ DenseStorageBase<Derived>::setConstant(Index rows, Index cols, const Scalar& val
   * \sa CwiseNullaryOp
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(const Scalar& low, const Scalar& high, Index size)
+EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(Index size, const Scalar& low, const Scalar& high)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
   return derived() = Derived::NullaryExpr(size, ei_linspaced_op<Scalar,false>(low,high,size));
diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h
index 612b28f1f..37c58223e 100644
--- a/Eigen/src/Core/CwiseUnaryView.h
+++ b/Eigen/src/Core/CwiseUnaryView.h
@@ -48,7 +48,7 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> >
   typedef typename MatrixType::Nested MatrixTypeNested;
   typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested;
   enum {
-    Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | DirectAccessBit)),
+    Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LvalueBit | LinearAccessBit | DirectAccessBit)),
     CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits<ViewOp>::Cost,
     MatrixTypeInnerStride =  ei_inner_stride_at_compile_time<MatrixType>::ret,
     // need to cast the sizeof's from size_t to int explicitly, otherwise:
diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h
index 8c1e24fec..c84544ec3 100644
--- a/Eigen/src/Core/DenseBase.h
+++ b/Eigen/src/Core/DenseBase.h
@@ -327,10 +327,14 @@ template<typename Derived> class DenseBase
     const RowsBlockXpr topRows(Index n) const;
     RowsBlockXpr       bottomRows(Index n);
     const RowsBlockXpr bottomRows(Index n) const;
+	RowsBlockXpr       middleRows(Index startRow, Index numRows);
+	const RowsBlockXpr middleRows(Index startRow, Index numRows) const;
     ColsBlockXpr       leftCols(Index n);
     const ColsBlockXpr leftCols(Index n) const;
     ColsBlockXpr       rightCols(Index n);
     const ColsBlockXpr rightCols(Index n) const;
+	ColsBlockXpr       middleCols(Index startCol, Index numCols);
+	const ColsBlockXpr middleCols(Index startCol, Index numCols) const;
 
     template<int CRows, int CCols> Block<Derived, CRows, CCols>       topLeftCorner();
     template<int CRows, int CCols> const Block<Derived, CRows, CCols> topLeftCorner() const;
@@ -345,10 +349,14 @@ template<typename Derived> class DenseBase
     template<int NRows> const typename NRowsBlockXpr<NRows>::Type topRows() const;
     template<int NRows> typename NRowsBlockXpr<NRows>::Type       bottomRows();
     template<int NRows> const typename NRowsBlockXpr<NRows>::Type bottomRows() const;
+	template<int NRows> typename NRowsBlockXpr<NRows>::Type       middleRows(Index startRow);
+	template<int NRows> const typename NRowsBlockXpr<NRows>::Type middleRows(Index startRow) const;
     template<int NCols> typename NColsBlockXpr<NCols>::Type       leftCols();
     template<int NCols> const typename NColsBlockXpr<NCols>::Type leftCols() const;
     template<int NCols> typename NColsBlockXpr<NCols>::Type       rightCols();
     template<int NCols> const typename NColsBlockXpr<NCols>::Type rightCols() const;
+	template<int NCols> typename NColsBlockXpr<NCols>::Type       middleCols(Index startCol);
+	template<int NCols> const typename NColsBlockXpr<NCols>::Type middleCols(Index startCol) const;
 
     template<int BlockRows, int BlockCols>
     Block<Derived, BlockRows, BlockCols> block(Index startRow, Index startCol);
@@ -390,9 +398,13 @@ template<typename Derived> class DenseBase
     Constant(const Scalar& value);
 
     static const SequentialLinSpacedReturnType
-    LinSpaced(Sequential_t, const Scalar& low, const Scalar& high, Index size);
+    LinSpaced(Sequential_t, Index size, const Scalar& low, const Scalar& high);
     static const RandomAccessLinSpacedReturnType
-    LinSpaced(const Scalar& low, const Scalar& high, Index size);
+    LinSpaced(Index size, const Scalar& low, const Scalar& high);
+    static const SequentialLinSpacedReturnType
+    LinSpaced(Sequential_t, const Scalar& low, const Scalar& high);
+    static const RandomAccessLinSpacedReturnType
+    LinSpaced(const Scalar& low, const Scalar& high);
 
     template<typename CustomNullaryOp>
     static const CwiseNullaryOp<CustomNullaryOp, Derived>
@@ -413,7 +425,8 @@ template<typename Derived> class DenseBase
 
     void fill(const Scalar& value);
     Derived& setConstant(const Scalar& value);
-    Derived& setLinSpaced(const Scalar& low, const Scalar& high, Index size);
+    Derived& setLinSpaced(Index size, const Scalar& low, const Scalar& high);
+    Derived& setLinSpaced(const Scalar& low, const Scalar& high);
     Derived& setZero();
     Derived& setOnes();
     Derived& setRandom();
diff --git a/Eigen/src/Core/DenseCoeffsBase.h b/Eigen/src/Core/DenseCoeffsBase.h
index 6802cea24..0b348c81d 100644
--- a/Eigen/src/Core/DenseCoeffsBase.h
+++ b/Eigen/src/Core/DenseCoeffsBase.h
@@ -25,15 +25,15 @@
 #ifndef EIGEN_DENSECOEFFSBASE_H
 #define EIGEN_DENSECOEFFSBASE_H
 
-template<typename Derived, bool EnableDirectAccessAPI>
-class DenseCoeffsBase : public EigenBase<Derived>
+template<typename Derived>
+class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
 {
   public:
     typedef typename ei_traits<Derived>::StorageKind StorageKind;
     typedef typename ei_traits<Derived>::Index Index;
     typedef typename ei_traits<Derived>::Scalar Scalar;
     typedef typename ei_packet_traits<Scalar>::type PacketScalar;
-    typedef typename ei_meta_if<ei_has_direct_access<Derived>::ret,
+    typedef typename ei_meta_if<bool(ei_traits<Derived>::Flags&LvalueBit),
                                 const Scalar&,
                                 typename ei_meta_if<ei_is_arithmetic<Scalar>::ret, Scalar, const Scalar>::ret
                                >::ret CoeffReturnType;
@@ -239,11 +239,11 @@ class DenseCoeffsBase : public EigenBase<Derived>
 };
 
 template<typename Derived>
-class DenseCoeffsBase<Derived, true> : public DenseCoeffsBase<Derived, false>
+class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived, ReadOnlyAccessors>
 {
   public:
 
-    typedef DenseCoeffsBase<Derived, false> Base;
+    typedef DenseCoeffsBase<Derived, ReadOnlyAccessors> Base;
 
     typedef typename ei_traits<Derived>::StorageKind StorageKind;
     typedef typename ei_traits<Derived>::Index Index;
@@ -512,6 +512,23 @@ class DenseCoeffsBase<Derived, true> : public DenseCoeffsBase<Derived, false>
     }
 #endif
 
+};
+
+template<typename Derived>
+class DenseCoeffsBase<Derived, DirectAccessors> : public DenseCoeffsBase<Derived, WriteAccessors>
+{
+  public:
+
+    typedef DenseCoeffsBase<Derived, WriteAccessors> Base;
+    typedef typename ei_traits<Derived>::Index Index;
+    typedef typename ei_traits<Derived>::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+
+    using Base::rows;
+    using Base::cols;
+    using Base::size;
+    using Base::derived;
+
     /** \returns the pointer increment between two consecutive elements within a slice in the inner direction.
       *
       * \sa outerStride(), rowStride(), colStride()
@@ -531,6 +548,7 @@ class DenseCoeffsBase<Derived, true> : public DenseCoeffsBase<Derived, false>
       return derived().outerStride();
     }
 
+    // FIXME shall we remove it ?
     inline Index stride() const
     {
       return Derived::IsVectorAtCompileTime ? innerStride() : outerStride();
diff --git a/Eigen/src/Core/DenseStorageBase.h b/Eigen/src/Core/DenseStorageBase.h
index 16e0a86f1..e818d40a1 100644
--- a/Eigen/src/Core/DenseStorageBase.h
+++ b/Eigen/src/Core/DenseStorageBase.h
@@ -483,8 +483,8 @@ class DenseStorageBase : public ei_dense_xpr_base<Derived>::type
     template<typename T0, typename T1>
     EIGEN_STRONG_INLINE void _init2(Index rows, Index cols, typename ei_enable_if<Base::SizeAtCompileTime!=2,T0>::type* = 0)
     {
-      ei_assert(rows > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
-             && cols > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
+      ei_assert(rows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
+             && cols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
       m_storage.resize(rows*cols,rows,cols);
       EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
     }
diff --git a/Eigen/src/Core/Diagonal.h b/Eigen/src/Core/Diagonal.h
index 8c2eacd96..0b7d14179 100644
--- a/Eigen/src/Core/Diagonal.h
+++ b/Eigen/src/Core/Diagonal.h
@@ -62,7 +62,7 @@ struct ei_traits<Diagonal<MatrixType,DiagIndex> >
                                                                     MatrixType::MaxColsAtCompileTime)
                          : (EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime) - AbsDiagIndex),
     MaxColsAtCompileTime = 1,
-    Flags = (unsigned int)_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | DirectAccessBit) & ~RowMajorBit,
+    Flags = (unsigned int)_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | LvalueBit | DirectAccessBit) & ~RowMajorBit,
     CoeffReadCost = _MatrixTypeNested::CoeffReadCost,
     MatrixTypeOuterStride = ei_outer_stride_at_compile_time<MatrixType>::ret,
     InnerStrideAtCompileTime = MatrixTypeOuterStride == Dynamic ? Dynamic : MatrixTypeOuterStride+1,
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 8c4200a6f..0c0525028 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -105,6 +105,9 @@ struct ei_traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime>
   typedef Matrix<_Scalar,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1> DiagonalVectorType;
   typedef Dense StorageKind;
   typedef DenseIndex Index;
+  enum {
+    Flags = LvalueBit
+  };
 };
 
 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
@@ -213,7 +216,7 @@ struct ei_traits<DiagonalWrapper<_DiagonalVectorType> >
     ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
     MaxRowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
     MaxColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
-    Flags = 0
+    Flags =  ei_traits<DiagonalVectorType>::Flags & LvalueBit
   };
 };
 
diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h
index cde91ff97..41ae4af42 100644
--- a/Eigen/src/Core/Functors.h
+++ b/Eigen/src/Core/Functors.h
@@ -577,7 +577,7 @@ template <typename Scalar, bool RandomAccess> struct ei_linspaced_op
   template<typename Index>
   EIGEN_STRONG_INLINE const Packet packetOp(Index i, Index = 0) const { return impl.packetOp(i); }
   // This proxy object handles the actual required temporaries, the different
-  // implementations (random vs. sequential access) as well as the piping
+  // implementations (random vs. sequential access) as well as the
   // correct piping to size 2/4 packet operations.
   const ei_linspaced_op_impl<Scalar,RandomAccess> impl;
 };
diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h
index b1dd7d727..504a51229 100644
--- a/Eigen/src/Core/Redux.h
+++ b/Eigen/src/Core/Redux.h
@@ -183,7 +183,7 @@ struct ei_redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>
   typedef typename Derived::Index Index;
   static EIGEN_STRONG_INLINE Scalar run(const Derived& mat, const Func& func)
   {
-    ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix");
+    ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
     Scalar res;
     res = mat.coeffByOuterInner(0, 0);
     for(Index i = 1; i < mat.innerSize(); ++i)
@@ -210,6 +210,7 @@ struct ei_redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling>
   static Scalar run(const Derived& mat, const Func& func)
   {
     const Index size = mat.size();
+    ei_assert(size && "you are using an empty matrix");
     const Index packetSize = ei_packet_traits<Scalar>::size;
     const Index alignedStart = ei_first_aligned(mat);
     enum {
@@ -253,6 +254,7 @@ struct ei_redux_impl<Func, Derived, SliceVectorizedTraversal, NoUnrolling>
 
   static Scalar run(const Derived& mat, const Func& func)
   {
+    ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
     const Index innerSize = mat.innerSize();
     const Index outerSize = mat.outerSize();
     enum {
@@ -294,6 +296,7 @@ struct ei_redux_impl<Func, Derived, LinearVectorizedTraversal, CompleteUnrolling
   };
   EIGEN_STRONG_INLINE static Scalar run(const Derived& mat, const Func& func)
   {
+    ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
     Scalar res = func.predux(ei_redux_vec_unroller<Func, Derived, 0, Size / PacketSize>::run(mat,func));
     if (VectorizedSize != Size)
       res = func(res,ei_redux_novec_unroller<Func, Derived, VectorizedSize, Size-VectorizedSize>::run(mat,func));
@@ -345,6 +348,8 @@ template<typename Derived>
 EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar
 DenseBase<Derived>::sum() const
 {
+  if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
+    return Scalar(0);
   return this->redux(Eigen::ei_scalar_sum_op<Scalar>());
 }
 
@@ -370,6 +375,8 @@ template<typename Derived>
 EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar
 DenseBase<Derived>::prod() const
 {
+  if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
+    return Scalar(1);
   return this->redux(Eigen::ei_scalar_product_op<Scalar>());
 }
 
diff --git a/Eigen/src/Core/Reverse.h b/Eigen/src/Core/Reverse.h
index 5a96aeeb3..6051f7366 100644
--- a/Eigen/src/Core/Reverse.h
+++ b/Eigen/src/Core/Reverse.h
@@ -59,7 +59,7 @@ struct ei_traits<Reverse<MatrixType, Direction> >
     LinearAccess = ( (Direction==BothDirections) && (int(_MatrixTypeNested::Flags)&PacketAccessBit) )
                  ? LinearAccessBit : 0,
 
-    Flags = int(_MatrixTypeNested::Flags) & (HereditaryBits | PacketAccessBit | LinearAccess),
+    Flags = int(_MatrixTypeNested::Flags) & (HereditaryBits | LvalueBit | PacketAccessBit | LinearAccess),
 
     CoeffReadCost = _MatrixTypeNested::CoeffReadCost
   };
@@ -109,6 +109,11 @@ template<typename MatrixType, int Direction> class Reverse
     inline Index rows() const { return m_matrix.rows(); }
     inline Index cols() const { return m_matrix.cols(); }
 
+    inline Index innerStride() const
+    {
+      return -m_matrix.innerStride();
+    }
+
     inline Scalar& operator()(Index row, Index col)
     {
       ei_assert(row >= 0 && row < rows() && col >= 0 && col < cols());
diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
index e3b9d25ef..7e2d496fe 100644
--- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h
+++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h
@@ -126,7 +126,7 @@ void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrd
 
   ei_manage_caching_sizes(GetAction, &l1, &l2);
   k = std::min<std::ptrdiff_t>(k, l1/kdiv);
-  std::ptrdiff_t _m = l2/(4 * sizeof(LhsScalar) * k);
+  std::ptrdiff_t _m = k>0 ? l2/(4 * sizeof(LhsScalar) * k) : 0;
   if(_m<m) m = _m & mr_mask;
   n = n;
 }
diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h
index 8222123d9..44986dc16 100644
--- a/Eigen/src/Core/products/GeneralMatrixVector.h
+++ b/Eigen/src/Core/products/GeneralMatrixVector.h
@@ -319,6 +319,7 @@ EIGEN_DONT_INLINE static void run(
   ResScalar* res, Index resIncr,
   ResScalar alpha)
 {
+  EIGEN_UNUSED_VARIABLE(rhsIncr);
   ei_internal_assert(rhsIncr==1);
   #ifdef _EIGEN_ACCUMULATE_PACKETS
   #error _EIGEN_ACCUMULATE_PACKETS has already been defined
diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
index 0fbdd024c..ede8b77bf 100644
--- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
+++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h
@@ -135,7 +135,7 @@ struct ei_symm_pack_rhs
         for (Index w=0 ; w<h; ++w)
           blockB[count+w] = rhs(k,j2+w);
 
-        blockB[count+h] = rhs(k,k);
+        blockB[count+h] = ei_real(rhs(k,k));
 
         // transpose
         for (Index w=h+1 ; w<nr; ++w)
diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index 0f90b5683..cef5eeba1 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -75,7 +75,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,LhsIsTriangular,
     Scalar alpha)
   {
     ei_product_triangular_matrix_matrix<Scalar, Index,
-      (Mode&UnitDiag) | ((Mode&Upper) ? Lower : Upper),
+      (Mode&(UnitDiag|ZeroDiag)) | ((Mode&Upper) ? Lower : Upper),
       (!LhsIsTriangular),
       RhsStorageOrder==RowMajor ? ColMajor : RowMajor,
       ConjugateRhs,
@@ -108,7 +108,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
     typedef ei_gebp_traits<Scalar,Scalar> Traits;
     enum {
       SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower
+      IsLower = (Mode&Lower) == Lower,
+      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
     };
 
     Index kc = depth; // cache block size along the K direction
@@ -124,7 +125,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
     triangularBuffer.setZero();
-    triangularBuffer.diagonal().setOnes();
+    if((Mode&ZeroDiag)==ZeroDiag)
+      triangularBuffer.diagonal().setZero();
+    else
+      triangularBuffer.diagonal().setOnes();
 
     ei_gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
     ei_gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
@@ -166,7 +170,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,true,
           // To this end we do an extra triangular copy to a small temporary buffer
           for (Index k=0;k<actualPanelWidth;++k)
           {
-            if (!(Mode&UnitDiag))
+            if (SetDiag)
               triangularBuffer.coeffRef(k,k) = lhs(startBlock+k,startBlock+k);
             for (Index i=IsLower ? k+1 : 0; IsLower ? i<actualPanelWidth : i<k; ++i)
               triangularBuffer.coeffRef(i,k) = lhs(startBlock+i,startBlock+k);
@@ -231,7 +235,8 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
     typedef ei_gebp_traits<Scalar,Scalar> Traits;
     enum {
       SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower
+      IsLower = (Mode&Lower) == Lower,
+      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
     };
 
     Index kc = depth; // cache block size along the K direction
@@ -247,7 +252,10 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
     triangularBuffer.setZero();
-    triangularBuffer.diagonal().setOnes();
+    if((Mode&ZeroDiag)==ZeroDiag)
+      triangularBuffer.diagonal().setZero();
+    else
+      triangularBuffer.diagonal().setOnes();
 
     ei_gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
     ei_gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
@@ -295,7 +303,7 @@ struct ei_product_triangular_matrix_matrix<Scalar,Index,Mode,false,
           // append the triangular part via a temporary buffer
           for (Index j=0;j<actualPanelWidth;++j)
           {
-            if (!(Mode&UnitDiag))
+            if (SetDiag)
               triangularBuffer.coeffRef(j,j) = rhs(actual_j2+j,actual_j2+j);
             for (Index k=IsLower ? j+1 : 0; IsLower ? k<actualPanelWidth : k<j; ++k)
               triangularBuffer.coeffRef(k,j) = rhs(actual_j2+k,actual_j2+j);
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index eb70981e4..60fdcf2e2 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -139,6 +139,14 @@ const unsigned int DirectAccessBit = 0x20;
   * means the first coefficient packet is guaranteed to be aligned */
 const unsigned int AlignedBit = 0x40;
 
+/** \ingroup flags
+  *
+  * Means the expression is writable. Note that DirectAccessBit implies LvalueBit.
+  * Internaly, it is mainly used to enable the writable coeff accessors, and makes
+  * the read-only coeff accessors to return by const reference.
+  */
+const unsigned int LvalueBit = 0x80;
+
 const unsigned int NestByRefBit = 0x100;
 
 // list of flags that are inherited by default
@@ -236,6 +244,10 @@ enum {
   IsSparse
 };
 
+enum AccessorLevels {
+  ReadOnlyAccessors, WriteAccessors, DirectAccessors
+};
+
 enum DecompositionOptions {
   Pivoting            = 0x01, // LDLT,
   NoPivoting          = 0x02, // LDLT,
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 2c0e9c3e2..7f626c62b 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -36,7 +36,10 @@ template<typename Derived> struct ei_has_direct_access
 
 template<typename Derived> struct EigenBase;
 template<typename Derived> class DenseBase;
-template<typename Derived, bool EnableDirectAccessAPI = ei_has_direct_access<Derived>::ret>
+template<typename Derived,
+         AccessorLevels Level = (ei_traits<Derived>::Flags & DirectAccessBit) ? DirectAccessors
+                              : (ei_traits<Derived>::Flags & LvalueBit) ? WriteAccessors
+                              : ReadOnlyAccessors>
 class DenseCoeffsBase;
 
 template<typename _Scalar, int _Rows, int _Cols,
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index 1b89525a8..75c3e5ec4 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -315,7 +315,8 @@ template<typename T> inline T* ei_construct_elements_of_array(T *ptr, size_t siz
 template<typename T> inline void ei_destruct_elements_of_array(T *ptr, size_t size)
 {
   // always destruct an array starting from the end.
-  while(size) ptr[--size].~T();
+  if(ptr)
+    while(size) ptr[--size].~T();
 }
 
 /*****************************************************************************
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h
index 9d3169967..bf686975b 100644
--- a/Eigen/src/Core/util/StaticAssert.h
+++ b/Eigen/src/Core/util/StaticAssert.h
@@ -60,6 +60,7 @@
         YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES,
         THIS_METHOD_IS_ONLY_FOR_VECTORS_OF_A_SPECIFIC_SIZE,
         THIS_METHOD_IS_ONLY_FOR_MATRICES_OF_A_SPECIFIC_SIZE,
+        THIS_METHOD_IS_ONLY_FOR_OBJECTS_OF_A_SPECIFIC_SIZE,
         YOU_MADE_A_PROGRAMMING_MISTAKE,
         EIGEN_INTERNAL_COMPILATION_ERROR_OR_YOU_MADE_A_PROGRAMMING_MISTAKE,
         YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR,
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index c93398092..3c89a7f3f 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -155,7 +155,7 @@ class ei_compute_matrix_flags
     };
 
   public:
-    enum { ret = LinearAccessBit | DirectAccessBit | NestByRefBit | packet_access_bit | row_major_bit | aligned_bit };
+    enum { ret = LinearAccessBit | LvalueBit | DirectAccessBit | NestByRefBit | packet_access_bit | row_major_bit | aligned_bit };
 };
 
 template<int _Rows, int _Cols> struct ei_size_at_compile_time
@@ -355,7 +355,7 @@ template<typename T, int n=1, typename PlainObject = typename ei_eval<T>::type>
 
 template<unsigned int Flags> struct ei_are_flags_consistent
 {
-  enum { ret = true };
+  enum { ret = EIGEN_IMPLIES(bool(Flags&DirectAccessBit), bool(Flags&LvalueBit)) };
 };
 
 template<typename Derived, typename XprKind = typename ei_traits<Derived>::XprKind>
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index 08276c5c0..7bf1d140e 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -291,8 +291,8 @@ void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm
       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
       if(z==ComplexScalar(0))
       {
-	// If the i-th and k-th eigenvalue are equal, then z equals 0.
-	// Use a small value instead, to prevent division by zero.
+        // If the i-th and k-th eigenvalue are equal, then z equals 0.
+        // Use a small value instead, to prevent division by zero.
         ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
       }
       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index 0f2b78e27..79554187a 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -130,7 +130,7 @@ template<typename _MatrixType> class HessenbergDecomposition
     {
       if(matrix.rows()<2)
       {
-	m_isInitialized = true;
+        m_isInitialized = true;
         return;
       }
       m_hCoeffs.resize(matrix.rows()-1,1);
@@ -160,7 +160,7 @@ template<typename _MatrixType> class HessenbergDecomposition
       m_matrix = matrix;
       if(matrix.rows()<2)
       {
-	m_isInitialized = true;
+        m_isInitialized = true;
         return *this;
       }
       m_hCoeffs.resize(matrix.rows()-1,1);
@@ -360,7 +360,7 @@ template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
       result = m_hess.packedMatrix();
       Index n = result.rows();
       if (n>2)
-	result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
+        result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
     }
 
     Index rows() const { return m_hess.packedMatrix().rows(); }
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index 611b89730..4211981af 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -384,7 +384,9 @@ void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
 }
 
 // forward declaration, implementation at the end of this file
-template<typename MatrixType, int Size=MatrixType::ColsAtCompileTime>
+template<typename MatrixType,
+         int Size=MatrixType::ColsAtCompileTime,
+         bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
 struct ei_tridiagonalization_inplace_selector;
 
 /** \brief Performs a full tridiagonalization in place
@@ -439,7 +441,7 @@ void ei_tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiago
 /** \internal
   * General full tridiagonalization
   */
-template<typename MatrixType, int Size>
+template<typename MatrixType, int Size, bool IsComplex>
 struct ei_tridiagonalization_inplace_selector
 {
   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
@@ -458,11 +460,11 @@ struct ei_tridiagonalization_inplace_selector
 };
 
 /** \internal
-  * Specialization for 3x3 matrices.
+  * Specialization for 3x3 real matrices.
   * Especially useful for plane fitting.
   */
 template<typename MatrixType>
-struct ei_tridiagonalization_inplace_selector<MatrixType,3>
+struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
 {
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::RealScalar RealScalar;
@@ -470,14 +472,14 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3>
   template<typename DiagonalType, typename SubDiagonalType>
   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
   {
-    diag[0] = ei_real(mat(0,0));
+    diag[0] = mat(0,0);
     RealScalar v1norm2 = ei_abs2(mat(2,0));
-    if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
+    if(v1norm2 == RealScalar(0))
     {
-      diag[1] = ei_real(mat(1,1));
-      diag[2] = ei_real(mat(2,2));
-      subdiag[0] = ei_real(mat(1,0));
-      subdiag[1] = ei_real(mat(2,1));
+      diag[1] = mat(1,1);
+      diag[2] = mat(2,2);
+      subdiag[0] = mat(1,0);
+      subdiag[1] = mat(2,1);
       if (extractQ)
         mat.setIdentity();
     }
@@ -485,18 +487,18 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3>
     {
       RealScalar beta = ei_sqrt(ei_abs2(mat(1,0)) + v1norm2);
       RealScalar invBeta = RealScalar(1)/beta;
-      Scalar m01 = ei_conj(mat(1,0)) * invBeta;
-      Scalar m02 = ei_conj(mat(2,0)) * invBeta;
-      Scalar q = RealScalar(2)*m01*ei_conj(mat(2,1)) + m02*(mat(2,2) - mat(1,1));
-      diag[1] = ei_real(mat(1,1) + m02*q);
-      diag[2] = ei_real(mat(2,2) - m02*q);
+      Scalar m01 = mat(1,0) * invBeta;
+      Scalar m02 = mat(2,0) * invBeta;
+      Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
+      diag[1] = mat(1,1) + m02*q;
+      diag[2] = mat(2,2) - m02*q;
       subdiag[0] = beta;
-      subdiag[1] = ei_real(ei_conj(mat(2,1)) - m01 * q);
+      subdiag[1] = mat(2,1) - m01 * q;
       if (extractQ)
       {
         mat << 1,   0,    0,
-              0, m01,  m02,
-              0, m02, -m01;
+               0, m01,  m02,
+               0, m02, -m01;
       }
     }
   }
@@ -505,8 +507,8 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3>
 /** \internal
   * Trivial specialization for 1x1 matrices
   */
-template<typename MatrixType>
-struct ei_tridiagonalization_inplace_selector<MatrixType,1>
+template<typename MatrixType, bool IsComplex>
+struct ei_tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
 {
   typedef typename MatrixType::Scalar Scalar;
 
diff --git a/Eigen/src/Geometry/arch/Geometry_SSE.h b/Eigen/src/Geometry/arch/Geometry_SSE.h
index 080b87d4b..0078e4aab 100644
--- a/Eigen/src/Geometry/arch/Geometry_SSE.h
+++ b/Eigen/src/Geometry/arch/Geometry_SSE.h
@@ -64,4 +64,58 @@ struct ei_cross3_impl<Architecture::SSE,VectorLhs,VectorRhs,float,true>
   }
 };
 
+
+
+
+template<class Derived, class OtherDerived>
+struct ei_quat_product<Architecture::SSE, Derived, OtherDerived, double, Aligned>
+{
+  inline static Quaternion<double> run(const QuaternionBase<Derived>& _a, const QuaternionBase<OtherDerived>& _b)
+  {
+  const Packet2d mask = _mm_castsi128_pd(_mm_set_epi32(0x0,0x0,0x80000000,0x0));
+
+  Quaternion<double> res;
+
+  const double* a = _a.coeffs().data();
+  Packet2d b_xy = _b.coeffs().template packet<Aligned>(0);
+  Packet2d b_zw = _b.coeffs().template packet<Aligned>(2);
+  Packet2d a_xx = ei_pset1(a[0]);
+  Packet2d a_yy = ei_pset1(a[1]);
+  Packet2d a_zz = ei_pset1(a[2]);
+  Packet2d a_ww = ei_pset1(a[3]);
+
+  // two temporaries:
+  Packet2d t1, t2;
+
+  /*
+   * t1 = ww*xy + yy*zw
+   * t2 = zz*xy - xx*zw
+   * res.xy = t1 +/- swap(t2)
+   */
+  t1 = ei_padd(ei_pmul(a_ww, b_xy), ei_pmul(a_yy, b_zw));
+  t2 = ei_psub(ei_pmul(a_zz, b_xy), ei_pmul(a_xx, b_zw));
+#ifdef __SSE3__
+  ei_pstore(&res.x(), _mm_addsub_pd(t1, ei_preverse(t2)));
+#else
+  ei_pstore(&res.x(), ei_padd(t1, ei_pxor(mask,ei_preverse(t2))));
+#endif
+  
+  /*
+   * t1 = ww*zw - yy*xy
+   * t2 = zz*zw + xx*xy
+   * res.zw = t1 -/+ swap(t2) = swap( swap(t1) +/- t2)
+   */
+  t1 = ei_psub(ei_pmul(a_ww, b_zw), ei_pmul(a_yy, b_xy));
+  t2 = ei_padd(ei_pmul(a_zz, b_zw), ei_pmul(a_xx, b_xy));
+#ifdef __SSE3__
+  ei_pstore(&res.z(), ei_preverse(_mm_addsub_pd(ei_preverse(t1), t2)));
+#else
+  ei_pstore(&res.z(), ei_psub(t1, ei_pxor(mask,ei_preverse(t2))));
+#endif
+
+  return res;
+}
+};
+
+
 #endif // EIGEN_GEOMETRY_SSE_H
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
index f943a0852..c45e6469d 100644
--- a/Eigen/src/Householder/Householder.h
+++ b/Eigen/src/Householder/Householder.h
@@ -65,7 +65,7 @@ void MatrixBase<Derived>::makeHouseholder(
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
   VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
   
-  RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm();
+  RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
   Scalar c0 = coeff(0);
 
   if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0))
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index d0b70a31c..ea7db9c0f 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -47,6 +47,8 @@ template<typename Derived,
 {
   static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
   {
+    if(Derived::ColsAtCompileTime==Dynamic && m.rows()==0)
+      return typename ei_traits<Derived>::Scalar(1);
     return m.partialPivLu().determinant();
   }
 };
diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h
index 69e1f55ba..620f09289 100644
--- a/Eigen/src/Sparse/DynamicSparseMatrix.h
+++ b/Eigen/src/Sparse/DynamicSparseMatrix.h
@@ -54,7 +54,7 @@ struct ei_traits<DynamicSparseMatrix<_Scalar, _Flags, _Index> >
     ColsAtCompileTime = Dynamic,
     MaxRowsAtCompileTime = Dynamic,
     MaxColsAtCompileTime = Dynamic,
-    Flags = _Flags | NestByRefBit,
+    Flags = _Flags | NestByRefBit | LvalueBit,
     CoeffReadCost = NumTraits<Scalar>::ReadCost,
     SupportedAccessPatterns = OuterRandomAccessPattern
   };
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index ecf8efbe7..820cf2884 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -54,7 +54,7 @@ struct ei_traits<SparseMatrix<_Scalar, _Options, _Index> >
     ColsAtCompileTime = Dynamic,
     MaxRowsAtCompileTime = Dynamic,
     MaxColsAtCompileTime = Dynamic,
-    Flags = _Options | NestByRefBit,
+    Flags = _Options | NestByRefBit | LvalueBit,
     CoeffReadCost = NumTraits<Scalar>::ReadCost,
     SupportedAccessPatterns = InnerRandomAccessPattern
   };
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index dfe79a3b2..c5d0a6981 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/Eigen/src/Sparse/SparseVector.h
@@ -48,7 +48,7 @@ struct ei_traits<SparseVector<_Scalar, _Options, _Index> >
     ColsAtCompileTime = IsColVector ? 1 : Dynamic,
     MaxRowsAtCompileTime = RowsAtCompileTime,
     MaxColsAtCompileTime = ColsAtCompileTime,
-    Flags = _Options | NestByRefBit,
+    Flags = _Options | NestByRefBit | LvalueBit,
     CoeffReadCost = NumTraits<Scalar>::ReadCost,
     SupportedAccessPatterns = InnerRandomAccessPattern
   };
diff --git a/Eigen/src/Sparse/SparseView.h b/Eigen/src/Sparse/SparseView.h
index 3b656d5c8..5a152b255 100644
--- a/Eigen/src/Sparse/SparseView.h
+++ b/Eigen/src/Sparse/SparseView.h
@@ -31,18 +31,23 @@ struct ei_traits<SparseView<MatrixType> > : ei_traits<MatrixType>
 {
   typedef int Index;
   typedef Sparse StorageKind;
+  enum {
+    Flags = int(ei_traits<MatrixType>::Flags) & (RowMajorBit)
+  };
 };
 
 template<typename MatrixType>
 class SparseView : public SparseMatrixBase<SparseView<MatrixType> >
 {
   typedef typename MatrixType::Nested MatrixTypeNested;
+  typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested;
 public:
   EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView)
 
   SparseView(const MatrixType& mat, const Scalar& m_reference = Scalar(0),
              typename NumTraits<Scalar>::Real m_epsilon = NumTraits<Scalar>::dummy_precision()) : 
     m_matrix(mat), m_reference(m_reference), m_epsilon(m_epsilon) {}
+
   class InnerIterator;
 
   inline Index rows() const { return m_matrix.rows(); }
@@ -58,10 +63,10 @@ protected:
 };
 
 template<typename MatrixType>
-class SparseView<MatrixType>::InnerIterator : public MatrixTypeNested::InnerIterator
+class SparseView<MatrixType>::InnerIterator : public _MatrixTypeNested::InnerIterator
 {
 public:
-  typedef typename MatrixTypeNested::InnerIterator IterBase;
+  typedef typename _MatrixTypeNested::InnerIterator IterBase;
   InnerIterator(const SparseView& view, Index outer) :
   IterBase(view.m_matrix, outer), m_view(view)
   {