sync with default branch

2025-02-23 18:20:47 +08:00 · 2010-07-22 16:29:35 +02:00 · 2010-07-22 16:29:35 +02:00 · 7020f30da3
commit 7020f30da3
parent b9edd6fb85 96ba7cd655
91 changed files with 5391 additions and 602 deletions
--- 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
  };
 };

--- 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 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 top columns of *this.
+
+
+/** \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.
  *
--- a/Eigen/src/Core/CwiseNullaryOp.h
+++ b/Eigen/src/Core/CwiseNullaryOp.h
@ -240,16 +240,29 @@ 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.
  *
@ -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));
--- 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:
--- 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();
--- 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();
--- 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
    }
--- 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,
--- 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
  };
 };

--- 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;
 };
--- 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>());
 }

--- 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());
--- 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;
 }
--- 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
--- 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)
--- 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);
--- 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,
--- 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,
--- 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();
 }

 /*****************************************************************************
--- 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,
--- 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>
--- 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;
--- 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(); }
--- 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;

--- 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
--- 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))
--- 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();
  }
 };
--- 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
  };
--- 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
  };
--- 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
  };
--- 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)
  {
--- a/bench/bench_gemm.cpp
+++ b/bench/bench_gemm.cpp
@ -10,8 +10,8 @@ using namespace std;
 using namespace Eigen;

 #ifndef SCALAR
-#define SCALAR std::complex<float>
-//  #define SCALAR float
+// #define SCALAR std::complex<float>
+#define SCALAR float
 #endif

 typedef SCALAR Scalar;
--- a/bench/btl/libs/C_BLAS/blas.h
+++ b/bench/btl/libs/C_BLAS/blas.h
@ -11,7 +11,7 @@ typedef long BLASLONG;
 typedef unsigned long BLASULONG;
 #endif

-int    BLASFUNC(xerbla)(char *, int *info, int);
+int    BLASFUNC(xerbla)(const char *, int *info, int);

 float  BLASFUNC(sdot)  (int *, float  *, int *, float  *, int *);
 float  BLASFUNC(sdsdot)(int *, float  *,        float  *, int *, float  *, int *);
@ -38,10 +38,10 @@ void  BLASFUNC(zdotc)  (double *, int *, double  *, int *, double  *, int *);
 void  BLASFUNC(xdotu)  (double *, int *, double  *, int *, double  *, int *);
 void  BLASFUNC(xdotc)  (double *, int *, double  *, int *, double  *, int *);
 #else
-float   BLASFUNC(cdotu)  (int *, float  *, int *, float  *, int *);
-float   BLASFUNC(cdotc)  (int *, float  *, int *, float  *, int *);
-double  BLASFUNC(zdotu)  (int *, double  *, int *, double  *, int *);
-double  BLASFUNC(zdotc)  (int *, double  *, int *, double  *, int *);
+std::complex<float>   BLASFUNC(cdotu)  (int *, float  *, int *, float  *, int *);
+std::complex<float>   BLASFUNC(cdotc)  (int *, float  *, int *, float  *, int *);
+std::complex<double>  BLASFUNC(zdotu)  (int *, double  *, int *, double  *, int *);
+std::complex<double>  BLASFUNC(zdotc)  (int *, double  *, int *, double  *, int *);
 double  BLASFUNC(xdotu)  (int *, double  *, int *, double  *, int *);
 double  BLASFUNC(xdotc)  (int *, double  *, int *, double  *, int *);
 #endif
--- a/bench/eig33.cpp
+++ b/bench/eig33.cpp
@ -0,0 +1,139 @@
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+#include <Eigen/Geometry>
+#include <bench/BenchTimer.h>
+
+using namespace Eigen;
+using namespace std;
+
+template<typename Matrix, typename Roots>
+inline void computeRoots (const Matrix& rkA, Roots& adRoot)
+{
+  typedef typename Matrix::Scalar Scalar;
+  const Scalar msInv3 = 1.0/3.0;
+  const Scalar msRoot3 = ei_sqrt(Scalar(3.0));
+
+  Scalar dA00 = rkA(0,0);
+  Scalar dA01 = rkA(0,1);
+  Scalar dA02 = rkA(0,2);
+  Scalar dA11 = rkA(1,1);
+  Scalar dA12 = rkA(1,2);
+  Scalar dA22 = rkA(2,2);
+
+  // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
+  // eigenvalues are the roots to this equation, all guaranteed to be
+  // real-valued, because the matrix is symmetric.
+  Scalar dC0 = dA00*dA11*dA22 + Scalar(2)*dA01*dA02*dA12 - dA00*dA12*dA12 - dA11*dA02*dA02 - dA22*dA01*dA01;
+  Scalar dC1 = dA00*dA11 - dA01*dA01 + dA00*dA22 - dA02*dA02 + dA11*dA22 - dA12*dA12;
+  Scalar dC2 = dA00 + dA11 + dA22;
+
+  // Construct the parameters used in classifying the roots of the equation
+  // and in solving the equation for the roots in closed form.
+  Scalar dC2Div3 = dC2*msInv3;
+  Scalar dADiv3 = (dC1 - dC2*dC2Div3)*msInv3;
+  if (dADiv3 > Scalar(0))
+    dADiv3 = Scalar(0);
+
+  Scalar dMBDiv2 = Scalar(0.5)*(dC0 + dC2Div3*(Scalar(2)*dC2Div3*dC2Div3 - dC1));
+
+  Scalar dQ = dMBDiv2*dMBDiv2 + dADiv3*dADiv3*dADiv3;
+  if (dQ > Scalar(0))
+    dQ = Scalar(0);
+
+  // Compute the eigenvalues by solving for the roots of the polynomial.
+  Scalar dMagnitude = ei_sqrt(-dADiv3);
+  Scalar dAngle = std::atan2(ei_sqrt(-dQ),dMBDiv2)*msInv3;
+  Scalar dCos = ei_cos(dAngle);
+  Scalar dSin = ei_sin(dAngle);
+  adRoot(0) = dC2Div3 + 2.f*dMagnitude*dCos;
+  adRoot(1) = dC2Div3 - dMagnitude*(dCos + msRoot3*dSin);
+  adRoot(2) = dC2Div3 - dMagnitude*(dCos - msRoot3*dSin);
+
+  // Sort in increasing order.
+  if (adRoot(0) >= adRoot(1))
+    std::swap(adRoot(0),adRoot(1));
+  if (adRoot(1) >= adRoot(2))
+  {
+    std::swap(adRoot(1),adRoot(2));
+    if (adRoot(0) >= adRoot(1))
+      std::swap(adRoot(0),adRoot(1));
+  }
+}
+
+template<typename Matrix, typename Vector>
+void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
+{
+  typedef typename Matrix::Scalar Scalar;
+    // Scale the matrix so its entries are in [-1,1].  The scaling is applied
+    // only when at least one matrix entry has magnitude larger than 1.
+
+    Scalar scale = mat.cwiseAbs()/*.template triangularView<Lower>()*/.maxCoeff();
+    scale = std::max(scale,Scalar(1));
+    Matrix scaledMat = mat / scale;
+
+    // Compute the eigenvalues
+//     scaledMat.setZero();
+    computeRoots(scaledMat,evals);
+
+    // compute the eigen vectors
+    // here we assume 3 differents eigenvalues
+
+    // "optimized version" which appears to be slower with gcc!
+//     Vector base;
+//     Scalar alpha, beta;
+//     base <<   scaledMat(1,0) * scaledMat(2,1),
+//               scaledMat(1,0) * scaledMat(2,0),
+//              -scaledMat(1,0) * scaledMat(1,0);
+//     for(int k=0; k<2; ++k)
+//     {
+//       alpha = scaledMat(0,0) - evals(k);
+//       beta  = scaledMat(1,1) - evals(k);
+//       evecs.col(k) = (base + Vector(-beta*scaledMat(2,0), -alpha*scaledMat(2,1), alpha*beta)).normalized();
+//     }
+//     evecs.col(2) = evecs.col(0).cross(evecs.col(1)).normalized();
+
+    // naive version
+    Matrix tmp;
+    tmp = scaledMat;
+    tmp.diagonal().array() -= evals(0);
+    evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized();
+
+    tmp = scaledMat;
+    tmp.diagonal().array() -= evals(1);
+    evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized();
+
+    tmp = scaledMat;
+    tmp.diagonal().array() -= evals(2);
+    evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
+    
+    // Rescale back to the original size.
+    evals *= scale;
+}
+
+int main()
+{
+  BenchTimer t;
+  int tries = 10;
+  int rep = 400000;
+  typedef Matrix3f Mat;
+  typedef Vector3f Vec;
+  Mat A = Mat::Random(3,3);
+  A = A.adjoint() * A;
+
+  SelfAdjointEigenSolver<Mat> eig(A);
+  BENCH(t, tries, rep, eig.compute(A));
+  std::cout << "Eigen:  " << t.best() << "s\n";
+
+  Mat evecs;
+  Vec evals;
+  BENCH(t, tries, rep, eigen33(A,evecs,evals));
+  std::cout << "Direct: " << t.best() << "s\n\n";
+
+  std::cerr << "Eigenvalue/eigenvector diffs:\n";
+  std::cerr << (evals - eig.eigenvalues()).transpose() << "\n";
+  for(int k=0;k<3;++k)
+    if(evecs.col(k).dot(eig.eigenvectors().col(k))<0)
+      evecs.col(k) = -evecs.col(k);
+  std::cerr << evecs - eig.eigenvectors() << "\n\n";
+}
--- a/bench/quatmul.cpp
+++ b/bench/quatmul.cpp
@ -0,0 +1,47 @@
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <bench/BenchTimer.h>
+
+using namespace Eigen; 
+
+template<typename Quat>
+EIGEN_DONT_INLINE void quatmul_default(const Quat& a, const Quat& b, Quat& c)
+{
+  c = a * b;
+}
+
+template<typename Quat>
+EIGEN_DONT_INLINE void quatmul_novec(const Quat& a, const Quat& b, Quat& c)
+{
+  c = ei_quat_product<0, Quat, Quat, typename Quat::Scalar, Aligned>::run(a,b);
+}
+
+template<typename Quat> void bench(const std::string& label)
+{
+  int tries = 10;
+  int rep = 1000000;
+  BenchTimer t;
+  
+  Quat a(4, 1, 2, 3);
+  Quat b(2, 3, 4, 5);
+  Quat c;
+  
+  std::cout.precision(3);
+  
+  BENCH(t, tries, rep, quatmul_default(a,b,c));
+  std::cout << label << " default " << 1e3*t.best(CPU_TIMER) << "ms  \t" << 1e-6*double(rep)/(t.best(CPU_TIMER)) << " M mul/s\n";
+  
+  BENCH(t, tries, rep, quatmul_novec(a,b,c));
+  std::cout << label << " novec   " << 1e3*t.best(CPU_TIMER) << "ms  \t" << 1e-6*double(rep)/(t.best(CPU_TIMER)) << " M mul/s\n";
+}
+
+int main()
+{
+  bench<Quaternionf>("float ");
+  bench<Quaterniond>("double");
+
+  return 0;
+
+}
+
--- a/blas/common.h
+++ b/blas/common.h
@ -26,6 +26,7 @@
 #define EIGEN_BLAS_COMMON_H

 #include <iostream>
+#include <complex>

 #ifndef SCALAR
 #error the token SCALAR must be defined to compile this file
@ -86,15 +87,15 @@ enum
  Conj = IsComplex
 };

-typedef Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<Dynamic> > MatrixType;
+typedef Map<Matrix<Scalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > MatrixType;
 typedef Map<Matrix<Scalar,Dynamic,1>, 0, InnerStride<Dynamic> > StridedVectorType;
 typedef Map<Matrix<Scalar,Dynamic,1> > CompactVectorType;

 template<typename T>
-Map<Matrix<T,Dynamic,Dynamic>, 0, OuterStride<Dynamic> >
+Map<Matrix<T,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> >
 matrix(T* data, int rows, int cols, int stride)
 {
-  return Map<Matrix<T,Dynamic,Dynamic>, 0, OuterStride<Dynamic> >(data, rows, cols, OuterStride<Dynamic>(stride));
+  return Map<Matrix<T,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> >(data, rows, cols, OuterStride<>(stride));
 }

 template<typename T>
--- a/blas/level1_impl.h
+++ b/blas/level1_impl.h
@ -30,8 +30,6 @@ int EIGEN_BLAS_FUNC(axpy)(int *n, RealScalar *palpha, RealScalar *px, int *incx,
  Scalar* y = reinterpret_cast<Scalar*>(py);
  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);

-//   std::cerr << "axpy " << *n << " " << alpha << " " << *incx << " " << *incy << "\n";
-
  if(*incx==1 && *incy==1)    vector(y,*n) += alpha * vector(x,*n);
  else if(*incx>0 && *incy>0) vector(y,*n,*incy) += alpha * vector(x,*n,*incx);
  else if(*incx>0 && *incy<0) vector(y,*n,-*incy).reverse() += alpha * vector(x,*n,*incx);
@ -126,7 +124,7 @@ int EIGEN_CAT(EIGEN_CAT(i,SCALAR_SUFFIX),amax_)(int *n, RealScalar *px, int *inc
  if(*n<=0)
    return 0;

-  int ret;
+  DenseIndex ret;

  if(*incx==1)  vector(x,*n).cwiseAbs().maxCoeff(&ret);
  else          vector(x,*n,std::abs(*incx)).cwiseAbs().maxCoeff(&ret);
@ -155,44 +153,36 @@ Scalar EIGEN_BLAS_FUNC(sdot)(int *n, RealScalar *px, int *incx, RealScalar *py,
 #if ISCOMPLEX

 // computes a dot product of a conjugated vector with another vector.
-void EIGEN_BLAS_FUNC(dotc)(RealScalar* dot, int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+Scalar EIGEN_BLAS_FUNC(dotc)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
 {
-
-  std::cerr << "Eigen BLAS: _dotc is not implemented yet\n";
-
-  return;
-
-  // TODO: find how to return a complex to fortran
-
 //   std::cerr << "_dotc " << *n << " " << *incx << " " << *incy << "\n";

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);

-  if(*incx==1 && *incy==1)
-    *reinterpret_cast<Scalar*>(dot) = vector(x,*n).dot(vector(y,*n));
-  else
-    *reinterpret_cast<Scalar*>(dot) = vector(x,*n,*incx).dot(vector(y,*n,*incy));
+  Scalar res;
+  if(*incx==1 && *incy==1)    res = (vector(x,*n).dot(vector(y,*n)));
+  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).dot(vector(y,*n,*incy)));
+  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,*incy)));
+  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).dot(vector(y,*n,-*incy).reverse()));
+  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().dot(vector(y,*n,-*incy).reverse()));
+  return res;
 }

 // computes a vector-vector dot product without complex conjugation.
-void EIGEN_BLAS_FUNC(dotu)(RealScalar* dot, int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
+Scalar EIGEN_BLAS_FUNC(dotu)(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
 {
-  std::cerr << "Eigen BLAS: _dotu is not implemented yet\n";
-
-  return;
-
-  // TODO: find how to return a complex to fortran
-
 //   std::cerr << "_dotu " << *n << " " << *incx << " " << *incy << "\n";

  Scalar* x = reinterpret_cast<Scalar*>(px);
  Scalar* y = reinterpret_cast<Scalar*>(py);
-
-  if(*incx==1 && *incy==1)
-    *reinterpret_cast<Scalar*>(dot) = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
-  else
-    *reinterpret_cast<Scalar*>(dot) = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
+  Scalar res;
+  if(*incx==1 && *incy==1)    res = (vector(x,*n).cwiseProduct(vector(y,*n))).sum();
+  else if(*incx>0 && *incy>0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx<0 && *incy>0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,*incy))).sum();
+  else if(*incx>0 && *incy<0) res = (vector(x,*n,*incx).cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  else if(*incx<0 && *incy<0) res = (vector(x,*n,-*incx).reverse().cwiseProduct(vector(y,*n,-*incy).reverse())).sum();
+  return res;
 }

 #endif // ISCOMPLEX
@ -253,15 +243,60 @@ int EIGEN_BLAS_FUNC(rot)(int *n, RealScalar *px, int *incx, RealScalar *py, int

 int EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar *ps)
 {
-  Scalar a = *reinterpret_cast<Scalar*>(pa);
-  Scalar b = *reinterpret_cast<Scalar*>(pb);
-  Scalar* c = reinterpret_cast<Scalar*>(pc);
+  Scalar& a = *reinterpret_cast<Scalar*>(pa);
+  Scalar& b = *reinterpret_cast<Scalar*>(pb);
+  RealScalar* c = pc;
  Scalar* s = reinterpret_cast<Scalar*>(ps);

-  PlanarRotation<Scalar> r;
-  r.makeGivens(a,b);
-  *c = r.c();
-  *s = r.s();
+  #if !ISCOMPLEX
+  Scalar r,z;
+  Scalar aa = ei_abs(a);
+  Scalar ab = ei_abs(b);
+  if((aa+ab)==Scalar(0))
+  {
+    *c = 1;
+    *s = 0;
+    r = 0;
+    z = 0;
+  }
+  else
+  {
+    r = ei_sqrt(a*a + b*b);
+    Scalar amax = aa>ab ? a : b;
+    r = amax>0 ? r : -r;
+    *c = a/r;
+    *s = b/r;
+    z = 1;
+    if (aa > ab) z = *s;
+    if (ab > aa && *c!=RealScalar(0))
+      z = Scalar(1)/ *c;
+  }
+  *pa = r;
+  *pb = z;
+  #else
+  Scalar alpha;
+  RealScalar norm,scale;
+  if(ei_abs(a)==RealScalar(0))
+  {
+    *c = RealScalar(0);
+    *s = Scalar(1);
+    a = b;
+  }
+  else
+  {
+    scale = ei_abs(a) + ei_abs(b);
+    norm = scale*ei_sqrt((ei_abs2(a/scale))+ (ei_abs2(b/scale)));
+    alpha = a/ei_abs(a);
+    *c = ei_abs(a)/norm;
+    *s = alpha*ei_conj(b)/norm;
+    a = alpha*norm;
+  }
+  #endif
+
+//   PlanarRotation<Scalar> r;
+//   r.makeGivens(a,b);
+//   *c = r.c();
+//   *s = r.s();

  return 0;
 }
--- a/blas/level3_impl.h
+++ b/blas/level3_impl.h
@ -27,7 +27,7 @@
 int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
 {
 //   std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
-  typedef void (*functype)(int, int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar, Eigen::GemmParallelInfo<Scalar>*);
+  typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, ei_level3_blocking<Scalar,Scalar>&, Eigen::GemmParallelInfo<DenseIndex>*);
  static functype func[12];

  static bool init = false;
@ -35,15 +35,15 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScal
  {
    for(int k=0; k<12; ++k)
      func[k] = 0;
-    func[NOTR  | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,ColMajor,false,ColMajor>::run);
-    func[ADJ   | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, ColMajor,false,ColMajor>::run);
-    func[NOTR  | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,false,ColMajor>::run);
-    func[TR    | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,false,ColMajor>::run);
-    func[ADJ   | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,false,ColMajor>::run);
-    func[NOTR  | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,ColMajor,false,RowMajor,Conj, ColMajor>::run);
-    func[TR    | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,false,RowMajor,Conj, ColMajor>::run);
-    func[ADJ   | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,RowMajor,Conj, RowMajor,Conj, ColMajor>::run);
+    func[NOTR  | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,false,ColMajor,false,ColMajor>::run);
+    func[ADJ   | (NOTR << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,Conj, ColMajor,false,ColMajor>::run);
+    func[NOTR  | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,ColMajor,false,RowMajor,false,ColMajor>::run);
+    func[TR    | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,false,RowMajor,false,ColMajor>::run);
+    func[ADJ   | (TR   << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,Conj, RowMajor,false,ColMajor>::run);
+    func[NOTR  | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,ColMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[TR    | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[ADJ   | (ADJ  << 2)] = (ei_general_matrix_matrix_product<Scalar,DenseIndex,RowMajor,Conj, RowMajor,Conj, ColMajor>::run);
    init = true;
  }

@ -62,19 +62,21 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScal
  }

  if(beta!=Scalar(1))
-    if(beta==Scalar(0))
-      matrix(c, *m, *n, *ldc).setZero();
-    else
-      matrix(c, *m, *n, *ldc) *= beta;
+  {
+    if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
+    else                matrix(c, *m, *n, *ldc) *= beta;
+  }

-  func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, 0);
+  ei_gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k);
+
+  func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0);
  return 0;
 }

 int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
 {
 //   std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
-  typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int);
+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex);
  static functype func[32];

  static bool init = false;
@ -83,38 +85,38 @@ int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m,
    for(int k=0; k<32; ++k)
      func[k] = 0;

-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|0,          false,ColMajor,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|0,          false,RowMajor,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,ColMajor,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,RowMajor,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|0,          false,ColMajor,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|0,          false,RowMajor,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,ColMajor,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,RowMajor,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|0,          false,ColMajor,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|0,          false,RowMajor,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,ColMajor,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,RowMajor,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|0,          false,ColMajor,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|0,          false,RowMajor,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,ColMajor,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,RowMajor,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor>::run);


-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);

    init = true;
  }
@ -148,7 +150,7 @@ int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m,
 int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
 {
 //   std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n";
-  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
+  typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
  static functype func[32];
  static bool init = false;
  if(!init)
@ -156,37 +158,37 @@ int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m,
    for(int k=0; k<32; ++k)
      func[k] = 0;

-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);

-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);

-    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);

-    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
-    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
+    func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
+    func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);

-    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
-    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
-    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
+    func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
+    func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
+    func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (ei_product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);

    init = true;
  }
@ -203,14 +205,17 @@ int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m,
    return 0;
  }

+  if(*m==0 || *n==0)
+    return 1;
+
  // FIXME find a way to avoid this copy
-  Matrix<Scalar,Dynamic,Dynamic> tmp = matrix(b,*m,*n,*ldb);
+  Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb);
  matrix(b,*m,*n,*ldb).setZero();

  if(SIDE(*side)==LEFT)
-    func[code](*m, *n, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha);
+    func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha);
  else
-    func[code](*n, *m, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha);
+    func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha);
  return 1;
 }

@ -233,19 +238,46 @@ int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *pa
  }

  if(beta!=Scalar(1))
+  {
    if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
    else                matrix(c, *m, *n, *ldc) *= beta;
+  }

+  if(*m==0 || *n==0)
+  {
+    return 1;
+  }
+
+  #if ISCOMPLEX
+  // FIXME add support for symmetric complex matrix
+  int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
+  Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size);
+  if(UPLO(*uplo)==UP)
+  {
+    matA.triangularView<Upper>() = matrix(a,size,size,*lda);
+    matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose();
+  }
+  else if(UPLO(*uplo)==LO)
+  {
+    matA.triangularView<Lower>() = matrix(a,size,size,*lda);
+    matA.triangularView<Upper>() = matrix(a,size,size,*lda).transpose();
+  }
  if(SIDE(*side)==LEFT)
-    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
-    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
+    matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb);
+  else if(SIDE(*side)==RIGHT)
+    matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA;
+  #else
+  if(SIDE(*side)==LEFT)
+    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
+    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
    else                      return 0;
  else if(SIDE(*side)==RIGHT)
-    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
-    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
+    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
+    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
    else                      return 0;
  else
    return 0;
+  #endif

  return 0;
 }
@ -255,7 +287,7 @@ int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *pa
 int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
 {
 //   std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
-  typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar);
+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
  static functype func[8];

  static bool init = false;
@ -264,13 +296,13 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
    for(int k=0; k<8; ++k)
      func[k] = 0;

-    func[NOTR  | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, Upper>::run);
-    func[TR    | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Upper>::run);
-    func[ADJ   | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Upper>::run);
+    func[NOTR  | (UP << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,ColMajor,ColMajor,true, Upper>::run);
+    func[TR    | (UP << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Upper>::run);
+    func[ADJ   | (UP << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Upper>::run);

-    func[NOTR  | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, Lower>::run);
-    func[TR    | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Lower>::run);
-    func[ADJ   | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Lower>::run);
+    func[NOTR  | (LO << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,ColMajor,ColMajor,true, Lower>::run);
+    func[TR    | (LO << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Lower>::run);
+    func[ADJ   | (LO << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Lower>::run);

    init = true;
  }
@ -289,10 +321,30 @@ int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
  }

  if(beta!=Scalar(1))
+  {
    if(UPLO(*uplo)==UP) matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
    else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
+  }

+  #if ISCOMPLEX
+  // FIXME add support for symmetric complex matrix
+  if(UPLO(*uplo)==UP)
+  {
+    if(OP(*op)==NOTR)
+      matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
+    else
+      matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
+  }
+  else
+  {
+    if(OP(*op)==NOTR)
+      matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
+    else
+      matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
+  }
+  #else
  func[code](*n, *k, a, *lda, c, *ldc, alpha);
+  #endif

  return 0;
 }
@ -307,8 +359,44 @@ int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *pal
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);

-  // TODO
-  std::cerr << "Eigen BLAS: _syr2k is not implemented yet\n";
+  if(*n<=0 || *k<0)
+  {
+    return 0;
+  }
+
+  if(beta!=Scalar(1))
+  {
+    if(UPLO(*uplo)==UP) matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
+    else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
+  }
+
+  if(*k==0)
+    return 1;
+
+  if(OP(*op)==NOTR)
+  {
+    if(UPLO(*uplo)==UP)
+    {
+      matrix(c, *n, *n, *ldc).triangularView<Upper>()
+        += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
+        +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
+    }
+    else if(UPLO(*uplo)==LO)
+      matrix(c, *n, *n, *ldc).triangularView<Lower>()
+        += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
+        +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
+  }
+  else if(OP(*op)==TR || OP(*op)==ADJ)
+  {
+    if(UPLO(*uplo)==UP)
+      matrix(c, *n, *n, *ldc).triangularView<Upper>()
+        += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
+        +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
+    else if(UPLO(*uplo)==LO)
+      matrix(c, *n, *n, *ldc).triangularView<Lower>()
+        += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
+        +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
+  }

  return 0;
 }
@ -333,17 +421,32 @@ int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *pa
    return 0;
  }

-  if(beta!=Scalar(1))
+  if(beta==Scalar(0))
+    matrix(c, *m, *n, *ldc).setZero();
+  else if(beta!=Scalar(1))
    matrix(c, *m, *n, *ldc) *= beta;
+  
+  if(*m==0 || *n==0)
+  {
+    return 1;
+  }

  if(SIDE(*side)==LEFT)
-    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, RowMajor,true,Conj,  ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
-    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
+  {
+    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj,  ColMajor,false,false, ColMajor>
+                                ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
+    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor>
+                                ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
    else                      return 0;
+  }
  else if(SIDE(*side)==RIGHT)
-    if(UPLO(*uplo)==UP)       ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, RowMajor,true,Conj,  ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
-    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
+  {
+    if(UPLO(*uplo)==UP)       matrix(c,*m,*n,*ldc) += alpha * matrix(b,*m,*n,*ldb) * matrix(a,*n,*n,*lda).selfadjointView<Upper>();/*ei_product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, RowMajor,true,Conj,  ColMajor>
+                                ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/
+    else if(UPLO(*uplo)==LO)  ei_product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor>
+                                ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
    else                      return 0;
+  }
  else
  {
    return 0;
@ -356,7 +459,7 @@ int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *pa
 // c = alpha*conj(a')*a + beta*c  for op  = 'C'or'c'
 int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
 {
-  typedef void (*functype)(int, int, const Scalar *, int, Scalar *, int, Scalar);
+  typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
  static functype func[8];

  static bool init = false;
@ -365,11 +468,11 @@ int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
    for(int k=0; k<8; ++k)
      func[k] = 0;

-    func[NOTR  | (UP << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, Upper>::run);
-    func[ADJ   | (UP << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Upper>::run);
+    func[NOTR  | (UP << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,ColMajor,ColMajor,true, Upper>::run);
+    func[ADJ   | (UP << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Upper>::run);

-    func[NOTR  | (LO << 2)] = (ei_selfadjoint_product<Scalar,ColMajor,ColMajor,true, Lower>::run);
-    func[ADJ   | (LO << 2)] = (ei_selfadjoint_product<Scalar,RowMajor,ColMajor,false,Lower>::run);
+    func[NOTR  | (LO << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,ColMajor,ColMajor,true, Lower>::run);
+    func[ADJ   | (LO << 2)] = (ei_selfadjoint_product<Scalar,DenseIndex,RowMajor,ColMajor,false,Lower>::run);

    init = true;
  }
@ -408,24 +511,60 @@ int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palp
 }

 // c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c,  for op = 'N'or'n'
-// c = alpha*conj(b')*a + conj(alpha)*conj(a')*b + beta*c,  for op = 'C'or'c'
+// c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c,  for op = 'C'or'c'
 int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
 {
  Scalar* a = reinterpret_cast<Scalar*>(pa);
  Scalar* b = reinterpret_cast<Scalar*>(pb);
  Scalar* c = reinterpret_cast<Scalar*>(pc);
  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
-  Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
+  RealScalar beta  = *pbeta;

-  if(*n<0 || *k<0)
+  if(*n<=0 || *k<0)
  {
    return 0;
  }

-  // TODO
-  std::cerr << "Eigen BLAS: _her2k is not implemented yet\n";
+  if(beta!=RealScalar(1))
+  {
+    if(UPLO(*uplo)==UP) matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
+    else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;

-  return 0;
+    matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
+    matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
+  }
+  else if(*k>0 && alpha!=Scalar(0))
+    matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
+
+  if(*k==0)
+    return 1;
+
+  if(OP(*op)==NOTR)
+  {
+    if(UPLO(*uplo)==UP)
+    {
+      matrix(c, *n, *n, *ldc).triangularView<Upper>()
+        +=         alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
+        +  ei_conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
+    }
+    else if(UPLO(*uplo)==LO)
+      matrix(c, *n, *n, *ldc).triangularView<Lower>()
+        += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
+        +  ei_conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
+  }
+  else if(OP(*op)==ADJ)
+  {
+    if(UPLO(*uplo)==UP)
+      matrix(c, *n, *n, *ldc).triangularView<Upper>()
+        += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
+        +  ei_conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
+    else if(UPLO(*uplo)==LO)
+      matrix(c, *n, *n, *ldc).triangularView<Lower>()
+        += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
+        +  ei_conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
+  }
+
+  return 1;
 }

 #endif // ISCOMPLEX
--- a/blas/xerbla.cpp
+++ b/blas/xerbla.cpp
@ -6,7 +6,7 @@ extern "C"
 {
 #endif

-int xerbla_(char * msg, int *info, int)
+int xerbla_(const char * msg, int *info, int)
 {
  std::cerr << "Eigen BLAS ERROR #" << *info << ": " << msg << "\n";
  return 0;
--- a/cmake/FindGMP.cmake
+++ b/cmake/FindGMP.cmake
@ -0,0 +1,21 @@
+# Try to find the GNU Multiple Precision Arithmetic Library (GMP)
+# See http://gmplib.org/
+
+if (GMP_INCLUDES AND GMP_LIBRARIES)
+  set(GMP_FIND_QUIETLY TRUE)
+endif (GMP_INCLUDES AND GMP_LIBRARIES)
+
+find_path(GMP_INCLUDES
+  NAMES
+  gmp.h
+  PATHS
+  $ENV{GMPDIR}
+  ${INCLUDE_INSTALL_DIR}
+)
+
+find_library(GMP_LIBRARIES gmp PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR})
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(GMP DEFAULT_MSG
+                                  GMP_INCLUDES GMP_LIBRARIES)
+mark_as_advanced(GMP_INCLUDES GMP_LIBRARIES)
--- a/cmake/FindMPFR.cmake
+++ b/cmake/FindMPFR.cmake
@ -0,0 +1,83 @@
+# Try to find the MPFR library
+# See http://www.mpfr.org/
+#
+# This module supports requiring a minimum version, e.g. you can do
+#   find_package(MPFR 2.3.0)
+# to require version 2.3.0 to newer of MPFR.
+#
+# Once done this will define
+#
+#  MPFR_FOUND - system has MPFR lib with correct version
+#  MPFR_INCLUDES - the MPFR include directory
+#  MPFR_LIBRARIES - the MPFR library
+#  MPFR_VERSION - MPFR version
+
+# Copyright (c) 2006, 2007 Montel Laurent, <montel@kde.org>
+# Copyright (c) 2008, 2009 Gael Guennebaud, <g.gael@free.fr>
+# Copyright (c) 2010 Jitse Niesen, <jitse@maths.leeds.ac.uk>
+# Redistribution and use is allowed according to the terms of the BSD license.
+
+# Set MPFR_INCLUDES
+
+find_path(MPFR_INCLUDES
+  NAMES
+  mpfr.h
+  PATHS
+  $ENV{GMPDIR}
+  ${INCLUDE_INSTALL_DIR}
+)
+
+# Set MPFR_FIND_VERSION to 1.0.0 if no minimum version is specified
+
+if(NOT MPFR_FIND_VERSION)
+  if(NOT MPFR_FIND_VERSION_MAJOR)
+    set(MPFR_FIND_VERSION_MAJOR 1)
+  endif(NOT MPFR_FIND_VERSION_MAJOR)
+  if(NOT MPFR_FIND_VERSION_MINOR)
+    set(MPFR_FIND_VERSION_MINOR 0)
+  endif(NOT MPFR_FIND_VERSION_MINOR)
+  if(NOT MPFR_FIND_VERSION_PATCH)
+    set(MPFR_FIND_VERSION_PATCH 0)
+  endif(NOT MPFR_FIND_VERSION_PATCH)
+
+  set(MPFR_FIND_VERSION "${MPFR_FIND_VERSION_MAJOR}.${MPFR_FIND_VERSION_MINOR}.${MPFR_FIND_VERSION_PATCH}")
+endif(NOT MPFR_FIND_VERSION)
+
+
+if(MPFR_INCLUDES)
+
+  # Set MPFR_VERSION
+  
+  file(READ "${MPFR_INCLUDES}/mpfr.h" _mpfr_version_header)
+  
+  string(REGEX MATCH "define[ \t]+MPFR_VERSION_MAJOR[ \t]+([0-9]+)" _mpfr_major_version_match "${_mpfr_version_header}")
+  set(MPFR_MAJOR_VERSION "${CMAKE_MATCH_1}")
+  string(REGEX MATCH "define[ \t]+MPFR_VERSION_MINOR[ \t]+([0-9]+)" _mpfr_minor_version_match "${_mpfr_version_header}")
+  set(MPFR_MINOR_VERSION "${CMAKE_MATCH_1}")
+  string(REGEX MATCH "define[ \t]+MPFR_VERSION_PATCHLEVEL[ \t]+([0-9]+)" _mpfr_patchlevel_version_match "${_mpfr_version_header}")
+  set(MPFR_PATCHLEVEL_VERSION "${CMAKE_MATCH_1}")
+  
+  set(MPFR_VERSION ${MPFR_MAJOR_VERSION}.${MPFR_MINOR_VERSION}.${MPFR_PATCHLEVEL_VERSION})
+  
+  # Check whether found version exceeds minimum version
+  
+  if(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})
+    set(MPFR_VERSION_OK FALSE)
+    message(STATUS "MPFR version ${MPFR_VERSION} found in ${MPFR_INCLUDES}, "
+                   "but at least version ${MPFR_FIND_VERSION} is required")
+  else(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})
+    set(MPFR_VERSION_OK TRUE)
+  endif(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})
+
+endif(MPFR_INCLUDES)
+
+# Set MPFR_LIBRARIES
+
+find_library(MPFR_LIBRARIES mpfr PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR})
+
+# Epilogue
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(MPFR DEFAULT_MSG
+                                  MPFR_INCLUDES MPFR_LIBRARIES MPFR_VERSION_OK)
+mark_as_advanced(MPFR_INCLUDES MPFR_LIBRARIES)
--- a/demos/mandelbrot/mandelbrot.cpp
+++ b/demos/mandelbrot/mandelbrot.cpp
@ -23,6 +23,7 @@
 // Eigen. If not, see <http://www.gnu.org/licenses/>.

 #include "mandelbrot.h"
+#include <iostream>
 #include<QtGui/QPainter>
 #include<QtGui/QImage>
 #include<QtGui/QMouseEvent>
@ -45,7 +46,7 @@ template<> struct iters_before_test<double> { enum { ret = 16 }; };
 template<typename Real> void MandelbrotThread::render(int img_width, int img_height)
 {
  enum { packetSize = Eigen::ei_packet_traits<Real>::size }; // number of reals in a Packet
-  typedef Eigen::Matrix<Real, packetSize, 1> Packet; // wrap a Packet as a vector
+  typedef Eigen::Array<Real, packetSize, 1> Packet; // wrap a Packet as a vector

  enum { iters_before_test = iters_before_test<Real>::ret };
  max_iter = (max_iter / iters_before_test) * iters_before_test;
@ -54,7 +55,7 @@ template<typename Real> void MandelbrotThread::render(int img_width, int img_hei
  const double xradius = widget->xradius;
  const double yradius = xradius * img_height / img_width;
  const int threadcount = widget->threadcount;
-  typedef Eigen::Matrix<Real, 2, 1> Vector2;
+  typedef Eigen::Array<Real, 2, 1> Vector2;
  Vector2 start(widget->center.x() - widget->xradius, widget->center.y() - yradius);
  Vector2 step(2*widget->xradius/img_width, 2*yradius/img_height);
  total_iter = 0;
@ -87,16 +88,16 @@ template<typename Real> void MandelbrotThread::render(int img_width, int img_hei
        {
 #         define ITERATE \
            pzr_buf = pzr; \
-            pzr = pzr.cwise().square(); \
-            pzr -= pzi.cwise().square(); \
+            pzr = pzr.square(); \
+            pzr -= pzi.square(); \
            pzr += pcr; \
-            pzi = (2*pzr_buf).cwise()*pzi; \
+            pzi = (2*pzr_buf)*pzi; \
            pzi += pci;
          ITERATE ITERATE ITERATE ITERATE
        }
-        pix_dont_diverge = ((pzr.cwise().square() + pzi.cwise().square())
+        pix_dont_diverge = ((pzr.square() + pzi.square())
                           .eval() // temporary fix as what follows is not yet vectorized by Eigen
-                           .cwise() <= Packet::Constant(4))
+                           <= Packet::Constant(4))
                                // the 4 here is not a magic value, it's a math fact that if
                                // the square modulus is >4 then divergence is inevitable.
                           .template cast<int>();
--- a/demos/mandelbrot/mandelbrot.h
+++ b/demos/mandelbrot/mandelbrot.h
@ -25,7 +25,7 @@
 #ifndef MANDELBROT_H
 #define MANDELBROT_H

-#include <Eigen/Array>
+#include <Eigen/Core>
 #include <QtGui/QApplication>
 #include <QtGui/QWidget>
 #include <QtCore/QThread>
--- a/demos/opengl/quaternion_demo.cpp
+++ b/demos/opengl/quaternion_demo.cpp
@ -25,10 +25,11 @@
 #include "quaternion_demo.h"
 #include "icosphere.h"

-#include <Eigen/Array>
+#include <Eigen/Geometry>
 #include <Eigen/QR>
 #include <Eigen/LU>

+#include <iostream>
 #include <QEvent>
 #include <QMouseEvent>
 #include <QInputDialog>
@ -187,8 +188,8 @@ public:

  Matrix3 toRotationMatrix(void) const
  {
-    Vector3 c = m_angles.cwise().cos();
-    Vector3 s = m_angles.cwise().sin();
+    Vector3 c = m_angles.array().cos();
+    Vector3 s = m_angles.array().sin();
    Matrix3 res;
    res <<  c.y()*c.z(),                    -c.y()*s.z(),                   s.y(),
            c.z()*s.x()*s.y()+c.x()*s.z(),  c.x()*c.z()-s.x()*s.y()*s.z(),  -c.y()*s.x(),
--- a/doc/C01_TutorialMatrixClass.dox
+++ b/doc/C01_TutorialMatrixClass.dox
@ -132,13 +132,15 @@ Vector4d c(5.0, 6.0, 7.0, 8.0);
 The primary coefficient accessors and mutators in Eigen are the overloaded parenthesis operators.
 For matrices, the row index is always passed first. For vectors, just pass one index.
 The numbering starts at 0. This example is self-explanatory:
-\include tut_matrix_coefficient_accessors.cpp
-Output: \verbinclude tut_matrix_coefficient_accessors.out

-Note that the syntax
-\code
-m(index)
-\endcode
+<table class="tutorial_code"><tr><td>
+Example: \include tut_matrix_coefficient_accessors.cpp
+</td>
+<td>
+Output: \verbinclude tut_matrix_coefficient_accessors.out
+</td></tr></table>
+
+Note that the syntax <tt> m(index) </tt>
 is not restricted to vectors, it is also available for general matrices, meaning index-based access
 in the array of coefficients. This however depends on the matrix's storage order. All Eigen matrices default to
 column-major storage order, but this can be changed to row-major, see \ref TopicStorageOrders "Storage orders".
@ -149,17 +151,29 @@ would make matrix[i,j] compile to the same thing as matrix[j] !

 \section TutorialMatrixCommaInitializer Comma-initialization

-Matrix and vector coefficients can be conveniently set using the so-called \em comma-initializer syntax.
+%Matrix and vector coefficients can be conveniently set using the so-called \em comma-initializer syntax.
 For now, it is enough to know this example:
-\include Tutorial_commainit_01.cpp
+
+<table class="tutorial_code"><tr><td>
+Example: \include Tutorial_commainit_01.cpp
+</td>
+<td>
 Output: \verbinclude Tutorial_commainit_01.out
+</td></tr></table>
+
+
 The right-hand side can also contain matrix expressions as discussed in \ref TutorialAdvancedInitialization "this page".

 \section TutorialMatrixSizesResizing Resizing

 The current size of a matrix can be retrieved by \link EigenBase::rows() rows()\endlink, \link EigenBase::cols() cols() \endlink and \link EigenBase::size() size()\endlink. These methods return the number of rows, the number of columns and the number of coefficients, respectively. Resizing a dynamic-size matrix is done by the \link DenseStorageBase::resize(Index,Index) resize() \endlink method.
-For example: \include tut_matrix_resize.cpp
+
+<table class="tutorial_code"><tr><td>
+Example: \include tut_matrix_resize.cpp
+</td>
+<td>
 Output: \verbinclude tut_matrix_resize.out
+</td></tr></table>

 The resize() method is a no-operation if the actual matrix size doesn't change; otherwise it is destructive: the values of the coefficients may change.
 If you want a conservative variant of resize() which does not change the coefficients, use \link DenseStorageBase::conservativeResize() conservativeResize()\endlink, see \ref TopicResizing "this page" for more details.
@ -167,14 +181,25 @@ If you want a conservative variant of resize() which does not change the coeffic
 All these methods are still available on fixed-size matrices, for the sake of API uniformity. Of course, you can't actually
 resize a fixed-size matrix. Trying to change a fixed size to an actually different value will trigger an assertion failure;
 but the following code is legal:
-\include tut_matrix_resize_fixed_size.cpp
+
+<table class="tutorial_code"><tr><td>
+Example: \include tut_matrix_resize_fixed_size.cpp
+</td>
+<td>
 Output: \verbinclude tut_matrix_resize_fixed_size.out
+</td></tr></table>
+

 \section TutorialMatrixAssignment Assignment and resizing

 Assignment is the action of copying a matrix into another, using \c operator=. Eigen resizes the matrix on the left-hand side automatically so that it matches the size of the matrix on the right-hand size. For example:
-\include tut_matrix_assignment_resizing.cpp
+
+<table class="tutorial_code"><tr><td>
+Example: \include tut_matrix_assignment_resizing.cpp
+</td>
+<td>
 Output: \verbinclude tut_matrix_assignment_resizing.out
+</td></tr></table>

 Of course, if the left-hand side is of fixed size, resizing it is not allowed.

--- a/doc/C02_TutorialMatrixArithmetic.dox
+++ b/doc/C02_TutorialMatrixArithmetic.dox
@ -43,7 +43,7 @@ also have the same \c Scalar type, as Eigen doesn't do automatic type promotion.
 Example: \include tut_arithmetic_add_sub.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_add_sub.out
+Output: \verbinclude tut_arithmetic_add_sub.out
 </td></tr></table>

 \section TutorialArithmeticScalarMulDiv Scalar multiplication and division
@ -59,7 +59,7 @@ Multiplication and division by a scalar is very simple too. The operators at han
 Example: \include tut_arithmetic_scalar_mul_div.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_scalar_mul_div.out
+Output: \verbinclude tut_arithmetic_scalar_mul_div.out
 </td></tr></table>


@ -93,7 +93,7 @@ The transpose \f$ a^T \f$, conjugate \f$ \bar{a} \f$, and adjoint (i.e., conjuga
 Example: \include tut_arithmetic_transpose_conjugate.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_transpose_conjugate.out
+Output: \verbinclude tut_arithmetic_transpose_conjugate.out
 </td></tr></table>

 For real matrices, \c conjugate() is a no-operation, and so \c adjoint() is 100% equivalent to \c transpose().
@ -103,7 +103,7 @@ As for basic arithmetic operators, \c transpose() and \c adjoint() simply return
 Example: \include tut_arithmetic_transpose_aliasing.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_transpose_aliasing.out
+Output: \verbinclude tut_arithmetic_transpose_aliasing.out
 </td></tr></table>
 This is the so-called \ref TopicAliasing "aliasing issue". In "debug mode", i.e., when \ref TopicAssertions "assertions" have not been disabled, such common pitfalls are automatically detected. 

@ -112,7 +112,7 @@ For \em in-place transposition, as for instance in <tt>a = a.transpose()</tt>, s
 Example: \include tut_arithmetic_transpose_inplace.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_transpose_inplace.out
+Output: \verbinclude tut_arithmetic_transpose_inplace.out
 </td></tr></table>
 There is also the \link MatrixBase::adjointInPlace() adjointInPlace()\endlink function for complex matrices.

@ -129,7 +129,7 @@ two operators:
 Example: \include tut_arithmetic_matrix_mul.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_matrix_mul.out
+Output: \verbinclude tut_arithmetic_matrix_mul.out
 </td></tr></table>

 Note: if you read the above paragraph on expression templates and are worried that doing \c m=m*m might cause
@ -154,7 +154,7 @@ The above-discussed \c operator* cannot be used to compute dot and cross product
 Example: \include tut_arithmetic_dot_cross.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_dot_cross.out
+Output: \verbinclude tut_arithmetic_dot_cross.out
 </td></tr></table>

 Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
@ -168,7 +168,7 @@ Eigen also provides some reduction operations to reduce a given matrix or vector
 Example: \include tut_arithmetic_redux_basic.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_redux_basic.out
+Output: \verbinclude tut_arithmetic_redux_basic.out
 </td></tr></table>

 The \em trace of a matrix, as returned by the function \link MatrixBase::trace() trace()\endlink, is the sum of the diagonal coefficients and can also be computed as efficiently using <tt>a.diagonal().sum()</tt>, as we will see later on.
@ -179,7 +179,7 @@ There also exist variants of the \c minCoeff and \c maxCoeff functions returning
 Example: \include tut_arithmetic_redux_minmax.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_redux_minmax.out
+Output: \verbinclude tut_arithmetic_redux_minmax.out
 </td></tr></table>


--- a/doc/C03_TutorialArrayClass.dox
+++ b/doc/C03_TutorialArrayClass.dox
@ -7,7 +7,7 @@ namespace Eigen {
 \li \b Next: \ref TutorialBlockOperations

 This tutorial aims to provide an overview and explanations on how to use
-Eigen's \link ArrayBase Array \endlink class
+Eigen's Array class.

 \b Table \b of \b contents
  - \ref TutorialArrayClassIntro
@ -15,6 +15,7 @@ Eigen's \link ArrayBase Array \endlink class
  - \ref TutorialArrayClassAccess
  - \ref TutorialArrayClassAddSub
  - \ref TutorialArrayClassMult
+  - \ref TutorialArrayClassCwiseOther
  - \ref TutorialArrayClassConvert

 \section TutorialArrayClassIntro What is the Array class?
@ -27,17 +28,17 @@ such as adding a constant to every coefficient in the array or multiplying two a

 \section TutorialArrayClassTypes Array types
 Array is a class template taking the same template parameters as Matrix.
-As with Matrix, the first 3 template parameters are mandatory:
+As with Matrix, the first three template parameters are mandatory:
 \code
 Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
 \endcode
-And the last 3 template parameters are optional. Since this is exactly the same as for Matrix,
-we won't reexplain it and just refer to \ref TutorialMatrixClass "this page".
+The last three template parameters are optional. Since this is exactly the same as for Matrix,
+we won't explain it again here and just refer to \ref TutorialMatrixClass.

 Eigen also provides typedefs for some common cases, in a way that is similar to the Matrix typedefs
 but with some slight differences, as the word "array" is used for both 1-dimensional and 2-dimensional arrays.
 We adopt that convention that typedefs of the form ArrayNt stand for 1-dimensional arrays, where N and t are
-as in the Matrix typedefs explained on \ref TutorialMatrixClass "this page". For 2-dimensional arrays, we
+the size and the scalar type, as in the Matrix typedefs explained on \ref TutorialMatrixClass "this page". For 2-dimensional arrays, we
 use typedefs of the form ArrayNNt. Some examples are shown in the following table:

 <table class="tutorial_code" align="center">
@ -71,69 +72,103 @@ use typedefs of the form ArrayNNt. Some examples are shown in the following tabl


 \section TutorialArrayClassAccess Accessing values inside an Array
-Write and read access to coefficients of an array expression is done in the same way as with matrix expressions.
-For example:

-\include Tutorial_ArrayClass_accessors.cpp
-Output:
-\verbinclude Tutorial_ArrayClass_accessors.out
+The parenthesis operator is overloaded to provide write and read access to the coefficients of an array, just as with matrices.
+Furthermore, the \c << operator can be used to initialize arrays (via the comma initializer) or to print them.

-For more information about the comma initializer, refer to \ref TutorialAdvancedInitialization "this page".
+<table class="tutorial_code"><tr><td>
+Example: \include Tutorial_ArrayClass_accessors.cpp
+</td>
+<td>
+Output: \verbinclude Tutorial_ArrayClass_accessors.out
+</td></tr></table>
+
+For more information about the comma initializer, see \ref TutorialAdvancedInitialization.


-\section TutorialArrayClassAddSub Addition and substraction
+\section TutorialArrayClassAddSub Addition and subtraction
+
+Adding and subtracting two arrays is the same as for matrices.
+The operation is valid if both arrays have the same size, and the addition or subtraction is done coefficient-wise.
+
+Arrays also support expressions of the form <tt>array + scalar</tt> which add a scalar to each coefficient in the array.
+This provides a functionality that is not directly available for Matrix objects.
+
+<table class="tutorial_code"><tr><td>
+Example: \include Tutorial_ArrayClass_addition.cpp
+</td>
+<td>
+Output: \verbinclude Tutorial_ArrayClass_addition.out
+</td></tr></table>


-Adding and substracting two arrays has the same result as for Matrices.
-This is valid as long as both arrays have the same sizes:
-
-\include Tutorial_ArrayClass_addition.cpp
-Output:
-\verbinclude Tutorial_ArrayClass_addition.out
-
-It is also allowed to add a scalar to each coefficient in an Array,
-providing a functionality that was not available for Matrix objects:
-
-\include Tutorial_ArrayClass_addition_scalar.cpp
-Output:
-\verbinclude Tutorial_ArrayClass_addition_scalar.out
-
-
-\subsection TutorialArrayClassMult Array multiplication
+\section TutorialArrayClassMult Array multiplication

 First of all, of course you can multiply an array by a scalar, this works in the same way as matrices. Where arrays
-are fundamentally different from matrices, is when you multiply two arrays together. While Matrices interprete
-multiplication as matrix product, arrays interprete multiplication as coefficient-wise product. For example:
+are fundamentally different from matrices, is when you multiply two together. Matrices interpret
+multiplication as the matrix product and arrays interpret multiplication as the coefficient-wise product. Thus, two 
+arrays can be multiplied if they have the same size.

-\include Tutorial_ArrayClass_mult.cpp
-Output:
-\verbinclude Tutorial_ArrayClass_mult.out
+<table class="tutorial_code"><tr><td>
+Example: \include Tutorial_ArrayClass_mult.cpp
+</td>
+<td>
+Output: \verbinclude Tutorial_ArrayClass_mult.out
+</td></tr></table>


+\section TutorialArrayClassCwiseOther Other coefficient-wise operations
+
+The Array class defined other coefficient-wise operations besides the addition, subtraction and multiplication
+operators described about. For example, the \link ArrayBase::abs() .abs() \endlink method takes the absolute
+value of each coefficient, while \link ArrayBase::sqrt() .sqrt() \endlink computes the square root of the
+coefficients. If you have two arrays of the same size, you can call \link ArrayBase::min() .min() \endlink to
+construct the array whose coefficients are the minimum of the corresponding coefficients of the two given
+arrays. These operations are illustrated in the following example.
+
+<table class="tutorial_code"><tr><td>
+Example: \include Tutorial_ArrayClass_cwise_other.cpp
+</td>
+<td>
+Output: \verbinclude Tutorial_ArrayClass_cwise_other.out
+</td></tr></table>
+
+More coefficient-wise operations can be found in the \ref QuickRefPage.
+

 \section TutorialArrayClassConvert Converting between array and matrix expressions

-It is possible to wrap a matrix expression as an array expression and conversely. This gives access to all operations
-regardless of the choice of declaring objects as arrays or as matrices.
+When should you use objects of the Matrix class and when should you use objects of the Array class? You cannot
+apply Matrix operations on arrays, or Array operations on matrices. Thus, if you need to do linear algebraic
+operations such as matrix multiplication, then you should use matrices; if you need to do coefficient-wise
+operations, then you should use arrays. However, sometimes it is not that simple, but you need to use both
+Matrix and Array operations. In that case, you need to convert a matrix to an array or reversely. This gives
+access to all operations regardless of the choice of declaring objects as arrays or as matrices.

-\link MatrixBase Matrix expressions \endlink have a \link MatrixBase::array() .array() \endlink method that
-'converts' them into \link ArrayBase array expressions \endlink, so that coefficient-wise operations
+\link MatrixBase Matrix expressions \endlink have an \link MatrixBase::array() .array() \endlink method that
+'converts' them into \link ArrayBase array expressions\endlink, so that coefficient-wise operations
 can be applied easily. Conversely, \link ArrayBase array expressions \endlink
 have a \link ArrayBase::matrix() .matrix() \endlink method. As with all Eigen expression abstractions,
 this doesn't have any runtime cost (provided that you let your compiler optimize).
-
 Both \link MatrixBase::array() .array() \endlink and \link ArrayBase::matrix() .matrix() \endlink 
 can be used as rvalues and as lvalues.

-Mixing matrices and arrays in an expression is forbidden with Eigen. However,
+Mixing matrices and arrays in an expression is forbidden with Eigen. For instance, you cannot add a amtrix and
+array directly; the operands of a \c + operator should either both be matrices or both be arrays. However,
 it is easy to convert from one to the other with \link MatrixBase::array() .array() \endlink and 
-\link ArrayBase::matrix() .matrix() \endlink.
+\link ArrayBase::matrix() .matrix()\endlink. The exception to this rule is the assignment operator: it is
+allowed to assign a matrix expression to an array variable, or to assign an array expression to a matrix
+variable.

-On the other hand, assigning a matrix (resp. array) expression to an array (resp. matrix) expression is allowed.
+The following example shows how to use array operations on a Matrix object by employing the 
+\link MatrixBase::array() .array() \endlink method. For example, the statement 
+<tt>result = m.array() * n.array()</tt> takes two matrices \c m and \c n, converts them both to an array, uses
+* to multiply them coefficient-wise and assigns the result to the matrix variable \c result (this is legal
+because Eigen allows assigning array expressions to matrix variables). 

-\subsection TutorialArrayClassInteropMatrix Matrix to array example
-The following example shows how to use array operations on a Matrix object by employing
-the \link MatrixBase::array() .array() \endlink method:
+As a matter of fact, this usage case is so common that Eigen provides a \link MatrixBase::cwiseProduct()
+.cwiseProduct() \endlink method for matrices to compute the coefficient-wise product. This is also shown in
+the example program.

 <table class="tutorial_code"><tr><td>
 \include Tutorial_ArrayClass_interop_matrix.cpp
@ -143,21 +178,13 @@ Output:
 \verbinclude Tutorial_ArrayClass_interop_matrix.out
 </td></tr></table>

+Similarly, if \c array1 and \c array2 are arrays, then the expression <tt>array1.matrix() * array2.matrix()</tt>
+computes their matrix product.

-\subsection TutorialArrayClassInteropArray Array to matrix example
-The following example shows how to use matrix operations with an Array
-object by employing the \link ArrayBase::matrix() .matrix() \endlink method:
-
-<table class="tutorial_code"><tr><td>
-\include Tutorial_ArrayClass_interop_array.cpp
-</td>
-<td>
-Output:
-\verbinclude Tutorial_ArrayClass_interop_array.out
-</td></tr></table>
-
-\subsection TutorialArrayClassInteropCombination Example with combinations of .array() and .matrix()
-Here is a more advanced example:
+Here is a more advanced example. The expression <tt>(m.array() + 4).matrix() * m</tt> adds 4 to every
+coefficient in the matrix \c m and then computes the matrix product of the result with \c m. Similarly, the
+expression <tt>(m.array() * n.array()).matrix() * m</tt> computes the coefficient-wise product of the matrices
+\c m and \c n and then the matrix product of the result with \c m.

 <table class="tutorial_code"><tr><td>
 \include Tutorial_ArrayClass_interop.cpp
--- a/doc/C04_TutorialBlockOperations.dox
+++ b/doc/C04_TutorialBlockOperations.dox
@ -13,21 +13,22 @@ provided that you let your compiler optimize.

 \b Table \b of \b contents
  - \ref TutorialBlockOperationsUsing
-  - \ref TutorialBlockOperationsSyntax
-     - \ref TutorialBlockOperationsSyntaxColumnRows
-     - \ref TutorialBlockOperationsSyntaxCorners
+  - \ref TutorialBlockOperationsSyntaxColumnRows
+  - \ref TutorialBlockOperationsSyntaxCorners
+  - \ref TutorialBlockOperationsSyntaxVectors
+

 \section TutorialBlockOperationsUsing Using block operations

 The most general block operation in Eigen is called \link DenseBase::block() .block() \endlink.
-This function returns a block of size <tt>(p,q)</tt> whose origin is at <tt>(i,j)</tt> by using 
-the following syntax:
+This function returns a block of size <tt>(p,q)</tt> whose origin is at <tt>(i,j)</tt>.
+There are two versions, whose syntax is as follows:

 <table class="tutorial_code" align="center">
-<tr><td align="center">\b Block \b operation</td>
-<td align="center">Default \b version</td>
+<tr><td align="center">\b %Block \b operation</td>
+<td align="center">Default version</td>
 <td align="center">Optimized version when the<br>size is known at compile time</td></tr>
-<tr><td>Block of size <tt>(p,q)</tt>, starting at <tt>(i,j)</tt></td>
+<tr><td>%Block of size <tt>(p,q)</tt>, starting at <tt>(i,j)</tt></td>
    <td>\code
 matrix.block(i,j,p,q);\endcode </td>
    <td>\code 
@ -35,7 +36,15 @@ matrix.block<p,q>(i,j);\endcode </td>
 </tr>
 </table>

-Therefore, if we want to print the values of a block inside a matrix we can simply write:
+The default version is a method which takes four arguments. It can always be used. The optimized version
+takes two template arguments (the size of the block) and two normal arguments (the position of the block).
+It can only be used if the size of the block is known at compile time, but it may be faster than the
+non-optimized version, especially if the size of the block is small. Both versions can be used on fixed-size
+and dynamic-size matrices and arrays.
+
+The following program uses the default and optimized versions to print the values of several blocks inside a
+matrix.
+
 <table class="tutorial_code"><tr><td>
 \include Tutorial_BlockOperations_print_block.cpp
 </td>
@ -44,10 +53,15 @@ Output:
 \verbinclude Tutorial_BlockOperations_print_block.out
 </td></tr></table>

+In the above example the \link DenseBase::block() .block() \endlink function was employed 
+to read the values inside matrix \p m . However, blocks can also be used as lvalues, meaning that you can
+assign to a block. 

-In the previous example the \link DenseBase::block() .block() \endlink function was employed 
-to read the values inside matrix \p m . Blocks can also be used to perform operations and 
-assignments within matrices or arrays of different size:
+This is illustrated in the following example, which uses arrays instead of matrices.  The coefficients of the
+5-by-5 array \c n are first all set to 0.6, but then the 3-by-3 block in the middle is set to the values in 
+\c m . The penultimate line shows that blocks can be combined with matrices and arrays to create more complex
+expressions. Blocks of an array are an array expression, and thus the multiplication here is coefficient-wise
+multiplication.

 <table class="tutorial_code"><tr><td>
 \include Tutorial_BlockOperations_block_assignment.cpp
@ -57,55 +71,38 @@ Output:
 \verbinclude Tutorial_BlockOperations_block_assignment.out
 </td></tr></table>

+The \link DenseBase::block() .block() \endlink method is used for general block operations, but there are
+other methods for special cases. These are described in the rest of this page.

-Blocks can also be combined with matrices and arrays to create more complex expressions:

-\code
-  MatrixXf m(3,3), n(2,2);
-  MatrixXf p(3,3);
-  
-  m.block(0,0,2,2) = m.block(0,0,2,2) * n + p.block(1,1,2,2);
-\endcode
+\section TutorialBlockOperationsSyntaxColumnRows Columns and rows

-It is important to point out that \link DenseBase::block() .block() \endlink is the 
-general case for a block operation but there are many other useful block operations, 
-as described in the next section.
-
-\section TutorialBlockOperationsSyntax Block operation syntax
-The following tables show a summary of Eigen's block operations and how they are applied to 
-fixed- and dynamic-sized Eigen objects.
-
-\subsection TutorialBlockOperationsSyntaxColumnRows Columns and rows
-Other extremely useful block operations are \link DenseBase::col() .col() \endlink and 
-\link DenseBase::row() .row() \endlink which provide access to a 
-specific row or column. This is a special case in the sense that the syntax for fixed- and 
-dynamic-sized objects is exactly the same:
+Individual columns and rows are special cases of blocks. Eigen provides methods to easily access them:
+\link DenseBase::col() .col() \endlink and \link DenseBase::row() .row()\endlink. There is no syntax variant
+for an optimized version.

 <table class="tutorial_code" align="center">
-<tr><td align="center">\b Block \b operation</td>
+<tr><td align="center">\b %Block \b operation</td>
 <td align="center">Default version</td>
 <td align="center">Optimized version when the<br>size is known at compile time</td></tr>
 <tr><td>i<sup>th</sup> row
                    \link DenseBase::row() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.row(i);\endcode </td>
+matrix.row(i);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.row(i);\endcode </td>
+matrix.row(i);\endcode </td>
 </tr>
 <tr><td>j<sup>th</sup> column
                    \link DenseBase::col() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.col(j);\endcode </td>
+matrix.col(j);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.col(j);\endcode </td>
+matrix.col(j);\endcode </td>
 </tr>
 </table>

-A simple example demonstrating these feature follows:
+The argument for \p col() and \p row() is the index of the column or row to be accessed, starting at
+0. Therefore, \p col(0) will access the first column and \p col(1) the second one.

 <table class="tutorial_code"><tr><td>
 C++ code:
@ -113,94 +110,83 @@ C++ code:
 </td>
 <td>
 Output:
-\include Tutorial_BlockOperations_colrow.out
+\verbinclude Tutorial_BlockOperations_colrow.out
 </td></tr></table>


-\b NOTE: the argument for \p col() and \p row() is the index of the column or row to be accessed, 
-starting at 0. Therefore, \p col(0) will access the first column and \p col(1) the second one.
+\section TutorialBlockOperationsSyntaxCorners Corner-related operations

+Eigen also provides special methods for blocks that are flushed against one of the corners or sides of a
+matrix or array. For instance, \link DenseBase::topLeftCorner() .topLeftCorner() \endlink can be used to refer
+to a block in the top-left corner of a matrix. Use <tt>matrix.topLeftCorner(p,q)</tt> to access the block
+consisting of the coefficients <tt>matrix(i,j)</tt> with \c i &lt; \c p and \c j &lt; \c q. As an other
+example, blocks consisting of whole rows flushed against the top side of the matrix can be accessed by
+\link DenseBase::topRows() .topRows() \endlink. 
+
+The different possibilities are summarized in the following table:

-\subsection TutorialBlockOperationsSyntaxCorners Corner-related operations
 <table class="tutorial_code" align="center">
-<tr><td align="center">\b Block \b operation</td>
+<tr><td align="center">\b %Block \b operation</td>
 <td align="center">Default version</td>
 <td align="center">Optimized version when the<br>size is known at compile time</td></tr>
 <tr><td>Top-left p by q block \link DenseBase::topLeftCorner() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.topLeftCorner(p,q);\endcode </td>
+matrix.topLeftCorner(p,q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.topLeftCorner<p,q>();\endcode </td>
+matrix.topLeftCorner<p,q>();\endcode </td>
 </tr>
 <tr><td>Bottom-left p by q block
              \link DenseBase::bottomLeftCorner() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.bottomLeftCorner(p,q);\endcode </td>
+matrix.bottomLeftCorner(p,q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.bottomLeftCorner<p,q>();\endcode </td>
+matrix.bottomLeftCorner<p,q>();\endcode </td>
 </tr>
 <tr><td>Top-right p by q block
              \link DenseBase::topRightCorner() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.topRightCorner(p,q);\endcode </td>
+matrix.topRightCorner(p,q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.topRightCorner<p,q>();\endcode </td>
+matrix.topRightCorner<p,q>();\endcode </td>
 </tr>
 <tr><td>Bottom-right p by q block
               \link DenseBase::bottomRightCorner() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.bottomRightCorner(p,q);\endcode </td>
+matrix.bottomRightCorner(p,q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.bottomRightCorner<p,q>();\endcode </td>
+matrix.bottomRightCorner<p,q>();\endcode </td>
 </tr>
-<tr><td>Block containing the first q rows
+<tr><td>%Block containing the first q rows
                   \link DenseBase::topRows() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.topRows(q);\endcode </td>
+matrix.topRows(q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.topRows<q>();\endcode </td>
+matrix.topRows<q>();\endcode </td>
 </tr>
-<tr><td>Block containing the last q rows
+<tr><td>%Block containing the last q rows
                    \link DenseBase::bottomRows() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.bottomRows(q);\endcode </td>
+matrix.bottomRows(q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.bottomRows<q>();\endcode </td>
+matrix.bottomRows<q>();\endcode </td>
 </tr>
-<tr><td>Block containing the first p columns
+<tr><td>%Block containing the first p columns
                    \link DenseBase::leftCols() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.leftCols(p);\endcode </td>
+matrix.leftCols(p);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.leftCols<p>();\endcode </td>
+matrix.leftCols<p>();\endcode </td>
 </tr>
-<tr><td>Block containing the last q columns
+<tr><td>%Block containing the last q columns
                    \link DenseBase::rightCols() * \endlink</td>
    <td>\code
-MatrixXf m;
-std::cout << m.rightCols(q);\endcode </td>
+matrix.rightCols(q);\endcode </td>
    <td>\code 
-Matrix3f m;
-std::cout << m.rightCols<q>();\endcode </td>
+matrix.rightCols<q>();\endcode </td>
 </tr>
 </table>

-
-Here is a simple example showing the power of the operations presented above:
+Here is a simple example illustrating the use of the operations presented above:

 <table class="tutorial_code"><tr><td>
 C++ code:
@ -208,49 +194,38 @@ C++ code:
 </td>
 <td>
 Output:
-\include Tutorial_BlockOperations_corner.out
+\verbinclude Tutorial_BlockOperations_corner.out
 </td></tr></table>


+\section TutorialBlockOperationsSyntaxVectors Block operations for vectors

-
-
-
-
-
-\subsection TutorialBlockOperationsSyntaxVectors Block operations for vectors
-Eigen also provides a set of block operations designed specifically for vectors:
+Eigen also provides a set of block operations designed specifically for vectors and one-dimensional arrays:

 <table class="tutorial_code" align="center">
-<tr><td align="center">\b Block \b operation</td>
+<tr><td align="center">\b %Block \b operation</td>
 <td align="center">Default version</td>
 <td align="center">Optimized version when the<br>size is known at compile time</td></tr>
-<tr><td>Block containing the first \p n elements 
+<tr><td>%Block containing the first \p n elements 
                    \link DenseBase::head() * \endlink</td>
    <td>\code
-VectorXf v;
-std::cout << v.head(n);\endcode </td>
+vector.head(n);\endcode </td>
    <td>\code 
-Vector3f v;
-std::cout << v.head<n>();\endcode </td>
+vector.head<n>();\endcode </td>
 </tr>
-<tr><td>Block containing the last \p n elements
+<tr><td>%Block containing the last \p n elements
                    \link DenseBase::tail() * \endlink</td>
    <td>\code
-VectorXf v;
-std::cout << v.tail(n);\endcode </td>
+vector.tail(n);\endcode </td>
    <td>\code 
-Vector3f m;
-std::cout << v.tail<n>();\endcode </td>
+vector.tail<n>();\endcode </td>
 </tr>
-<tr><td>Block containing \p n elements, starting at position \p i
+<tr><td>%Block containing \p n elements, starting at position \p i
                    \link DenseBase::segment() * \endlink</td>
    <td>\code
-VectorXf v;
-std::cout << v.segment(i,n);\endcode </td>
+vector.segment(i,n);\endcode </td>
    <td>\code 
-Vector3f m;
-std::cout << v.segment<n>(i);\endcode </td>
+vector.segment<n>(i);\endcode </td>
 </tr>
 </table>

@ -262,7 +237,7 @@ C++ code:
 </td>
 <td>
 Output:
-\include Tutorial_BlockOperations_vector.out
+\verbinclude Tutorial_BlockOperations_vector.out
 </td></tr></table>

 \li \b Next: \ref TutorialAdvancedInitialization
--- a/doc/C07_TutorialReductionsVisitorsBroadcasting.dox
+++ b/doc/C07_TutorialReductionsVisitorsBroadcasting.dox
@ -30,7 +30,7 @@ which returns the addition of all the coefficients inside a given matrix or arra
 Example: \include tut_arithmetic_redux_basic.cpp
 </td>
 <td>
-Output: \include tut_arithmetic_redux_basic.out
+Output: \verbinclude tut_arithmetic_redux_basic.out
 </td></tr></table>

 The \em trace of a matrix, as returned by the function \c trace(), is the sum of the diagonal coefficients and can also be computed as efficiently using <tt>a.diagonal().sum()</tt>, as we will see later on.
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@ -1217,12 +1217,14 @@ PREDEFINED             = EIGEN_EMPTY_STRUCT \
 # The macro definition that is found in the sources will be used.
 # Use the PREDEFINED tag if you want to use a different macro definition.

-EXPAND_AS_DEFINED      = EIGEN_MAKE_SCALAR_OPS \
-                         EIGEN_MAKE_TYPEDEFS \
+EXPAND_AS_DEFINED      = EIGEN_MAKE_TYPEDEFS \
                         EIGEN_MAKE_FIXED_TYPEDEFS \
                         EIGEN_MAKE_TYPEDEFS_ALL_SIZES \
+                         EIGEN_MAKE_CWISE_BINARY_OP \
                         EIGEN_CWISE_UNOP_RETURN_TYPE \
                         EIGEN_CWISE_BINOP_RETURN_TYPE \
+                         EIGEN_CWISE_PRODUCT_RETURN_TYPE \
+                         EIGEN_CURRENT_STORAGE_BASE_CLASS \
                         _EIGEN_GENERIC_PUBLIC_INTERFACE

 # If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then
--- a/doc/QuickReference.dox
+++ b/doc/QuickReference.dox
@ -218,11 +218,13 @@ x = FixedXD::Zero();
 x = FixedXD::Ones();
 x = FixedXD::Constant(value);
 x = FixedXD::Random();
+x = FixedXD::LinSpaced(low, high, size);

 x.setZero();
 x.setOnes();
 x.setConstant(value);
 x.setRandom();
+x.setLinSpaced(low, high, size);
 \endcode
  </td>
  <td>
@ -234,11 +236,13 @@ x = Dynamic2D::Zero(rows, cols);
 x = Dynamic2D::Ones(rows, cols);
 x = Dynamic2D::Constant(rows, cols, value);
 x = Dynamic2D::Random(rows, cols);
+N/A

 x.setZero(rows, cols);
 x.setOnes(rows, cols);
 x.setConstant(rows, cols, value);
 x.setRandom(rows, cols);
+N/A
 \endcode
  </td>
  <td>
@ -250,11 +254,13 @@ x = Dynamic1D::Zero(size);
 x = Dynamic1D::Ones(size);
 x = Dynamic1D::Constant(size, value);
 x = Dynamic1D::Random(size);
+x = Dynamic1D::LinSpaced(low, high, size);

 x.setZero(size);
 x.setOnes(size);
 x.setConstant(size, value);
 x.setRandom(size);
+x.setLinSpaced(low, high, size);
 \endcode
  </td>
 </tr>
@ -595,7 +601,7 @@ m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<Eigen::Lower>() \endcode
 Solving linear equations:\n
 \f$ M_2 := L_1^{-1} M_2 \f$ \n
 \f$ M_3 := {L_1^*}^{-1} M_3 \f$ \n
-\f$ M_4 := M_3 U_1^{-1} \f$
+\f$ M_4 := M_4 U_1^{-1} \f$
 </td><td>\n \code
 L1.triangularView<Eigen::UnitLower>().solveInPlace(M2)
 L1.triangularView<Eigen::Lower>().adjoint().solveInPlace(M3)
--- a/doc/examples/DenseBase_middleCols_int.cpp
+++ b/doc/examples/DenseBase_middleCols_int.cpp
@ -0,0 +1,15 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main(void)
+{
+    int const N = 5;
+    MatrixXi A(N,N);
+    A.setRandom();
+    cout << "A =\n" << A << '\n' << endl;
+    cout << "A(1..3,:) =\n" << A.middleCols(1,3) << endl;
+    return 0;
+}
--- a/doc/examples/DenseBase_middleRows_int.cpp
+++ b/doc/examples/DenseBase_middleRows_int.cpp
@ -0,0 +1,15 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main(void)
+{
+    int const N = 5;
+    MatrixXi A(N,N);
+    A.setRandom();
+    cout << "A =\n" << A << '\n' << endl;
+    cout << "A(2..3,:) =\n" << A.middleRows(2,2) << endl;
+    return 0;
+}
--- a/doc/examples/DenseBase_template_int_middleCols.cpp
+++ b/doc/examples/DenseBase_template_int_middleCols.cpp
@ -0,0 +1,15 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main(void)
+{
+    int const N = 5;
+    MatrixXi A(N,N);
+    A.setRandom();
+    cout << "A =\n" << A << '\n' << endl;
+    cout << "A(:,1..3) =\n" << A.middleCols<3>(1) << endl;
+    return 0;
+}
--- a/doc/examples/DenseBase_template_int_middleRows.cpp
+++ b/doc/examples/DenseBase_template_int_middleRows.cpp
@ -0,0 +1,15 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main(void)
+{
+    int const N = 5;
+    MatrixXi A(N,N);
+    A.setRandom();
+    cout << "A =\n" << A << '\n' << endl;
+    cout << "A(1..3,:) =\n" << A.middleRows<3>(1) << endl;
+    return 0;
+}
--- a/doc/examples/Tutorial_ArrayClass_accessors.cpp
+++ b/doc/examples/Tutorial_ArrayClass_accessors.cpp
@ -8,17 +8,17 @@ int main()
 {
  ArrayXXf  m(2,2);
  
-  //assign some values coefficient by coefficient
+  // assign some values coefficient by coefficient
  m(0,0) = 1.0; m(0,1) = 2.0;
-  m(1,0) = 3.0; m(1,1) = 4.0;
+  m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);
  
-  //print values to standard output
+  // print values to standard output
  cout << m << endl << endl;
 
  // using the comma-initializer is also allowed
  m << 1.0,2.0,
       3.0,4.0;
     
-  //print values to standard output
+  // print values to standard output
  cout << m << endl;
 }
--- a/doc/examples/Tutorial_ArrayClass_addition.cpp
+++ b/doc/examples/Tutorial_ArrayClass_addition.cpp
@ -8,15 +8,16 @@ int main()
 {
  ArrayXXf a(3,3);
  ArrayXXf b(3,3);
-  
  a << 1,2,3,
       4,5,6,
       7,8,9;
-
  b << 1,2,3,
       1,2,3,
       1,2,3;
       
-  cout << "a + b = " << endl
-       << a + b << endl;
+  // Adding two arrays
+  cout << "a + b = " << endl << a + b << endl << endl;
+
+  // Subtracting a scalar from an array
+  cout << "a - 2 = " << endl << a - 2 << endl;
 }
--- a/doc/examples/Tutorial_ArrayClass_addition_scalar.cpp
+++ b/doc/examples/Tutorial_ArrayClass_addition_scalar.cpp
@ -1,17 +0,0 @@
-#include <Eigen/Dense>
-#include <iostream>
-
-using namespace Eigen;
-using namespace std;
-
-int main()
-{
-  ArrayXXf a(3,3);
-  
-  a << 1,2,3,
-       4,5,6,
-       7,8,9;
-
-  cout << "a + 2 = " << endl
-       << a + 2 << endl;
-}
--- a/doc/examples/Tutorial_ArrayClass_cwise_other.cpp
+++ b/doc/examples/Tutorial_ArrayClass_cwise_other.cpp
@ -0,0 +1,19 @@
+#include <Eigen/Dense>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  ArrayXf a = ArrayXf::Random(5);
+  a *= 2;
+  cout << "a =" << endl 
+       << a << endl;
+  cout << "a.abs() =" << endl 
+       << a.abs() << endl;
+  cout << "a.abs().sqrt() =" << endl 
+       << a.abs().sqrt() << endl;
+  cout << "a.min(a.abs().sqrt()) =" << endl 
+       << a.min(a.abs().sqrt()) << endl;
+}
--- a/doc/examples/Tutorial_ArrayClass_interop.cpp
+++ b/doc/examples/Tutorial_ArrayClass_interop.cpp
@ -8,31 +8,15 @@ int main()
 {
  MatrixXf m(2,2);
  MatrixXf n(2,2);
-  
  MatrixXf result(2,2);

-  //initialize matrices
  m << 1,2,
       3,4;
-
  n << 5,6,
       7,8;
  
-  // mix of array and matrix operations
-  //   first coefficient-wise addition
-  //   then the result is used with matrix multiplication
  result = (m.array() + 4).matrix() * m;
-
-  cout << "-- Combination 1: --" << endl
-    << result << endl << endl;
-
-
-  // mix of array and matrix operations
-  //   first coefficient-wise multiplication
-  //   then the result is used with matrix multiplication
+  cout << "-- Combination 1: --" << endl << result << endl << endl;
  result = (m.array() * n.array()).matrix() * m;
-
-  cout << "-- Combination 2: --" << endl
-    << result << endl << endl;
-
+  cout << "-- Combination 2: --" << endl << result << endl << endl;
 }
--- a/doc/examples/Tutorial_ArrayClass_interop_array.cpp
+++ b/doc/examples/Tutorial_ArrayClass_interop_array.cpp
@ -1,34 +0,0 @@
-#include <Eigen/Dense>
-#include <iostream>
-
-using namespace Eigen;
-using namespace std;
-
-int main()
-{
-  ArrayXXf m(2,2);
-  ArrayXXf n(2,2);
-  
-  ArrayXXf result(2,2);
-
-  //initialize arrays
-  m << 1,2,
-       3,4;
-
-  n << 5,6,
-       7,8;
-
-  
-  // --> array multiplication
-  result = m * n;
-
-  cout << "-- Array m*n: --" << endl
-    << result << endl << endl;
-
-
-  // --> Matrix multiplication
-  result = m.matrix() * n.matrix();
-  
-  cout << "-- Matrix m*n: --" << endl
-    << result << endl << endl;
-}
--- a/doc/examples/Tutorial_ArrayClass_interop_matrix.cpp
+++ b/doc/examples/Tutorial_ArrayClass_interop_matrix.cpp
@ -8,34 +8,19 @@ int main()
 {
  MatrixXf m(2,2);
  MatrixXf n(2,2);
-  
  MatrixXf result(2,2);

-  //initialize matrices
  m << 1,2,
       3,4;
-
  n << 5,6,
       7,8;

-  
-  // --> matrix multiplication
  result = m * n;
-
-  cout << "-- Matrix m*n: --" << endl
-    << result << endl << endl;
-
-
-  // --> coeff-wise multiplication
+  cout << "-- Matrix m*n: --" << endl << result << endl << endl;
  result = m.array() * n.array();
-  
-  cout << "-- Array m*n: --" << endl
-    << result << endl << endl;
-  
-  
-  // ->> coeff-wise addition of a scalar
+  cout << "-- Array m*n: --" << endl << result << endl << endl;
+  result = m.cwiseProduct(n);
+  cout << "-- With cwiseProduct: --" << endl << result << endl << endl;
  result = m.array() + 4;
-  
-  cout << "-- Array m + 4: --" << endl
-    << result << endl << endl;
+  cout << "-- Array m + 4: --" << endl << result << endl << endl;
 }
--- a/doc/examples/Tutorial_ArrayClass_mult.cpp
+++ b/doc/examples/Tutorial_ArrayClass_mult.cpp
@ -8,13 +8,9 @@ int main()
 {
  ArrayXXf a(2,2);
  ArrayXXf b(2,2);
-  
  a << 1,2,
       3,4;
-
  b << 5,6,
       7,8;
-
-  cout << "a * b = " << endl
-       << a * b << endl;
+  cout << "a * b = " << endl << a * b << endl;
 }
--- a/doc/examples/Tutorial_BlockOperations_block_assignment.cpp
+++ b/doc/examples/Tutorial_BlockOperations_block_assignment.cpp
@ -6,26 +6,13 @@ using namespace Eigen;

 int main()
 {
-  MatrixXf m(3,3), n(2,2);
-  
+  Array33f m;
  m << 1,2,3,
       4,5,6,
       7,8,9;
-       
-  // assignment through a block operation,
-  //  block as rvalue
-  n = m.block(0,0,2,2);
-  
-  //print n
+  Array<float,5,5> n = Array<float,5,5>::Constant(0.6);
+  n.block(1,1,3,3) = m;
  cout << "n = " << endl << n << endl << endl;
-  
-  
-  n << 1,1,
-       1,1;
-        
-  // block as lvalue
-  m.block(0,0,2,2) = n;
-  
-  //print m
-  cout << "m = " << endl << m << endl;
+  Array33f res = n.block(0,0,3,3) * m;
+  cout << "res =" << endl << res << endl;
 }
--- a/doc/examples/Tutorial_BlockOperations_colrow.cpp
+++ b/doc/examples/Tutorial_BlockOperations_colrow.cpp
@ -1,15 +1,14 @@
 #include <Eigen/Dense>
 #include <iostream>
-using namespace Eigen;

 int main()
 {
-  MatrixXf m(3,3);
-  
+  Eigen::MatrixXf m(3,3);
  m << 1,2,3,
       4,5,6,
       7,8,9;
-     
-  std::cout << "2nd Row: "
-    << m.row(1) << std::endl;
+  std::cout << "2nd Row: " << m.row(1) << std::endl;
+  m.col(0) += m.col(2);
+  std::cout << "m after adding third column to first:\n";
+  std::cout << m << std::endl;
 }
--- a/doc/examples/Tutorial_BlockOperations_corner.cpp
+++ b/doc/examples/Tutorial_BlockOperations_corner.cpp
@ -2,26 +2,16 @@
 #include <iostream>

 using namespace std;
-using namespace Eigen;

 int main()
 {
-  MatrixXf m(4,4);
-  
+  Eigen::Matrix4f m;
  m << 1, 2, 3, 4,
       5, 6, 7, 8,
       9, 10,11,12,
       13,14,15,16;
-
-  //print first two columns
-  cout << "-- leftCols(2) --" << endl
-    << m.leftCols(2) << endl << endl;
-  
-  //print last two rows
-  cout << "-- bottomRows(2) --" << endl
-    << m.bottomRows(2) << endl << endl;
-    
-  //print top-left 2x3 corner
-  cout << "-- topLeftCorner(2,3) --" << endl
-    << m.topLeftCorner(2,3) << endl;
+  cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl;
+  cout << "m.bottomRows<2>() =" << endl << m.bottomRows<2>() << endl << endl;
+  m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();
+  cout << "After assignment, m = " << endl << m << endl;
 }
--- a/doc/examples/Tutorial_BlockOperations_print_block.cpp
+++ b/doc/examples/Tutorial_BlockOperations_print_block.cpp
@ -1,14 +1,18 @@
 #include <Eigen/Dense>
 #include <iostream>
-using namespace Eigen;

 int main()
 {
-  MatrixXf m(3,3);
-  
-  m << 1,2,3,
-       4,5,6,
-       7,8,9;
-     
-  std::cout << m.block(0,0,2,2) << std::endl;
+  Eigen::MatrixXf m(4,4);
+  m <<  1, 2, 3, 4,
+        5, 6, 7, 8,
+        9,10,11,12,
+       13,14,15,16;
+  std::cout << "Block in the middle" << std::endl;
+  std::cout << m.block<2,2>(1,1) << std::endl << std::endl;
+  for (int i = 1; i < 4; ++i) 
+  {
+    std::cout << "Block of size " << i << std::endl;
+    std::cout << m.block(0,0,i,i) << std::endl << std::endl;
+  }
 }
--- a/doc/examples/Tutorial_BlockOperations_vector.cpp
+++ b/doc/examples/Tutorial_BlockOperations_vector.cpp
@ -2,23 +2,13 @@
 #include <iostream>

 using namespace std;
-using namespace Eigen;

 int main()
 {
-  VectorXf v(6);
-  
+  Eigen::ArrayXf v(6);
  v << 1, 2, 3, 4, 5, 6;
-
-  //print first three elements
-  cout << "-- head(3) --" << endl
-    << v.head(3) << endl << endl;
-  
-  //print last three elements
-  cout << "-- tail(3) --" << endl
-    << v.tail(3) << endl << endl;
-    
-  //print between 2nd and 5th elem. inclusive
-  cout << "-- segment(1,4) --" << endl
-    << v.segment(1,4) << endl;
+  cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
+  cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
+  v.segment(1,4) *= 2;
+  cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
 }
--- a/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_maxnorm.cpp
+++ b/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_maxnorm.cpp
@ -2,13 +2,14 @@
 #include <Eigen/Dense>

 using namespace std;
+using namespace Eigen;
 int main()
 {
-  Eigen::MatrixXf mat(2,4);
+  MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
-  int   maxIndex;
+  MatrixXf::Index   maxIndex;
  float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);
  
  std::cout << "Maximum sum at position " << maxIndex << std::endl;
--- a/doc/examples/tut_matrix_coefficient_accessors.cpp
+++ b/doc/examples/tut_matrix_coefficient_accessors.cpp
@ -15,5 +15,4 @@ int main()
  v(0) = 4;
  v(1) = v(0) - 1;
  std::cout << "Here is the vector v:\n" << v << std::endl;
-
 }
--- a/doc/snippets/DenseBase_LinSpaced.cpp
+++ b/doc/snippets/DenseBase_LinSpaced.cpp
@ -1,2 +1,2 @@
-cout << VectorXi::LinSpaced(7,10,4).transpose() << endl;
-cout << VectorXd::LinSpaced(0.0,1.0,5).transpose() << endl;
+cout << VectorXi::LinSpaced(4,7,10).transpose() << endl;
+cout << VectorXd::LinSpaced(5,0.0,1.0).transpose() << endl;
--- a/doc/snippets/DenseBase_LinSpaced_seq.cpp
+++ b/doc/snippets/DenseBase_LinSpaced_seq.cpp
@ -1,2 +1,2 @@
-cout << VectorXi::LinSpaced(Sequential,7,10,4).transpose() << endl;
-cout << VectorXd::LinSpaced(Sequential,0.0,1.0,5).transpose() << endl;
+cout << VectorXi::LinSpaced(Sequential,4,7,10).transpose() << endl;
+cout << VectorXd::LinSpaced(Sequential,5,0.0,1.0).transpose() << endl;
--- a/doc/snippets/DenseBase_setLinSpaced.cpp
+++ b/doc/snippets/DenseBase_setLinSpaced.cpp
@ -1,3 +1,3 @@
 VectorXf v;
-v.setLinSpaced(0.5f,1.5f,5).transpose();
+v.setLinSpaced(5,0.5f,1.5f).transpose();
 cout << v << endl;
--- a/doc/unsupported_modules.dox
+++ b/doc/unsupported_modules.dox
@ -61,4 +61,7 @@ namespace Eigen {
 /** \ingroup Unsupported_modules
  * \defgroup SparseExtra_Module */

+/** \ingroup Unsupported_modules
+  * \defgroup MpfrcxxSupport_Module */
+
 } // namespace Eigen
--- a/test/array_for_matrix.cpp
+++ b/test/array_for_matrix.cpp
@ -72,6 +72,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m)
  VERIFY_IS_APPROX(m3.rowwise() += rv1, m1.rowwise() + rv1);
  m3 = m1;
  VERIFY_IS_APPROX(m3.rowwise() -= rv1, m1.rowwise() - rv1);
+  
+  // empty objects
+  VERIFY_IS_APPROX(m1.block(0,0,0,cols).colwise().sum(),  RowVectorType::Zero(cols));
+  VERIFY_IS_APPROX(m1.block(0,0,rows,0).rowwise().prod(), ColVectorType::Ones(rows));
 }

 template<typename MatrixType> void comparisons(const MatrixType& m)
--- a/test/corners.cpp
+++ b/test/corners.cpp
@ -45,13 +45,20 @@ template<typename MatrixType> void corners(const MatrixType& m)
  COMPARE_CORNER(bottomLeftCorner(r,c), block(rows-r,0,r,c));
  COMPARE_CORNER(bottomRightCorner(r,c), block(rows-r,cols-c,r,c));

+  Index sr = ei_random<Index>(1,rows) - 1;
+  Index nr = ei_random<Index>(1,rows-sr);
+  Index sc = ei_random<Index>(1,cols) - 1;
+  Index nc = ei_random<Index>(1,cols-sc);
+
  COMPARE_CORNER(topRows(r), block(0,0,r,cols));
+  COMPARE_CORNER(middleRows(sr,nr), block(sr,0,nr,cols));
  COMPARE_CORNER(bottomRows(r), block(rows-r,0,r,cols));
  COMPARE_CORNER(leftCols(c), block(0,0,rows,c));
+  COMPARE_CORNER(middleCols(sc,nc), block(0,sc,rows,nc));
  COMPARE_CORNER(rightCols(c), block(0,cols-c,rows,c));
 }

-template<typename MatrixType, int CRows, int CCols> void corners_fixedsize()
+template<typename MatrixType, int CRows, int CCols, int SRows, int SCols> void corners_fixedsize()
 {
  MatrixType matrix = MatrixType::Random();
  const MatrixType const_matrix = MatrixType::Random();
@ -60,7 +67,9 @@ template<typename MatrixType, int CRows, int CCols> void corners_fixedsize()
    rows = MatrixType::RowsAtCompileTime,
    cols = MatrixType::ColsAtCompileTime,
    r = CRows,
-    c = CCols
+    c = CCols,
+	sr = SRows,
+	sc = SCols
  };

  VERIFY_IS_EQUAL((matrix.template topLeftCorner<r,c>()), (matrix.template block<r,c>(0,0)));
@ -69,8 +78,10 @@ template<typename MatrixType, int CRows, int CCols> void corners_fixedsize()
  VERIFY_IS_EQUAL((matrix.template bottomRightCorner<r,c>()), (matrix.template block<r,c>(rows-r,cols-c)));

  VERIFY_IS_EQUAL((matrix.template topRows<r>()), (matrix.template block<r,cols>(0,0)));
+  VERIFY_IS_EQUAL((matrix.template middleRows<r>(sr)), (matrix.template block<r,cols>(sr,0)));
  VERIFY_IS_EQUAL((matrix.template bottomRows<r>()), (matrix.template block<r,cols>(rows-r,0)));
  VERIFY_IS_EQUAL((matrix.template leftCols<c>()), (matrix.template block<rows,c>(0,0)));
+  VERIFY_IS_EQUAL((matrix.template middleCols<c>(sc)), (matrix.template block<rows,c>(0,sc)));
  VERIFY_IS_EQUAL((matrix.template rightCols<c>()), (matrix.template block<rows,c>(0,cols-c)));

  VERIFY_IS_EQUAL((const_matrix.template topLeftCorner<r,c>()), (const_matrix.template block<r,c>(0,0)));
@ -79,8 +90,10 @@ template<typename MatrixType, int CRows, int CCols> void corners_fixedsize()
  VERIFY_IS_EQUAL((const_matrix.template bottomRightCorner<r,c>()), (const_matrix.template block<r,c>(rows-r,cols-c)));

  VERIFY_IS_EQUAL((const_matrix.template topRows<r>()), (const_matrix.template block<r,cols>(0,0)));
+  VERIFY_IS_EQUAL((const_matrix.template middleRows<r>(sr)), (const_matrix.template block<r,cols>(sr,0)));
  VERIFY_IS_EQUAL((const_matrix.template bottomRows<r>()), (const_matrix.template block<r,cols>(rows-r,0)));
  VERIFY_IS_EQUAL((const_matrix.template leftCols<c>()), (const_matrix.template block<rows,c>(0,0)));
+  VERIFY_IS_EQUAL((const_matrix.template middleCols<c>(sc)), (const_matrix.template block<rows,c>(0,sc)));
  VERIFY_IS_EQUAL((const_matrix.template rightCols<c>()), (const_matrix.template block<rows,c>(0,cols-c)));
 }

@ -93,8 +106,8 @@ void test_corners()
    CALL_SUBTEST_4( corners(MatrixXcf(5, 7)) );
    CALL_SUBTEST_5( corners(MatrixXf(21, 20)) );

-    CALL_SUBTEST_1(( corners_fixedsize<Matrix<float, 1, 1>, 1, 1>() ));
-    CALL_SUBTEST_2(( corners_fixedsize<Matrix4d,2,2>() ));
-    CALL_SUBTEST_3(( corners_fixedsize<Matrix<int,10,12>,4,7>() ));
+    CALL_SUBTEST_1(( corners_fixedsize<Matrix<float, 1, 1>, 1, 1, 0, 0>() ));
+    CALL_SUBTEST_2(( corners_fixedsize<Matrix4d,2,2,1,1>() ));
+    CALL_SUBTEST_3(( corners_fixedsize<Matrix<int,10,12>,4,7,5,2>() ));
  }
 }
--- a/test/determinant.cpp
+++ b/test/determinant.cpp
@ -61,6 +61,9 @@ template<typename MatrixType> void determinant(const MatrixType& m)
  m2 = m1;
  m2.row(i) *= x;
  VERIFY_IS_APPROX(m2.determinant(), m1.determinant() * x);
+  
+  // check empty matrix
+  VERIFY_IS_APPROX(m2.block(0,0,0,0).determinant(), Scalar(1));
 }

 void test_determinant()
--- a/test/main.h
+++ b/test/main.h
@ -283,7 +283,7 @@ namespace Eigen

 namespace Eigen {

-template<typename T> inline typename NumTraits<T>::Real test_precision() { return T(0); }
+template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); }
 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>(); }
--- a/test/nullary.cpp
+++ b/test/nullary.cpp
@ -59,7 +59,7 @@ void testVectorType(const VectorType& base)

  // check whether the result yields what we expect it to do
  VectorType m(base);
-  m.setLinSpaced(low,high,size);
+  m.setLinSpaced(size,low,high);

  VectorType n(size);
  for (int i=0; i<size; ++i)
@ -68,7 +68,7 @@ void testVectorType(const VectorType& base)
  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );

  // random access version
-  m = VectorType::LinSpaced(low,high,size);
+  m = VectorType::LinSpaced(size,low,high);
  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );

  // These guys sometimes fail! This is not good. Any ideas how to fix them!?
@ -76,7 +76,7 @@ void testVectorType(const VectorType& base)
  //VERIFY( m(0) == low );

  // sequential access version
-  m = VectorType::LinSpaced(Sequential,low,high,size);
+  m = VectorType::LinSpaced(Sequential,size,low,high);
  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );

  // These guys sometimes fail! This is not good. Any ideas how to fix them!?
@ -86,12 +86,12 @@ void testVectorType(const VectorType& base)
  // check whether everything works with row and col major vectors
  Matrix<Scalar,Dynamic,1> row_vector(size);
  Matrix<Scalar,1,Dynamic> col_vector(size);
-  row_vector.setLinSpaced(low,high,size);
-  col_vector.setLinSpaced(low,high,size);
+  row_vector.setLinSpaced(size,low,high);
+  col_vector.setLinSpaced(size,low,high);
  VERIFY( row_vector.isApprox(col_vector.transpose(), NumTraits<Scalar>::epsilon()));

  Matrix<Scalar,Dynamic,1> size_changer(size+50);
-  size_changer.setLinSpaced(low,high,size);
+  size_changer.setLinSpaced(size,low,high);
  VERIFY( size_changer.size() == size );
 }

--- a/test/product_trmm.cpp
+++ b/test/product_trmm.cpp
@ -35,7 +35,7 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
  DenseIndex cols = ei_random<DenseIndex>(1,size);

  MatrixColMaj  triV(rows,cols), triH(cols,rows), upTri(cols,rows), loTri(rows,cols),
-                unitUpTri(cols,rows), unitLoTri(rows,cols);
+                unitUpTri(cols,rows), unitLoTri(rows,cols), strictlyUpTri(cols,rows), strictlyLoTri(rows,cols);
  MatrixColMaj  ge1(rows,cols), ge2(cols,rows), ge3;
  MatrixRowMaj  rge3;

@ -48,6 +48,8 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
  upTri = triH.template triangularView<Upper>();
  unitLoTri = triV.template triangularView<UnitLower>();
  unitUpTri = triH.template triangularView<UnitUpper>();
+  strictlyLoTri = triV.template triangularView<StrictlyLower>();
+  strictlyUpTri = triH.template triangularView<StrictlyUpper>();
  ge1.setRandom();
  ge2.setRandom();

@ -72,6 +74,11 @@ template<typename Scalar> void trmm(int size,int /*othersize*/)
  VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri);
  VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<UnitLower>(), ge2 * unitLoTri);
  VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<UnitUpper>() * ge2.adjoint(), ei_conj(s1) * unitLoTri.adjoint() * ge2.adjoint());
+
+  VERIFY_IS_APPROX( ge3 = triV.template triangularView<StrictlyLower>() * ge2, strictlyLoTri * ge2);
+  VERIFY_IS_APPROX( rge3.noalias() = ge2 * triV.template triangularView<StrictlyLower>(), ge2 * strictlyLoTri);
+  VERIFY_IS_APPROX( ge3 = ge2 * triV.template triangularView<StrictlyLower>(), ge2 * strictlyLoTri);
+  VERIFY_IS_APPROX( ge3 = (s1*triV).adjoint().template triangularView<StrictlyUpper>() * ge2.adjoint(), ei_conj(s1) * strictlyLoTri.adjoint() * ge2.adjoint());
 }

 void test_product_trmm()
--- a/test/redux.cpp
+++ b/test/redux.cpp
@ -64,6 +64,10 @@ template<typename MatrixType> void matrixRedux(const MatrixType& m)
  VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).prod(), m1.block(r0,c0,r1,c1).eval().prod());
  VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).real().minCoeff(), m1.block(r0,c0,r1,c1).real().eval().minCoeff());
  VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).real().maxCoeff(), m1.block(r0,c0,r1,c1).real().eval().maxCoeff());
+  
+  // test empty objects
+  VERIFY_IS_APPROX(m1.block(r0,c0,0,0).sum(),   Scalar(0));
+  VERIFY_IS_APPROX(m1.block(r0,c0,0,0).prod(),  Scalar(1));
 }

 template<typename VectorType> void vectorRedux(const VectorType& w)
@ -124,6 +128,13 @@ template<typename VectorType> void vectorRedux(const VectorType& w)
    VERIFY_IS_APPROX(minc, v.real().segment(i, size-2*i).minCoeff());
    VERIFY_IS_APPROX(maxc, v.real().segment(i, size-2*i).maxCoeff());
  }
+  
+  // test empty objects
+  VERIFY_IS_APPROX(v.head(0).sum(),   Scalar(0));
+  VERIFY_IS_APPROX(v.tail(0).prod(),  Scalar(1));
+  VERIFY_RAISES_ASSERT(v.head(0).mean());
+  VERIFY_RAISES_ASSERT(v.head(0).minCoeff());
+  VERIFY_RAISES_ASSERT(v.head(0).maxCoeff());
 }

 void test_redux()
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@ -250,6 +250,13 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
    VERIFY(countTrueNonZero==m2.nonZeros());
    VERIFY_IS_APPROX(m2, refM2);
  }
+  
+  // test sparseView
+  {
+    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
+    SparseMatrixType m2(rows, rows);
+    VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
+  }
 }

 void test_sparse_basic()
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@ -1,6 +1,6 @@
 set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials
                  CholmodSupport FFT NonLinearOptimization SparseExtra SuperLUSupport UmfPackSupport IterativeSolvers
-                  NumericalDiff Skyline TaucsSupport
+                  NumericalDiff Skyline TaucsSupport MPRealSupport
   )

 install(FILES
--- a/unsupported/Eigen/MPRealSupport
+++ b/unsupported/Eigen/MPRealSupport
@ -0,0 +1,160 @@
+// This file is part of a joint effort between Eigen, a lightweight C++ template library
+// for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/)
+//
+// Copyright (C) 2010 Pavel Holoborodko <pavel@holoborodko.com>
+// Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.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
+// this library. If not, see <http://www.gnu.org/licenses/>.
+// 
+// Contributors:
+// Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, Pere Constans
+
+#ifndef EIGEN_MPREALSUPPORT_MODULE_H
+#define EIGEN_MPREALSUPPORT_MODULE_H
+
+#include <mpreal.h>
+#include <Eigen/Core>
+
+namespace Eigen {
+  
+  /** \defgroup MPRealSupport_Module MPFRC++ Support module
+    *
+    * \code
+    * #include <Eigen/MPRealSupport>
+    * \endcode
+    *
+    * This module provides support for multi precision floating point numbers
+    * via the <a href="http://www.holoborodko.com/pavel/?page_id=12">MPFR C++</a>
+    * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>.
+    *
+    * Here is an example:
+    *
+\code
+#include <iostream>
+#include <Eigen/Mpfrc++Support>
+#include <Eigen/LU>
+using namespace mpfr;
+using namespace Eigen;
+int main()
+{
+  // set precision to 256 bits (double has only 53 bits)
+  mpreal::set_default_prec(256);
+  // Declare matrix and vector types with multi-precision scalar type
+  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
+  typedef Matrix<mpreal,Dynamic,1>        VectorXmp;
+
+  MatrixXmp A = MatrixXmp::Random(100,100);
+  VectorXmp b = VectorXmp::Random(100);
+
+  // Solve Ax=b using LU
+  VectorXmp x = A.lu().solve(b);
+  std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl;
+  return 0;
+}
+\endcode
+    *
+    */
+
+  template<> struct NumTraits<mpfr::mpreal>
+    : GenericNumTraits<mpfr::mpreal>
+  {
+    enum {
+      IsInteger = 0,
+      IsSigned = 1,
+      IsComplex = 0,
+      ReadCost = 10,
+      AddCost = 10,
+      MulCost = 40
+    };
+
+    typedef mpfr::mpreal Real;
+    typedef mpfr::mpreal NonInteger;
+
+    inline static mpfr::mpreal highest() { return  mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
+    inline static mpfr::mpreal lowest()  { return -mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
+
+    inline static Real epsilon()
+    {
+      return mpfr::machine_epsilon(mpfr::mpreal::get_default_prec());
+    }
+    inline static Real dummy_precision()
+    {
+      unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1)*90)/100;
+      return mpfr::machine_epsilon(weak_prec);
+    }
+  };
+
+  template<> mpfr::mpreal ei_random<mpfr::mpreal>()
+  {
+#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+    static gmp_randstate_t state;
+    static bool isFirstTime = true;
+
+    if(isFirstTime)
+    {
+      gmp_randinit_default(state);
+      gmp_randseed_ui(state,(unsigned)time(NULL));
+      isFirstTime = false;
+    }
+
+    return mpfr::urandom(state)*2-1;
+#else
+    return mpfr::mpreal(ei_random<double>());
+#endif
+  }
+
+  template<> mpfr::mpreal ei_random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
+  {
+    return a + (b-a) * ei_random<mpfr::mpreal>();
+  }
+}
+
+namespace mpfr {
+
+  inline const mpreal& ei_conj(const mpreal& x)  { return x; }
+  inline const mpreal& ei_real(const mpreal& x)  { return x; }
+  inline mpreal ei_imag(const mpreal&)    { return 0.0; }
+  inline mpreal ei_abs(const mpreal&  x)  { return fabs(x); }
+  inline mpreal ei_abs2(const mpreal& x)  { return x*x; }
+  inline mpreal ei_sqrt(const mpreal& x)  { return sqrt(x); }
+  inline mpreal ei_exp(const mpreal&  x)  { return exp(x); }
+  inline mpreal ei_log(const mpreal&  x)  { return log(x); }
+  inline mpreal ei_sin(const mpreal&  x)  { return sin(x); }
+  inline mpreal ei_cos(const mpreal&  x)  { return cos(x); }
+  inline mpreal ei_pow(const mpreal& x, mpreal& y)  { return pow(x, y); }
+
+  bool ei_isMuchSmallerThan(const mpreal& a, const mpreal& b, const mpreal& prec)
+  {
+    return ei_abs(a) <= abs(b) * prec;
+  }
+
+  inline bool ei_isApprox(const mpreal& a, const mpreal& b, const mpreal& prec)
+  {
+  return ei_abs(a - b) <= min(abs(a), abs(b)) * prec;
+  }
+
+  inline bool ei_isApproxOrLessThan(const mpreal& a, const mpreal& b, const mpreal& prec)
+  {
+    return a <= b || ei_isApprox(a, b, prec);
+  }
+}
+
+#endif // EIGEN_MPREALSUPPORT_MODULE_H
--- a/unsupported/Eigen/OpenGLSupport
+++ b/unsupported/Eigen/OpenGLSupport
@ -0,0 +1,190 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.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_OPENGL_MODULE
+#define EIGEN_OPENGL_MODULE
+
+#include <Eigen/Geometry>
+#include <GL/gl.h>
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+  * \defgroup OpenGLSUpport_Module OpenGL Support module
+  *
+  * This module provides wrapper functions for a couple of OpenGL functions
+  * which simplify the way to pass Eigen's object to openGL.
+  * Here is an exmaple:
+  * 
+  * \code
+  * // You need to add path_to_eigen/unsupported to your include path.
+  * #include <Eigen/OpenGLSupport>
+  * // ...
+  * Vector3f x, y;
+  * Matrix3f rot;
+  * 
+  * glVertex(y + x * rot);
+  * 
+  * Quaternion q;
+  * glRotate(q);
+  * 
+  * // ...
+  * \endcode
+  *
+  */
+//@{
+
+#define EIGEN_GL_FUNC_DECLARATION(FUNC)                                                                             \
+template< typename XprType,                                                                                         \
+          typename Scalar = typename XprType::Scalar,                                                               \
+          int Rows = XprType::RowsAtCompileTime,                                                                    \
+          int Cols = XprType::ColsAtCompileTime,                                                                    \
+          bool IsGLCompatible = bool(XprType::Flags&LinearAccessBit)                                                \
+                             && bool(XprType::Flags&DirectAccessBit)                                                \
+                             && (XprType::IsVectorAtCompileTime || (XprType::Flags&RowMajorBit)==0)>                \
+struct EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl);                                                                     \
+                                                                                                                    \
+template<typename XprType, typename Scalar, int Rows, int Cols>                                                     \
+struct EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<XprType,Scalar,Rows,Cols,false> {                                    \
+  inline static void run(const XprType& p) {                                                                        \
+    EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<typename ei_plain_matrix_type_column_major<XprType>::type>::run(p); }   \
+};                                                                                                                  \
+                                                                                                                    \
+template<typename Derived> inline void FUNC(const Eigen::DenseBase<Derived>& p) {                                   \
+  EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<Derived>::run(p.derived());                                               \
+}
+
+
+#define EIGEN_GL_FUNC_SPECIALIZATION_MAT(FUNC,SCALAR,ROWS,COLS,SUFFIX)                                              \
+  template< typename XprType> struct EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<XprType, SCALAR, ROWS, COLS, true> {   \
+    inline static void run(const XprType& p) { FUNC##SUFFIX(p.data()); }                                            \
+  };
+
+  
+#define EIGEN_GL_FUNC_SPECIALIZATION_VEC(FUNC,SCALAR,SIZE,SUFFIX)                                                   \
+  template< typename XprType> struct EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<XprType, SCALAR, SIZE, 1, true> {      \
+    inline static void run(const XprType& p) { FUNC##SUFFIX(p.data()); }                                            \
+  };                                                                                                                \
+  template< typename XprType> struct EIGEN_CAT(EIGEN_CAT(ei_gl_,FUNC),_impl)<XprType, SCALAR, 1, SIZE, true> {      \
+    inline static void run(const XprType& p) { FUNC##SUFFIX(p.data()); }                                            \
+  };
+  
+EIGEN_GL_FUNC_DECLARATION       (glVertex)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,int,    2,2iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,short,  2,2sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,float,  2,2fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,double, 2,2dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,int,    3,3iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,short,  3,3sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,double, 3,3dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,int,    4,4iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,short,  4,4sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,float,  4,4fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glVertex,double, 4,4dv)
+
+EIGEN_GL_FUNC_DECLARATION       (glTexCoord)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,int,    2,2iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,short,  2,2sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,float,  2,2fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,double, 2,2dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,int,    3,3iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,short,  3,3sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,double, 3,3dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,int,    4,4iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,short,  4,4sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,float,  4,4fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTexCoord,double, 4,4dv)
+
+EIGEN_GL_FUNC_DECLARATION       (glColor)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,int,    2,2iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,short,  2,2sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,float,  2,2fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,double, 2,2dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,int,    3,3iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,short,  3,3sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,double, 3,3dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,int,    4,4iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,short,  4,4sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,float,  4,4fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glColor,double, 4,4dv)
+
+EIGEN_GL_FUNC_DECLARATION       (glNormal)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glNormal,int,    3,3iv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glNormal,short,  3,3sv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glNormal,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glNormal,double, 3,3dv)
+
+inline void glScale2fv(const float*  v) { glScalef(v[0], v[1], 1.f);  }
+inline void glScale2dv(const double* v) { glScaled(v[0], v[1], 1.0);  }
+inline void glScale3fv(const float*  v) { glScalef(v[0], v[1], v[2]); }
+inline void glScale3dv(const double* v) { glScaled(v[0], v[1], v[2]); }
+
+EIGEN_GL_FUNC_DECLARATION       (glScale)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glScale,float,  2,2fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glScale,double, 2,2dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glScale,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glScale,double, 3,3dv)
+
+inline void glTranslate2fv(const float*  v) { glScalef(v[0], v[1], 1.f);  }
+inline void glTranslate2dv(const double* v) { glScaled(v[0], v[1], 1.0);  }
+inline void glTranslate3fv(const float*  v) { glScalef(v[0], v[1], v[2]); }
+inline void glTranslate3dv(const double* v) { glScaled(v[0], v[1], v[2]); }
+
+EIGEN_GL_FUNC_DECLARATION       (glTranslate)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTranslate,float,  2,2fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTranslate,double, 2,2dv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTranslate,float,  3,3fv)
+EIGEN_GL_FUNC_SPECIALIZATION_VEC(glTranslate,double, 3,3dv)
+
+EIGEN_GL_FUNC_DECLARATION       (glMultMatrix)
+EIGEN_GL_FUNC_SPECIALIZATION_MAT(glMultMatrix,float,  4,4,f)
+EIGEN_GL_FUNC_SPECIALIZATION_MAT(glMultMatrix,double, 4,4,d)
+
+EIGEN_GL_FUNC_DECLARATION       (glLoadMatrix)
+EIGEN_GL_FUNC_SPECIALIZATION_MAT(glLoadMatrix,float,  4,4,f)
+EIGEN_GL_FUNC_SPECIALIZATION_MAT(glLoadMatrix,double, 4,4,d)
+
+void glRotate(const Rotation2D<float>& rot)
+{
+  glRotatef(rot.angle()*180.f/float(M_PI), 0.f, 0.f, 1.f);
+}
+void glRotate(const Rotation2D<double>& rot)
+{
+  glRotated(rot.angle()*180.0/M_PI, 0.0, 0.0, 1.0);
+}
+
+template<typename Derived> void glRotate(const RotationBase<Derived,3>& rot)
+{  
+  Transform<typename Derived::Scalar,3,Projective> tr(rot);
+  glMultMatrix(tr.matrix());
+}
+
+//@}
+
+}
+
+#endif // EIGEN_OPENGL_MODULE
--- a/unsupported/doc/examples/CMakeLists.txt
+++ b/unsupported/doc/examples/CMakeLists.txt
@ -2,6 +2,8 @@ FILE(GLOB examples_SRCS "*.cpp")

 ADD_CUSTOM_TARGET(unsupported_examples)

+INCLUDE_DIRECTORIES(../../../unsupported ../../../unsupported/test)
+
 FOREACH(example_src ${examples_SRCS})
  GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
  ADD_EXECUTABLE(example_${example} ${example_src})
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@ -55,10 +55,10 @@ include_directories(../../test ../../unsupported ../../Eigen)

 if(ADOLC_FOUND)
  include_directories(${ADOLC_INCLUDES})
-  ei_add_property(EIGEN_TESTED_BACKENDS "Adolc")
+  ei_add_property(EIGEN_TESTED_BACKENDS "Adolc, ")
  ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES})
 else(ADOLC_FOUND)
-  ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
+  ei_add_property(EIGEN_MISSING_BACKENDS "Adolc, ")
 endif(ADOLC_FOUND)

 ei_add_test(NonLinearOptimization)
@ -70,6 +70,15 @@ ei_add_test(matrix_function)
 ei_add_test(alignedvector3)
 ei_add_test(FFT)

+find_package(MPFR 2.3.0)
+if(MPFR_FOUND)
+  include_directories(${MPFR_INCLUDES})
+  ei_add_property(EIGEN_TESTED_BACKENDS "MPFR C++, ")
+  ei_add_test(mpreal_support " " ${MPFR_LIBRARIES} )
+else()
+  ei_add_property(EIGEN_MISSING_BACKENDS "MPFR C++, ")
+endif()
+
 ei_add_test(sparse_llt  " " "${SPARSE_LIBS}")
 ei_add_test(sparse_ldlt " " "${SPARSE_LIBS}")
 ei_add_test(sparse_lu   " " "${SPARSE_LIBS}")
@ -77,8 +86,20 @@ ei_add_test(sparse_extra   " " " ")

 find_package(FFTW)
 if(FFTW_FOUND)
+  ei_add_property(EIGEN_TESTED_BACKENDS "fftw, ")
  ei_add_test(FFTW  "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" )
-endif(FFTW_FOUND)
+else()
+  ei_add_property(EIGEN_MISSING_BACKENDS "fftw, ")
+endif()
+
+find_package(OpenGL)
+find_package(GLUT)
+if(OPENGL_FOUND AND GLUT_FOUND)
+  ei_add_property(EIGEN_TESTED_BACKENDS "opengl, ")
+  ei_add_test(openglsupport  "" "${GLUT_LIBRARY}" )
+else()
+  ei_add_property(EIGEN_MISSING_BACKENDS "opengl, ")
+endif()

 find_package(GSL)
 if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
--- a/unsupported/test/mpreal.cpp
+++ b/unsupported/test/mpreal.cpp
@ -0,0 +1,408 @@
+/*
+	Multi-precision real number class. C++ wrapper fo MPFR library.
+	Project homepage: http://www.holoborodko.com/pavel/
+	Contact e-mail:   pavel@holoborodko.com
+
+	Copyright (c) 2008-2010 Pavel Holoborodko
+
+	This library 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 2.1 of the License, or (at your option) any later version.
+
+	This library 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 for more details.
+
+	You should have received a copy of the GNU Lesser General Public
+	License along with this library; if not, write to the Free Software
+	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+
+	Contributors:
+	Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, 
+	Heinz van Saanen, Pere Constans, Dmitriy Gubanov
+*/
+
+#include <cstring>
+#include "mpreal.h"
+
+using std::ws;
+using std::cerr;
+using std::endl;
+using std::string;
+using std::ostream;
+using std::istream;
+
+namespace mpfr{
+
+mp_rnd_t   mpreal::default_rnd  = mpfr_get_default_rounding_mode();	
+mp_prec_t  mpreal::default_prec = mpfr_get_default_prec();	
+int		   mpreal::default_base = 10;
+int        mpreal::double_bits = -1;
+
+// Default constructor: creates mp number and initializes it to 0.
+mpreal::mpreal() 
+{ 
+	mpfr_init2(mp,default_prec); 
+	mpfr_set_ui(mp,0,default_rnd);
+}
+
+mpreal::mpreal(const mpreal& u) 
+{
+	mpfr_init2(mp,mpfr_get_prec(u.mp));
+	mpfr_set(mp,u.mp,default_rnd);
+}
+
+mpreal::mpreal(const mpfr_t u)
+{
+	mpfr_init2(mp,mpfr_get_prec(u));
+	mpfr_set(mp,u,default_rnd);
+}
+
+mpreal::mpreal(const mpf_t u)
+{
+	mpfr_init2(mp,mpf_get_prec(u));
+	mpfr_set_f(mp,u,default_rnd);
+}
+
+mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+	mpfr_init2(mp,prec);
+	mpfr_set_z(mp,u,mode);
+}
+
+mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
+{
+	mpfr_init2(mp,prec);
+	mpfr_set_q(mp,u,mode);
+}
+
+mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
+{
+    if(double_bits == -1 || fits_in_bits(u, double_bits))
+    {
+    	mpfr_init2(mp,prec);
+	    mpfr_set_d(mp,u,mode);
+    }
+    else
+        throw conversion_overflow();
+}
+
+mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+    mpfr_init2(mp,prec);
+	mpfr_set_ld(mp,u,mode);
+}
+
+mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+	mpfr_init2(mp,prec);
+	mpfr_set_ui(mp,u,mode);
+}
+
+mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+	mpfr_init2(mp,prec);
+	mpfr_set_ui(mp,u,mode);
+}
+
+mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+	mpfr_init2(mp,prec);
+	mpfr_set_si(mp,u,mode);
+}
+
+mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
+{ 
+	mpfr_init2(mp,prec);
+	mpfr_set_si(mp,u,mode);
+}
+
+mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
+{
+	mpfr_init2(mp,prec);
+	mpfr_set_str(mp, s, base, mode); 
+}
+
+mpreal::~mpreal() 
+{ 
+	mpfr_clear(mp); 
+}                           
+
+// Operators - Assignment
+mpreal& mpreal::operator=(const char* s)
+{
+	mpfr_t t;
+	if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
+	{
+		mpfr_set(mp,t,mpreal::default_rnd);				
+		mpfr_clear(t);
+	}else{
+		mpfr_clear(t);
+		// cerr<<"fail to convert string"<<endl;
+	}
+
+	return *this;
+}
+
+const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mp_prec_t p1, p2, p3;
+
+	p1 = v1.get_prec(); 
+	p2 = v2.get_prec(); 
+	p3 = v3.get_prec(); 
+
+	a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
+
+	mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
+	return a;
+}
+
+const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mp_prec_t p1, p2, p3;
+
+	p1 = v1.get_prec(); 
+	p2 = v2.get_prec(); 
+	p3 = v3.get_prec(); 
+
+	a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
+
+	mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
+	return a;
+}
+
+const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mp_prec_t p1, p2;
+
+	p1 = v1.get_prec(); 
+	p2 = v2.get_prec(); 
+
+	a.set_prec(p1>p2?p1:p2);
+
+	mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
+
+	return a;
+}
+
+const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mp_prec_t yp, xp;
+
+	yp = y.get_prec(); 
+	xp = x.get_prec(); 
+
+	a.set_prec(yp>xp?yp:xp);
+
+	mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
+
+	return a;
+}
+
+const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
+{
+	mpreal x;
+	mpfr_ptr* t;
+	unsigned long int i;
+
+	t = new mpfr_ptr[n];
+	for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
+	mpfr_sum(x.mp,t,n,rnd_mode);
+	delete[] t;
+	return x;
+}
+
+const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+{	
+	mpreal a;
+	mp_prec_t yp, xp;
+
+	yp = y.get_prec(); 
+	xp = x.get_prec(); 
+
+	a.set_prec(yp>xp?yp:xp);
+
+	mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
+
+	return a;
+}
+
+const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mp_prec_t yp, xp;
+
+	yp = y.get_prec(); 
+	xp = x.get_prec(); 
+
+	a.set_prec(yp>xp?yp:xp);
+
+	mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
+
+	return a;
+}
+
+template <class T>
+std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
+{
+	std::ostringstream oss;
+	oss << f << t;
+	return oss.str();
+}
+
+mpreal::operator std::string() const
+{
+	return to_string();
+}
+
+string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
+{
+	char *s, *ns = NULL;	
+	size_t slen, nslen;
+	mp_exp_t exp;
+	string out;
+
+	if(mpfr_inf_p(mp))
+	{ 
+		if(mpfr_sgn(mp)>0) return "+@Inf@";
+		else			   return "-@Inf@";
+	}
+
+	if(mpfr_zero_p(mp)) return "0";
+	if(mpfr_nan_p(mp))  return "@NaN@";
+
+		
+	s  = mpfr_get_str(NULL,&exp,b,0,mp,mode);
+	ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
+
+	if(s!=NULL && ns!=NULL)
+	{
+		slen  = strlen(s);
+		nslen = strlen(ns);
+		if(nslen<=slen) 
+		{
+			mpfr_free_str(s);
+			s = ns;
+			slen = nslen;
+		}
+		else {
+			mpfr_free_str(ns);
+		}
+
+		// Make human eye-friendly formatting if possible
+		if (exp>0 && static_cast<size_t>(exp)<slen)
+		{
+			if(s[0]=='-')
+			{
+				// Remove zeros starting from right end
+				char* ptr = s+slen-1;
+				while (*ptr=='0' && ptr>s+exp) ptr--; 
+
+				if(ptr==s+exp) out = string(s,exp+1);
+				else		   out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1);
+
+				//out = string(s,exp+1)+'.'+string(s+exp+1);
+			}
+			else
+			{
+				// Remove zeros starting from right end
+				char* ptr = s+slen-1;
+				while (*ptr=='0' && ptr>s+exp-1) ptr--; 
+
+				if(ptr==s+exp-1) out = string(s,exp);
+				else		     out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
+
+				//out = string(s,exp)+'.'+string(s+exp);
+			}
+
+		}else{ // exp<0 || exp>slen
+			if(s[0]=='-')
+			{
+				// Remove zeros starting from right end
+				char* ptr = s+slen-1;
+				while (*ptr=='0' && ptr>s+1) ptr--; 
+
+				if(ptr==s+1) out = string(s,2);
+				else		 out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);
+
+				//out = string(s,2)+'.'+string(s+2);
+			}
+			else
+			{
+				// Remove zeros starting from right end
+				char* ptr = s+slen-1;
+				while (*ptr=='0' && ptr>s) ptr--; 
+
+				if(ptr==s) out = string(s,1);
+				else	   out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);
+
+				//out = string(s,1)+'.'+string(s+1);
+			}
+
+			// Make final string
+			if(--exp)
+			{
+				if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec);
+				else 	  out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
+			}
+		}
+
+		mpfr_free_str(s);
+		return out;
+	}else{
+		return "conversion error!";
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////
+// I/O
+ostream& operator<<(ostream& os, const mpreal& v)
+{
+	return os<<v.to_string(os.precision());
+}
+
+istream& operator>>(istream &is, mpreal& v)
+{
+	char c;	
+	string s = "";
+	mpfr_t t;
+
+	if(is.good())
+	{
+		is>>ws;
+		while ((c = is.get())!=EOF)
+		{
+			if(c ==' ' || c == '\t' || c == '\n' || c == '\r')
+			{
+				is.putback(c);
+				break;
+			}
+			s += c;
+		}
+
+		if(s.size() != 0)
+		{
+			// Protect current value from alternation in case of input error
+			// so some error handling(roll back) procedure can be used 
+			if(0==mpfr_init_set_str(t,s.c_str(),mpreal::default_base,mpreal::default_rnd))
+			{
+				mpfr_set(v.mp,t,mpreal::default_rnd);				
+				mpfr_clear(t);
+
+			}else{
+				mpfr_clear(t);
+				cerr<<"error reading from istream"<<endl;
+				// throw an exception
+			}
+		}
+	}
+	return is;
+}
+}
--- a/unsupported/test/mpreal.h
+++ b/unsupported/test/mpreal.h
--- a/unsupported/test/mpreal_support.cpp
+++ b/unsupported/test/mpreal_support.cpp
@ -0,0 +1,42 @@
+#include "main.h"
+#include <Eigen/MPRealSupport>
+
+using namespace mpfr;
+using namespace std;
+using namespace Eigen;
+
+void test_mpreal_support()
+{
+  // set precision to 256 bits (double has only 53 bits)
+  mpreal::set_default_prec(256);
+  typedef Matrix<mpreal,Eigen::Dynamic,Eigen::Dynamic> MatrixXmp;
+
+  std::cerr << "epsilon =         " << NumTraits<mpreal>::epsilon() << "\n";
+  std::cerr << "dummy_precision = " << NumTraits<mpreal>::dummy_precision() << "\n";
+  std::cerr << "highest =         " << NumTraits<mpreal>::highest() << "\n";
+  std::cerr << "lowest =          " << NumTraits<mpreal>::lowest() << "\n";
+
+  for(int i = 0; i < g_repeat; i++) {
+    int s = ei_random<int>(1,100);
+    MatrixXmp A = MatrixXmp::Random(s,s);
+    MatrixXmp B = MatrixXmp::Random(s,s);
+    MatrixXmp S = A.adjoint() * A;
+    MatrixXmp X;
+
+    // Cholesky
+    X = S.selfadjointView<Lower>().llt().solve(B);
+    VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
+
+    // partial LU
+    X = A.lu().solve(B);
+    VERIFY_IS_APPROX((A*X).eval(),B);
+
+    // symmetric eigenvalues
+    SelfAdjointEigenSolver<MatrixXmp> eig(S);
+    VERIFY_IS_EQUAL(eig.info(), Success);
+    VERIFY_IS_APPROX((S.selfadjointView<Lower>() * eig.eigenvectors()),
+                      eig.eigenvectors() * eig.eigenvalues().asDiagonal());
+  }
+}
+
+#include "mpreal.cpp"
--- a/unsupported/test/openglsupport.cpp
+++ b/unsupported/test/openglsupport.cpp
@ -0,0 +1,61 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include <iostream>
+#include <GL/glut.h>
+#include <Eigen/OpenGLSupport>
+using namespace Eigen;
+
+int main(int argc, char **argv)
+{
+  glutInit(&argc, argv);
+  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
+  glutInitWindowPosition (0,0);
+  glutInitWindowSize(10, 10);
+
+  if(glutCreateWindow("Eigen") <= 0)
+  {
+    std::cerr << "Unable to create GLUT Window.\n";
+    exit(1);
+  }
+
+  Vector3f v3f;
+  Matrix3f rot;
+  glBegin(GL_POINTS);
+  
+  glVertex(v3f);
+  glVertex(2*v3f+v3f);
+  glVertex(rot*v3f);
+  
+  glEnd();
+  
+  Quaterniond qd;
+  glRotate(qd);
+  
+  Matrix4f m44;
+  glLoadMatrix(m44);
+  glMultMatrix(m44);
+  
+  return 0;
+}