- make MatrixBase and all expressions aware of their preferred traversal order.

  Honor this preference in operator=.
- add several methods to the API
- rework API for diagonal matrices
- add benchmarking code

33 files changed, 354 insertions, 84 deletions

diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index 8456b2b78..19552308e 100644
--- a/Eigen/src/Core/Block.h
+++ b/Eigen/src/Core/Block.h
@@ -67,9 +67,10 @@ template<typename MatrixType, int BlockRows, int BlockCols> class Block
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = BlockRows,
                      _ColsAtCompileTime = BlockCols;
-
+    
     const Block& _ref() const { return *this; }
     int _rows() const { return BlockRows; }
     int _cols() const { return BlockCols; }
diff --git a/Eigen/src/Core/Cast.h b/Eigen/src/Core/Cast.h
index 513a070a9..126370487 100644
--- a/Eigen/src/Core/Cast.h
+++ b/Eigen/src/Core/Cast.h
@@ -57,6 +57,7 @@ template<typename NewScalar, typename MatrixType> class Cast : NoOperatorEquals,
     Cast(const MatRef& matrix) : m_matrix(matrix) {}
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
     const Cast& _ref() const { return *this; }
diff --git a/Eigen/src/Core/Column.h b/Eigen/src/Core/Column.h
index 30d8ae684..1ee63423b 100644
--- a/Eigen/src/Core/Column.h
+++ b/Eigen/src/Core/Column.h
@@ -63,6 +63,7 @@ template<typename MatrixType> class Column
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Column)
     
   private:
+    static const TraversalOrder _Order = ColumnMajor;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = 1;
     const Column& _ref() const { return *this; }
diff --git a/Eigen/src/Core/Conjugate.h b/Eigen/src/Core/Conjugate.h
index 00797c2ff..2a900a8e9 100644
--- a/Eigen/src/Core/Conjugate.h
+++ b/Eigen/src/Core/Conjugate.h
@@ -49,6 +49,7 @@ template<typename MatrixType> class Conjugate : NoOperatorEquals,
     Conjugate(const MatRef& matrix) : m_matrix(matrix) {}
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/DiagonalCoeffs.h b/Eigen/src/Core/DiagonalCoeffs.h
index 1b9920cdd..44904fbcd 100644
--- a/Eigen/src/Core/DiagonalCoeffs.h
+++ b/Eigen/src/Core/DiagonalCoeffs.h
@@ -51,6 +51,7 @@ template<typename MatrixType> class DiagonalCoeffs
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(DiagonalCoeffs)
     
   private:
+    static const TraversalOrder _Order = ColumnMajor;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = 1;
     
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 125fa58d9..67f10cfb7 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -26,27 +26,39 @@
 #ifndef EIGEN_DIAGONALMATRIX_H
 #define EIGEN_DIAGONALMATRIX_H
 
-template<typename MatrixType, typename CoeffsVectorType>
+/** \class DiagonalMatrix
+  *
+  * \brief Expression of a diagonal matrix
+  *
+  * \param CoeffsVectorType the type of the vector of diagonal coefficients
+  *
+  * This class is an expression of a diagonal matrix with given vector of diagonal
+  * coefficients. It is the return
+  * type of MatrixBase::diagonal(const OtherDerived&) and most of the time this is
+  * the only way it is used.
+  *
+  * \sa MatrixBase::diagonal(const OtherDerived&)
+  */
+template<typename CoeffsVectorType>
 class DiagonalMatrix : NoOperatorEquals,
-  public MatrixBase<typename MatrixType::Scalar,
-                    DiagonalMatrix<MatrixType, CoeffsVectorType> >
+  public MatrixBase<typename CoeffsVectorType::Scalar,
+                    DiagonalMatrix<CoeffsVectorType> >
 {
   public:
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename CoeffsVectorType::Scalar Scalar;
     typedef typename CoeffsVectorType::Ref CoeffsVecRef;
-    friend class MatrixBase<Scalar, DiagonalMatrix<MatrixType, CoeffsVectorType> >;
+    friend class MatrixBase<Scalar, DiagonalMatrix<CoeffsVectorType> >;
     
     DiagonalMatrix(const CoeffsVecRef& coeffs) : m_coeffs(coeffs)
     {
       assert(CoeffsVectorType::IsVectorAtCompileTime
-          && _RowsAtCompileTime == _ColsAtCompileTime
-          && _RowsAtCompileTime == CoeffsVectorType::SizeAtCompileTime
           && coeffs.size() > 0);
     }
     
   private:
-    static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-                     _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
+    static const TraversalOrder _Order = Indifferent;
+    static const int _RowsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
+                     _ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime;
 
     const DiagonalMatrix& _ref() const { return *this; }
     int _rows() const { return m_coeffs.size(); }
@@ -61,12 +73,20 @@ class DiagonalMatrix : NoOperatorEquals,
     CoeffsVecRef m_coeffs;
 };
 
+/** \returns an expression of a diagonal matrix with *this as vector of diagonal coefficients
+  *
+  * \only_for_vectors
+  *
+  * Example: \include MatrixBase_asDiagonal.cpp
+  * Output: \verbinclude MatrixBase_asDiagonal.out
+  *
+  * \sa class DiagonalMatrix
+  **/
 template<typename Scalar, typename Derived>
-template<typename OtherDerived>
-const DiagonalMatrix<Derived, OtherDerived>
-MatrixBase<Scalar, Derived>::diagonal(const OtherDerived& coeffs)
+const DiagonalMatrix<Derived>
+MatrixBase<Scalar, Derived>::asDiagonal() const
 {
-  return DiagonalMatrix<Derived, OtherDerived>(coeffs);
+  return DiagonalMatrix<Derived>(ref());
 }
 
 #endif // EIGEN_DIAGONALMATRIX_H
diff --git a/Eigen/src/Core/Difference.h b/Eigen/src/Core/Difference.h
index 7ed3bae4d..b249e4e98 100644
--- a/Eigen/src/Core/Difference.h
+++ b/Eigen/src/Core/Difference.h
@@ -42,6 +42,7 @@ template<typename Lhs, typename Rhs> class Difference : NoOperatorEquals,
     }
 
   private:
+    static const TraversalOrder _Order = Lhs::Order;
     static const int _RowsAtCompileTime = Lhs::RowsAtCompileTime,
                      _ColsAtCompileTime = Rhs::ColsAtCompileTime;
     
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index ae64a634a..322ff63a2 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -123,4 +123,27 @@ MatrixBase<Scalar, Derived>::normalized() const
   return (*this) / norm();
 }
 
+template<typename Scalar, typename Derived>
+template<typename OtherDerived>
+bool MatrixBase<Scalar, Derived>::isOrtho
+(const OtherDerived& other,
+ const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()) const
+{
+  return abs2(dot(other)) <= prec * prec * norm2() * other.norm2();
+}
+
+template<typename Scalar, typename Derived>
+bool MatrixBase<Scalar, Derived>::isOrtho
+(const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()) const
+{
+  for(int i = 0; i < cols(); i++)
+  {
+    if(!isApprox(col(i).norm2(), static_cast<Scalar>(1)))
+      return false;
+    for(int j = 0; j < i; j++)
+      if(!isMuchSmallerThan(col(i).dot(col(j)), static_cast<Scalar>(1)))
+        return false;
+  }
+  return true;
+}
 #endif // EIGEN_DOT_H
diff --git a/Eigen/src/Core/DynBlock.h b/Eigen/src/Core/DynBlock.h
index b46bed0eb..e47b09fb1 100644
--- a/Eigen/src/Core/DynBlock.h
+++ b/Eigen/src/Core/DynBlock.h
@@ -67,6 +67,7 @@ template<typename MatrixType> class DynBlock
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(DynBlock)
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int
       _RowsAtCompileTime = MatrixType::RowsAtCompileTime == 1 ? 1 : Dynamic,
       _ColsAtCompileTime = MatrixType::ColsAtCompileTime == 1 ? 1 : Dynamic;
diff --git a/Eigen/src/Core/Fuzzy.h b/Eigen/src/Core/Fuzzy.h
index bd3219568..6cf076cdc 100644
--- a/Eigen/src/Core/Fuzzy.h
+++ b/Eigen/src/Core/Fuzzy.h
@@ -30,7 +30,7 @@ template<typename Scalar, typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Scalar, Derived>::isApprox(
   const OtherDerived& other,
-  const typename NumTraits<Scalar>::Real& prec
+  const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
 ) const
 {
   assert(rows() == other.rows() && cols() == other.cols());
@@ -51,7 +51,7 @@ bool MatrixBase<Scalar, Derived>::isApprox(
 template<typename Scalar, typename Derived>
 bool MatrixBase<Scalar, Derived>::isMuchSmallerThan(
   const typename NumTraits<Scalar>::Real& other,
-  const typename NumTraits<Scalar>::Real& prec
+  const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
 ) const
 {
   if(IsVectorAtCompileTime)
@@ -71,7 +71,7 @@ template<typename Scalar, typename Derived>
 template<typename OtherDerived>
 bool MatrixBase<Scalar, Derived>::isMuchSmallerThan(
   const MatrixBase<Scalar, OtherDerived>& other,
-  const typename NumTraits<Scalar>::Real& prec
+  const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
 ) const
 {
   assert(rows() == other.rows() && cols() == other.cols());
diff --git a/Eigen/src/Core/Identity.h b/Eigen/src/Core/Identity.h
index a24a30224..068493125 100644
--- a/Eigen/src/Core/Identity.h
+++ b/Eigen/src/Core/Identity.h
@@ -39,6 +39,7 @@ template<typename MatrixType> class Identity : NoOperatorEquals,
     }
     
   private:
+    static const TraversalOrder _Order = Indifferent;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
@@ -61,4 +62,20 @@ const Identity<Derived> MatrixBase<Scalar, Derived>::identity(int rows)
   return Identity<Derived>(rows);
 }
 
+template<typename Scalar, typename Derived>
+bool MatrixBase<Scalar, Derived>::isIdentity
+(const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()) const
+{
+  for(int j = 0; j < col(); j++)
+  {
+    if(!isApprox(coeff(j, j), static_cast<Scalar>(1)))
+      return false;
+    for(int i = 0; i < j; i++)
+      if(!isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1)))
+        return false;
+  }
+  return true;
+}
+
+
 #endif // EIGEN_IDENTITY_H
diff --git a/Eigen/src/Core/Map.h b/Eigen/src/Core/Map.h
index 2e3a4838f..fa7a75a49 100644
--- a/Eigen/src/Core/Map.h
+++ b/Eigen/src/Core/Map.h
@@ -47,10 +47,10 @@ template<typename MatrixType> class Map
     friend class MatrixBase<Scalar, Map<MatrixType> >;
 
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
-    static const MatrixStorageOrder _StorageOrder = MatrixType::StorageOrder;
 
     const Map& _ref() const { return *this; }
     int _rows() const { return m_rows; }
@@ -58,17 +58,17 @@ template<typename MatrixType> class Map
     
     const Scalar& _coeff(int row, int col) const
     {
-      if(_StorageOrder == ColumnDominant)
+      if(_Order == ColumnMajor)
         return m_data[row + col * m_rows];
-      else // RowDominant
+      else // RowMajor
         return m_data[col + row * m_cols];
     }
     
     Scalar& _coeffRef(int row, int col)
     {
-      if(_StorageOrder == ColumnDominant)
+      if(_Order == ColumnMajor)
         return const_cast<Scalar*>(m_data)[row + col * m_rows];
-      else // RowDominant
+      else // RowMajor
         return const_cast<Scalar*>(m_data)[col + row * m_cols];
     }
   
@@ -89,7 +89,7 @@ template<typename MatrixType> class Map
 };
 
 /** This is the const version of map(Scalar*,int,int). */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int rows, int cols)
 {
@@ -97,7 +97,7 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int rows,
 }
 
 /** This is the const version of map(Scalar*,int). */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int size)
 {
@@ -109,7 +109,7 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int size)
 }
 
 /** This is the const version of map(Scalar*). */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data)
 {
@@ -124,8 +124,10 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data)
   *
   * Example: \include MatrixBase_map_int_int.cpp
   * Output: \verbinclude MatrixBase_map_int_int.out
+  *
+  * \sa map(const Scalar*, int, int), map(Scalar*, int), map(Scalar*), class Map
   */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int rows, int cols)
 {
@@ -135,14 +137,16 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int rows, int co
 /** \returns a expression of a vector mapping the given data.
   *
   * \param data The array of data to map
-  * \param rows The size (number of coefficients) of the expression to construct
+  * \param size The size (number of coefficients) of the expression to construct
   *
   * \only_for_vectors
   *
   * Example: \include MatrixBase_map_int.cpp
   * Output: \verbinclude MatrixBase_map_int.out
+  *
+  * \sa map(const Scalar*, int), map(Scalar*, int, int), map(Scalar*), class Map
   */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int size)
 {
@@ -159,15 +163,25 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int size)
   *
   * Example: \include MatrixBase_map.cpp
   * Output: \verbinclude MatrixBase_map.out
+  *
+  * \sa map(const Scalar*), map(Scalar*, int), map(Scalar*, int, int), class Map
   */
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data)
 {
   return Map<Matrix>(data, _Rows, _Cols);
 }
 
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+/** Constructor copying an existing array of data. Only useful for dynamic-size matrices:
+  * for fixed-size matrices, it is redundant to pass the \a rows and \a cols parameters.
+  * \param data The array of data to copy
+  * \param rows The number of rows of the matrix to construct
+  * \param cols The number of columns of the matrix to construct
+  *
+  * \sa Matrix(const Scalar *), Matrix::map(const Scalar *, int, int)
+  */
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data, int rows, int cols)
   : Storage(rows, cols)
@@ -175,7 +189,17 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   *this = map(data, rows, cols);
 }
 
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+/** Constructor copying an existing array of data. Only useful for dynamic-size vectors:
+  * for fixed-size vectors, it is redundant to pass the \a size parameter.
+  *
+  * \only_for_vectors
+  *
+  * \param data The array of data to copy
+  * \param size The size of the vector to construct
+  *
+  * \sa Matrix(const Scalar *), Matrix::map(const Scalar *, int)
+  */
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data, int size)
   : Storage(size)
@@ -183,7 +207,17 @@ Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   *this = map(data, size);
 }
 
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+/** Constructor copying an existing array of data.
+  * Only for fixed-size matrices and vectors.
+  * \param data The array of data to copy
+  *
+  * For dynamic-size matrices and vectors, see the variants taking additional int parameters
+  * for the dimensions.
+  *
+  * \sa Matrix(const Scalar *, int), Matrix(const Scalar *, int, int),
+  * Matrix::map(const Scalar *)
+  */
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data)
   : Storage()
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 149d2cdca..8148e2ac6 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -26,9 +26,54 @@
 #ifndef EIGEN_MATRIX_H
 #define EIGEN_MATRIX_H
 
-/** \class Matrix */
+/** \class Matrix
+  *
+  * \brief The matrix class, also used for vectors and row-vectors
+  *
+  * \param _Scalar the scalar type, i.e. the type of the coefficients
+  * \param _Rows the number of rows at compile-time. Use the special value \a Dynamic to specify that the number of rows is dynamic, i.e. is not fixed at compile-time.
+  * \param _Cols the number of columns at compile-time. Use the special value \a Dynamic to specify that the number of columns is dynamic, i.e. is not fixed at compile-time.
+  * \param _StorageOrder can be either \a RowMajor or \a ColumnMajor.
+  * This template parameter has a default value (EIGEN_DEFAULT_MATRIX_STORAGE_ORDER)
+  * which, if not predefined, is defined to \a ColumnMajor. You can override this behavior by
+  * predefining it before including Eigen headers.
+  *
+  * This single class template covers all kinds of matrix and vectors that Eigen can handle.
+  * All matrix and vector types are just typedefs to specializations of this class template.
+  *
+  * These typedefs are as follows:
+  * \li \c %Matrix##Size##Type for square matrices
+  * \li \c Vector##Size##Type for vectors (matrices with one column)
+  * \li \c RowVector##Size##Type for row-vectors (matrices with one row)
+  *
+  * where \c Size can be
+  * \li \c 2 for fixed size 2
+  * \li \c 3 for fixed size 3
+  * \li \c 4 for fixed size 4
+  * \li \c X for dynamic size
+  *
+  * and \c Type can be
+  * \li \c i for type \c int
+  * \li \c f for type \c float
+  * \li \c d for type \c double
+  * \li \c cf for type \c std::complex<float>
+  * \li \c cd for type \c std::complex<float>
+  *
+  * Examples:
+  * \li \c Matrix2d is a typedef for \c Matrix<double,2,2>
+  * \li \c VectorXf is a typedef for \c Matrix<float,Dynamic,1>
+  * \li \c RowVector3i is a typedef for \c Matrix<int,1,3>
+  *
+  * Of course these typedefs do not exhaust all the possibilities offered by the Matrix class
+  * template, they only address some of the most common cases. For instance, if you want a
+  * fixed-size matrix with 3 rows and 5 columns, there is no typedef for that, so you should use
+  * \c Matrix<double,3,5>.
+  *
+  * Note that most of the API is in the base class MatrixBase, and that the base class
+  * MatrixStorage also provides the MatrixStorage::resize() public method.
+  */
 template<typename _Scalar, int _Rows, int _Cols,
-         MatrixStorageOrder _StorageOrder = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER>
+         TraversalOrder _StorageOrder = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER>
 class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >,
                public MatrixStorage<_Scalar, _Rows, _Cols>
 {
@@ -46,30 +91,37 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
     Scalar* data()
     { return Storage::m_data; }
     
-    static const MatrixStorageOrder StorageOrder = _StorageOrder;
-    
   private:
+    static const TraversalOrder _Order = _StorageOrder;
     static const int _RowsAtCompileTime = _Rows, _ColsAtCompileTime = _Cols;
     
     Ref _ref() const { return Ref(*this); }
     
     const Scalar& _coeff(int row, int col) const
     {
-      if(_StorageOrder == ColumnDominant)
+      if(_Order == ColumnMajor)
         return (Storage::m_data)[row + col * Storage::_rows()];
-      else // RowDominant
+      else // RowMajor
         return (Storage::m_data)[col + row * Storage::_cols()];
     }
     
     Scalar& _coeffRef(int row, int col)
     {
-      if(_StorageOrder == ColumnDominant)
+      if(_Order == ColumnMajor)
         return (Storage::m_data)[row + col * Storage::_rows()];
-      else // RowDominant
+      else // RowMajor
         return (Storage::m_data)[col + row * Storage::_cols()];
     }
     
   public:
+    /** Copies the value of the expression \a other into *this.
+      *
+      * *this is resized (if possible) to match the dimensions of \a other.
+      *
+      * As a special exception, copying a row-vector into a vector (and conversely)
+      * is allowed. The resizing, if any, is then done in the appropriate way so that
+      * row-vectors remain row-vectors and vectors remain vectors.
+      */
     template<typename OtherDerived> 
     Matrix& operator=(const MatrixBase<Scalar, OtherDerived>& other)
     {
@@ -87,6 +139,9 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
       return Base::operator=(other);
     }
     
+    /** This is a special case of the templated operator=. Its purpose is to
+      * prevent a default operator= from hiding the templated operator=.
+      */
     Matrix& operator=(const Matrix& other)
     {
       return operator=<Matrix>(other);
@@ -104,10 +159,21 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
     static Map<Matrix> map(Scalar* array, int size);
     static Map<Matrix> map(Scalar* array);
     
+    /** Default constructor, does nothing. Only for fixed-size matrices.
+      * For dynamic-size matrices and vectors, this constructor is forbidden (guarded by
+      * an assertion) because it would leave the matrix without an allocated data buffer.
+      */
     explicit Matrix() : Storage()
     {
       assert(_RowsAtCompileTime > 0 && _ColsAtCompileTime > 0);
     }
+    
+    /** Constructs a vector or row-vector with given dimension. \only_for_vectors
+      *
+      * Note that this is only useful for dynamic-size vectors. For fixed-size vectors,
+      * it is redundant to pass the dimension here, so it makes more sense to use the default
+      * constructor Matrix() instead.
+      */
     explicit Matrix(int dim) : Storage(dim)
     {
       assert(dim > 0);
@@ -117,14 +183,16 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
               && (_RowsAtCompileTime == Dynamic || _RowsAtCompileTime == dim)));
     }
     
-    // this constructor is very tricky.
-    // When Matrix is a fixed-size vector type of size 2,
-    // Matrix(x,y) should mean "construct vector with coefficients x,y".
-    // Otherwise, Matrix(x,y) should mean "construct matrix with x rows and y cols".
-    // Note that in the case of fixed-size, Storage::Storage(int,int) does nothing,
-    // so it is harmless to call it and afterwards we just fill the m_data array
-    // with the two coefficients. In the case of dynamic size, Storage::Storage(int,int)
-    // does what we want to, so it only remains to add some asserts.
+    /** This constructor has two very different behaviors, depending on the type of *this.
+      *
+      * \li When Matrix is a fixed-size vector type of size 2, this constructor constructs
+      *     an initialized vector. The parameters \a x, \a y are copied into the first and second
+      *     coords of the vector respectively.
+      * \li Otherwise, this constructor constructs an uninitialized matrix with \a x rows and
+      *     \a y columns. This is useful for dynamic-size matrices. For fixed-size matrices,
+      *     it is redundant to pass these parameters, so one should use the default constructor
+      *     Matrix() instead.
+      */
     Matrix(int x, int y) : Storage(x, y)
     {
       if((_RowsAtCompileTime == 1 && _ColsAtCompileTime == 2)
@@ -139,6 +207,7 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
             && y > 0 && (_ColsAtCompileTime == Dynamic || _ColsAtCompileTime == y));
       }
     }
+    /** constructs an initialized 2D vector with given coefficients */
     Matrix(const float& x, const float& y)
     {
       assert((_RowsAtCompileTime == 1 && _ColsAtCompileTime == 2)
@@ -146,6 +215,7 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
       (Storage::m_data)[0] = x;
       (Storage::m_data)[1] = y;
     }
+    /** constructs an initialized 2D vector with given coefficients */
     Matrix(const double& x, const double& y)
     {
       assert((_RowsAtCompileTime == 1 && _ColsAtCompileTime == 2)
@@ -153,6 +223,7 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
       (Storage::m_data)[0] = x;
       (Storage::m_data)[1] = y;
     }
+    /** constructs an initialized 3D vector with given coefficients */
     Matrix(const Scalar& x, const Scalar& y, const Scalar& z)
     {
       assert((_RowsAtCompileTime == 1 && _ColsAtCompileTime == 3)
@@ -161,6 +232,7 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
       (Storage::m_data)[1] = y;
       (Storage::m_data)[2] = z;
     }
+    /** constructs an initialized 4D vector with given coefficients */
     Matrix(const Scalar& x, const Scalar& y, const Scalar& z, const Scalar& w)
     {
       assert((_RowsAtCompileTime == 1 && _ColsAtCompileTime == 4)
@@ -174,16 +246,19 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _Storage
     Matrix(const Scalar *data, int size);
     explicit Matrix(const Scalar *data);
     
+    /** Constructor copying the value of the expression \a other */
     template<typename OtherDerived>
     Matrix(const MatrixBase<Scalar, OtherDerived>& other)
              : Storage(other.rows(), other.cols())
     {
       *this = other;
     }
+    /** Copy constructor */
     Matrix(const Matrix& other) : Storage(other.rows(), other.cols())
     {
       *this = other;
     }
+    /** Destructor */
     ~Matrix() {}
 };
 
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 6f9ebe651..8f556388f 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -57,6 +57,8 @@
 template<typename Scalar, typename Derived> class MatrixBase
 {
   public:
+    static const TraversalOrder Order = Derived::_Order;
+  
     /** The number of rows at compile-time. This is just a copy of the value provided
       * by the \a Derived type. If a value is not known at compile-time,
       * it is set to the \a Dynamic constant.
@@ -154,6 +156,10 @@ template<typename Scalar, typename Derived> class MatrixBase
     RealScalar norm2() const;
     RealScalar norm()  const;
     const ScalarMultiple<RealScalar, Derived> normalized() const;
+    template<typename OtherDerived>
+    bool isOrtho(const OtherDerived& other,
+                 const typename NumTraits<Scalar>::Real& prec) const;
+    bool isOrtho(const typename NumTraits<Scalar>::Real& prec) const;
     
     static const Eval<Random<Derived> > random(int rows, int cols);
     static const Eval<Random<Derived> > random(int size);
@@ -166,9 +172,11 @@ template<typename Scalar, typename Derived> class MatrixBase
     static const Ones<Derived> ones();
     static const Identity<Derived> identity(int rows = RowsAtCompileTime);
     
-    template<typename OtherDerived>
-    static const DiagonalMatrix<Derived, OtherDerived>
-    diagonal(const OtherDerived& coeffs);
+    bool isZero(const typename NumTraits<Scalar>::Real& prec) const;
+    bool isOnes(const typename NumTraits<Scalar>::Real& prec) const;
+    bool isIdentity(const typename NumTraits<Scalar>::Real& prec) const;
+    
+    const DiagonalMatrix<Derived> asDiagonal() const;
     
     DiagonalCoeffs<Derived> diagonal();
     const DiagonalCoeffs<Derived> diagonal() const;
@@ -176,16 +184,16 @@ template<typename Scalar, typename Derived> class MatrixBase
     template<typename OtherDerived>
     bool isApprox(
       const OtherDerived& other,
-      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
+      const typename NumTraits<Scalar>::Real& prec
     ) const;
     bool isMuchSmallerThan(
       const typename NumTraits<Scalar>::Real& other,
-      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
+      const typename NumTraits<Scalar>::Real& prec
     ) const;
     template<typename OtherDerived>
     bool isMuchSmallerThan(
       const MatrixBase<Scalar, OtherDerived>& other,
-      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
+      const typename NumTraits<Scalar>::Real& prec
     ) const;
     
     template<typename OtherDerived>
diff --git a/Eigen/src/Core/MatrixRef.h b/Eigen/src/Core/MatrixRef.h
index 835782b32..57ae4a492 100644
--- a/Eigen/src/Core/MatrixRef.h
+++ b/Eigen/src/Core/MatrixRef.h
@@ -39,6 +39,7 @@ template<typename MatrixType> class MatrixRef
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(MatrixRef)
 
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
     
diff --git a/Eigen/src/Core/Minor.h b/Eigen/src/Core/Minor.h
index b4fedebe0..11d47d4ac 100644
--- a/Eigen/src/Core/Minor.h
+++ b/Eigen/src/Core/Minor.h
@@ -57,6 +57,7 @@ template<typename MatrixType> class Minor
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Minor)
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int
       _RowsAtCompileTime = (MatrixType::RowsAtCompileTime != Dynamic) ?
                            MatrixType::RowsAtCompileTime - 1 : Dynamic,
diff --git a/Eigen/src/Core/Ones.h b/Eigen/src/Core/Ones.h
index cf421f686..cea830dc8 100644
--- a/Eigen/src/Core/Ones.h
+++ b/Eigen/src/Core/Ones.h
@@ -40,6 +40,7 @@ template<typename MatrixType> class Ones : NoOperatorEquals,
     friend class MatrixBase<Scalar, Ones<MatrixType> >;
   
   private:
+    static const TraversalOrder _Order = Indifferent;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
@@ -125,4 +126,15 @@ const Ones<Derived> MatrixBase<Scalar, Derived>::ones()
   return Ones<Derived>(RowsAtCompileTime, ColsAtCompileTime);
 }
 
+template<typename Scalar, typename Derived>
+bool MatrixBase<Scalar, Derived>::isOnes
+(const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()) const
+{
+  for(int j = 0; j < col(); j++)
+    for(int i = 0; i < row(); i++)
+      if(!isApprox(coeff(i, j), static_cast<Scalar>(1)))
+        return false;
+  return true;
+}
+
 #endif // EIGEN_ONES_H
diff --git a/Eigen/src/Core/OperatorEquals.h b/Eigen/src/Core/OperatorEquals.h
index 30d42b670..e2cc5d559 100644
--- a/Eigen/src/Core/OperatorEquals.h
+++ b/Eigen/src/Core/OperatorEquals.h
@@ -27,28 +27,25 @@
 #ifndef EIGEN_OPERATOREQUALS_H
 #define EIGEN_OPERATOREQUALS_H
 
-template<typename Derived1, typename Derived2, int UnrollCount, int Rows>
+template<typename Derived1, typename Derived2, int UnrollCount, TraversalOrder Order>
 struct MatrixOperatorEqualsUnroller
 {
-  static const int col = (UnrollCount-1) / Rows;
-  static const int row = (UnrollCount-1) % Rows;
+  static const int col = (Order == ColumnMajor)
+                       ? (UnrollCount-1) / Derived1::RowsAtCompileTime
+                       : (UnrollCount-1) % Derived1::ColsAtCompileTime;
+  static const int row = (Order == ColumnMajor)
+                       ? (UnrollCount-1) % Derived1::RowsAtCompileTime
+                       : (UnrollCount-1) / Derived1::ColsAtCompileTime;
 
   static void run(Derived1 &dst, const Derived2 &src)
   {
-    MatrixOperatorEqualsUnroller<Derived1, Derived2, UnrollCount-1, Rows>::run(dst, src);
+    MatrixOperatorEqualsUnroller<Derived1, Derived2, UnrollCount-1, Order>::run(dst, src);
     dst.coeffRef(row, col) = src.coeff(row, col);
   }
 };
 
-// prevent buggy user code from causing an infinite recursion
-template<typename Derived1, typename Derived2, int UnrollCount>
-struct MatrixOperatorEqualsUnroller<Derived1, Derived2, UnrollCount, 0>
-{
-  static void run(Derived1 &, const Derived2 &) {}
-};
-
-template<typename Derived1, typename Derived2, int Rows>
-struct MatrixOperatorEqualsUnroller<Derived1, Derived2, 1, Rows>
+template<typename Derived1, typename Derived2, TraversalOrder Order>
+struct MatrixOperatorEqualsUnroller<Derived1, Derived2, 1, Order>
 {
   static void run(Derived1 &dst, const Derived2 &src)
   {
@@ -56,8 +53,15 @@ struct MatrixOperatorEqualsUnroller<Derived1, Derived2, 1, Rows>
   }
 };
 
-template<typename Derived1, typename Derived2, int Rows>
-struct MatrixOperatorEqualsUnroller<Derived1, Derived2, Dynamic, Rows>
+// prevent buggy user code from causing an infinite recursion
+template<typename Derived1, typename Derived2, TraversalOrder Order>
+struct MatrixOperatorEqualsUnroller<Derived1, Derived2, 0, Order>
+{
+  static void run(Derived1 &, const Derived2 &) {}
+};
+
+template<typename Derived1, typename Derived2, TraversalOrder Order>
+struct MatrixOperatorEqualsUnroller<Derived1, Derived2, Dynamic, Order>
 {
   static void run(Derived1 &, const Derived2 &) {}
 };
@@ -101,7 +105,8 @@ template<typename OtherDerived>
 Derived& MatrixBase<Scalar, Derived>
   ::operator=(const MatrixBase<Scalar, OtherDerived>& other)
 {
-  if(IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime) // copying a vector expression into a vector
+  if(IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime)
+    // copying a vector expression into a vector
   {
     assert(size() == other.size());
     if(EIGEN_UNROLLED_LOOPS && SizeAtCompileTime != Dynamic && SizeAtCompileTime <= 25)
@@ -113,17 +118,24 @@ Derived& MatrixBase<Scalar, Derived>
         coeffRef(i) = other.coeff(i);
     return *static_cast<Derived*>(this);
   }
-  else // all other cases (typically, but not necessarily, copying a matrix)
+  else // copying a matrix expression into a matrix
   {
     assert(rows() == other.rows() && cols() == other.cols());
     if(EIGEN_UNROLLED_LOOPS && SizeAtCompileTime != Dynamic && SizeAtCompileTime <= 25)
       MatrixOperatorEqualsUnroller
-        <Derived, OtherDerived, SizeAtCompileTime, RowsAtCompileTime>::run
+        <Derived, OtherDerived, SizeAtCompileTime, Order>::run
           (*static_cast<Derived*>(this), *static_cast<const OtherDerived*>(&other));
     else
-      for(int j = 0; j < cols(); j++) //traverse in column-dominant order
+    {
+      if(Order == ColumnMajor)
+        for(int j = 0; j < cols(); j++)
+          for(int i = 0; i < rows(); i++)
+            coeffRef(i, j) = other.coeff(i, j);
+      else // RowMajor
         for(int i = 0; i < rows(); i++)
-          coeffRef(i, j) = other.coeff(i, j);
+          for(int j = 0; j < cols(); j++)
+            coeffRef(i, j) = other.coeff(i, j);
+    }
     return *static_cast<Derived*>(this);
   }
 }
diff --git a/Eigen/src/Core/Opposite.h b/Eigen/src/Core/Opposite.h
index 110b63cf2..e0b6c1ac3 100644
--- a/Eigen/src/Core/Opposite.h
+++ b/Eigen/src/Core/Opposite.h
@@ -37,6 +37,7 @@ template<typename MatrixType> class Opposite : NoOperatorEquals,
     Opposite(const MatRef& matrix) : m_matrix(matrix) {}
     
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/Random.h b/Eigen/src/Core/Random.h
index ef07bfa43..1b870af74 100644
--- a/Eigen/src/Core/Random.h
+++ b/Eigen/src/Core/Random.h
@@ -40,6 +40,7 @@ template<typename MatrixType> class Random : NoOperatorEquals,
     friend class MatrixBase<Scalar, Random<MatrixType> >;
   
   private:
+    static const TraversalOrder _Order = Indifferent;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/Row.h b/Eigen/src/Core/Row.h
index 4387b2e6c..40f53de4c 100644
--- a/Eigen/src/Core/Row.h
+++ b/Eigen/src/Core/Row.h
@@ -69,6 +69,7 @@ template<typename MatrixType> class Row
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Row)
     
   private:
+    static const TraversalOrder _Order = RowMajor;
     static const int _RowsAtCompileTime = 1,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/ScalarMultiple.h b/Eigen/src/Core/ScalarMultiple.h
index 29c566977..76a4477cc 100644
--- a/Eigen/src/Core/ScalarMultiple.h
+++ b/Eigen/src/Core/ScalarMultiple.h
@@ -38,6 +38,7 @@ template<typename FactorType, typename MatrixType> class ScalarMultiple : NoOper
       : m_matrix(matrix), m_factor(factor) {}
 
   private:
+    static const TraversalOrder _Order = MatrixType::Order;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/Sum.h b/Eigen/src/Core/Sum.h
index 51dc66bfa..453033673 100644
--- a/Eigen/src/Core/Sum.h
+++ b/Eigen/src/Core/Sum.h
@@ -42,6 +42,7 @@ template<typename Lhs, typename Rhs> class Sum : NoOperatorEquals,
     }
 
   private:
+    static const TraversalOrder _Order = Lhs::Order;
     static const int _RowsAtCompileTime = Lhs::RowsAtCompileTime,
                      _ColsAtCompileTime = Rhs::ColsAtCompileTime;
 
diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index 777e86143..dd0f8cc3e 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -51,6 +51,8 @@ template<typename MatrixType> class Transpose
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose)
     
   private:
+    static const TraversalOrder _Order = (MatrixType::Order == ColumnMajor)
+                                       ? RowMajor : ColumnMajor;
     static const int _RowsAtCompileTime = MatrixType::ColsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::RowsAtCompileTime;
 
diff --git a/Eigen/src/Core/Util.h b/Eigen/src/Core/Util.h
index a38feaf9a..0202f37c9 100644
--- a/Eigen/src/Core/Util.h
+++ b/Eigen/src/Core/Util.h
@@ -33,7 +33,7 @@
 #endif
 
 #ifndef EIGEN_DEFAULT_MATRIX_STORAGE_ORDER
-#define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER ColumnDominant
+#define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER Indifferent
 #endif
 
 #undef minor
@@ -88,14 +88,15 @@ EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
 
 const int Dynamic = -1;
 
-enum MatrixStorageOrder
+enum TraversalOrder
 {
-  ColumnDominant,
-  RowDominant
+  ColumnMajor,
+  RowMajor,
+  Indifferent = ColumnMajor
 };
 
 //forward declarations
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
   class Matrix;
 template<typename MatrixType> class MatrixRef;
 template<typename NewScalar, typename MatrixType> class Cast;
@@ -114,7 +115,7 @@ template<typename FactorType, typename MatrixType> class ScalarMultiple;
 template<typename MatrixType> class Random;
 template<typename MatrixType> class Zero;
 template<typename MatrixType> class Ones;
-template<typename MatrixType, typename CoeffsVectorType> class DiagonalMatrix;
+template<typename CoeffsVectorType> class DiagonalMatrix;
 template<typename MatrixType> class DiagonalCoeffs;
 template<typename MatrixType> class Identity;
 template<typename ExpressionType> class Eval;
@@ -125,7 +126,7 @@ template<typename T> struct ForwardDecl
   typedef T Ref;
 };
 
-template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+template<typename _Scalar, int _Rows, int _Cols, TraversalOrder _StorageOrder>
 struct ForwardDecl<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 {
   typedef MatrixRef<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> > Ref;
diff --git a/Eigen/src/Core/Zero.h b/Eigen/src/Core/Zero.h
index 533e1e0f1..56ada635a 100644
--- a/Eigen/src/Core/Zero.h
+++ b/Eigen/src/Core/Zero.h
@@ -40,6 +40,7 @@ template<typename MatrixType> class Zero : NoOperatorEquals,
     friend class MatrixBase<Scalar, Zero<MatrixType> >;
   
   private:
+    static const TraversalOrder _Order = Indifferent;
     static const int _RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      _ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
@@ -125,4 +126,15 @@ const Zero<Derived> MatrixBase<Scalar, Derived>::zero()
   return Zero<Derived>(RowsAtCompileTime, ColsAtCompileTime);
 }
 
+template<typename Scalar, typename Derived>
+bool MatrixBase<Scalar, Derived>::isZero
+(const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()) const
+{
+  for(int j = 0; j < col(); j++)
+    for(int i = 0; i < row(); i++)
+      if(!isMuchSmallerThan(coeff(i, j), static_cast<Scalar>(1)))
+        return false;
+  return true;
+}
+
 #endif // EIGEN_ZERO_H
diff --git a/doc/benchmarkX.cpp b/doc/benchmarkX.cpp
new file mode 100644
index 000000000..1670341e9
--- /dev/null
+++ b/doc/benchmarkX.cpp
@@ -0,0 +1,22 @@
+// g++ -O3 -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX
+
+#include <Eigen/Core>
+
+using namespace std;
+USING_PART_OF_NAMESPACE_EIGEN
+
+int main(int argc, char *argv[])
+{
+	MatrixXd I = MatrixXd::identity(20);
+	MatrixXd m(20,20);
+	for(int i = 0; i < 20; i++) for(int j = 0; j < 20; j++)
+	{
+		m(i,j) = 0.1 * (i+20*j);
+	}
+	for(int a = 0; a < 1000000; a++)
+	{
+		m = I + 0.00005 * (m + m*m);
+	}
+	cout << m << endl;
+	return 0;
+}
diff --git a/doc/benchmark_suite b/doc/benchmark_suite
new file mode 100755
index 000000000..b75a8df67
--- /dev/null
+++ b/doc/benchmark_suite
@@ -0,0 +1,16 @@
+echo "Fixed size 3x3, ColumnMajor, -DNDEBUG"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=ColumnMajor -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+echo "Fixed size 3x3, ColumnMajor, with asserts"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=ColumnMajor benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+echo "Fixed size 3x3, RowMajor, -DNDEBUG"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=RowMajor -DNDEBUG benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+echo "Fixed size 3x3, RowMajor, with asserts"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=RowMajor benchmark.cpp -o benchmark && time ./benchmark >/dev/null
+echo "Dynamic size 20x20, ColumnMajor, -DNDEBUG"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=ColumnMajor -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+echo "Dynamic size 20x20, ColumnMajor, with asserts"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=ColumnMajor benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+echo "Dynamic size 20x20, RowMajor, -DNDEBUG"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=RowMajor -DNDEBUG benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
+echo "Dynamic size 20x20, RowMajor, with asserts"
+g++ -O3 -I .. -DEIGEN_DEFAULT_MATRIX_STORAGE_ORDER=RowMajor benchmarkX.cpp -o benchmarkX && time ./benchmarkX >/dev/null
diff --git a/doc/snippets/MatrixBase_asDiagonal.cpp b/doc/snippets/MatrixBase_asDiagonal.cpp
new file mode 100644
index 000000000..4cbff8231
--- /dev/null
+++ b/doc/snippets/MatrixBase_asDiagonal.cpp
@@ -0,0 +1 @@
+cout << Vector3i(2,5,6).asDiagonal() << endl;
diff --git a/doc/snippets/MatrixBase_block.cpp b/doc/snippets/MatrixBase_block.cpp
index d8a5ecc71..6bdc1ab03 100644
--- a/doc/snippets/MatrixBase_block.cpp
+++ b/doc/snippets/MatrixBase_block.cpp
@@ -1,4 +1,4 @@
-Matrix4d m = Matrix4d::diagonal(Vector4d(1,2,3,4));
+Matrix4d m = Vector4d(1,2,3,4).asDiagonal();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is m.block<2, 2>(2, 2):" << endl << m.block<2, 2>(2, 2) << endl;
 m.block<2, 2>(2, 0) = m.block<2, 2>(2, 2);
diff --git a/doc/snippets/MatrixBase_dynBlock.cpp b/doc/snippets/MatrixBase_dynBlock.cpp
index edea29da6..c04db2140 100644
--- a/doc/snippets/MatrixBase_dynBlock.cpp
+++ b/doc/snippets/MatrixBase_dynBlock.cpp
@@ -1,4 +1,4 @@
-Matrix3d m = Matrix3d::diagonal(Vector3d(1,2,3));
+Matrix3d m = Vector3d(1,2,3).asDiagonal();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is m.dynBlock(1, 1, 2, 1):" << endl << m.dynBlock(1, 1, 2, 1) << endl;
 m.dynBlock(1, 0, 2, 1) = m.dynBlock(1, 1, 2, 1);
diff --git a/test/main.h b/test/main.h
index c0c56ffa6..460000511 100644
--- a/test/main.h
+++ b/test/main.h
@@ -28,6 +28,7 @@
 
 #include <QtTest/QtTest>
 
+//#define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER RowMajor
 #define EIGEN_INTERNAL_DEBUGGING
 #include <Eigen/Core>
 
diff --git a/test/miscmatrices.cpp b/test/miscmatrices.cpp
index 052d5ec91..7126266fc 100644
--- a/test/miscmatrices.cpp
+++ b/test/miscmatrices.cpp
@@ -46,7 +46,7 @@ template<typename MatrixType> void miscMatrices(const MatrixType& m)
   VectorType v1 = VectorType::random(rows);
   v1[0];
   Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
-  square = MatrixType::diagonal(v1);
+  square = v1.asDiagonal();
   if(r==r2) VERIFY_IS_APPROX(square(r,r2), v1[r]);
   else VERIFY_IS_MUCH_SMALLER_THAN(square(r,r2), static_cast<Scalar>(1));
   square = MatrixType::zero(rows, rows);