matrix storage order can now also be row-dominant (choosable for each matrix separately)

map() moves from MatrixBase to Matrix
much more documentation/examples/snippets

25 files changed, 319 insertions, 67 deletions

diff --git a/Eigen/Core/Conjugate.h b/Eigen/Core/Conjugate.h
index 5254660ac..3fbe29402 100644
--- a/Eigen/Core/Conjugate.h
+++ b/Eigen/Core/Conjugate.h
@@ -26,6 +26,18 @@
 #ifndef EIGEN_CONJUGATE_H
 #define EIGEN_CONJUGATE_H
 
+/** \class Conjugate
+  *
+  * \brief Expression of the complex conjugate of a matrix
+  *
+  * \param MatrixType the type of the object of which we are taking the complex conjugate
+  *
+  * This class represents an expression of the complex conjugate of a matrix.
+  * It is the return type of MatrixBase::conjugate() and is also used by
+  * MatrixBase::adjoint() and most of the time these are the only ways it is used.
+  *
+  * \sa MatrixBase::conjugate(), MatrixBase::adjoint()
+  */
 template<typename MatrixType> class Conjugate : NoOperatorEquals,
   public MatrixBase<typename MatrixType::Scalar, Conjugate<MatrixType> >
 {
@@ -56,6 +68,9 @@ template<typename MatrixType> class Conjugate : NoOperatorEquals,
     MatRef m_matrix;
 };
 
+/** \returns an expression of the complex conjugate of *this.
+  *
+  * \sa adjoint(), class Conjugate */
 template<typename Scalar, typename Derived>
 const Conjugate<Derived>
 MatrixBase<Scalar, Derived>::conjugate() const
@@ -63,6 +78,12 @@ MatrixBase<Scalar, Derived>::conjugate() const
   return Conjugate<Derived>(static_cast<const Derived*>(this)->ref());
 }
 
+/** \returns an expression of the adjoint (i.e. conjugate transpose) of *this.
+  *
+  * Example: \include MatrixBase_adjoint.cpp
+  * Output: \verbinclude MatrixBase_adjoint.out
+  *
+  * \sa transpose(), conjugate(), class Transpose, class Conjugate */
 template<typename Scalar, typename Derived>
 const Transpose<Conjugate<Derived> >
 MatrixBase<Scalar, Derived>::adjoint() const
diff --git a/Eigen/Core/DiagonalCoeffs.h b/Eigen/Core/DiagonalCoeffs.h
index c8a0d43ed..1c5353df0 100644
--- a/Eigen/Core/DiagonalCoeffs.h
+++ b/Eigen/Core/DiagonalCoeffs.h
@@ -26,6 +26,18 @@
 #ifndef EIGEN_DIAGONALCOEFFS_H
 #define EIGEN_DIAGONALCOEFFS_H
 
+/** \class DiagonalCoeffs
+  *
+  * \brief Expression of the main diagonal of a square matrix
+  *
+  * \param MatrixType the type of the object in which we are taking the main diagonal
+  *
+  * This class represents an expression of the main diagonal of a square matrix.
+  * It is the return type of MatrixBase::diagonal() and most of the time this is
+  * the only way it is used.
+  *
+  * \sa MatrixBase::diagonal()
+  */
 template<typename MatrixType> class DiagonalCoeffs
   : public MatrixBase<typename MatrixType::Scalar, DiagonalCoeffs<MatrixType> >
 {
@@ -62,6 +74,12 @@ template<typename MatrixType> class DiagonalCoeffs
     MatRef m_matrix;
 };
 
+/** \returns an expression of the main diagonal of *this, which must be a square matrix.
+  *
+  * Example: \include MatrixBase_diagonal.cpp
+  * Output: \verbinclude MatrixBase_diagonal.out
+  *
+  * \sa class DiagonalCoeffs */
 template<typename Scalar, typename Derived>
 DiagonalCoeffs<Derived>
 MatrixBase<Scalar, Derived>::diagonal()
diff --git a/Eigen/Core/Difference.h b/Eigen/Core/Difference.h
index a03e3c937..815f5b2fb 100644
--- a/Eigen/Core/Difference.h
+++ b/Eigen/Core/Difference.h
@@ -62,6 +62,7 @@ template<typename Lhs, typename Rhs> class Difference : NoOperatorEquals,
     const RhsRef m_rhs;
 };
 
+/** \relates MatrixBase */
 template<typename Scalar, typename Derived1, typename Derived2>
 const Difference<Derived1, Derived2>
 operator-(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
diff --git a/Eigen/Core/Dot.h b/Eigen/Core/Dot.h
index 32cab5bbb..ae64a634a 100644
--- a/Eigen/Core/Dot.h
+++ b/Eigen/Core/Dot.h
@@ -58,6 +58,16 @@ struct DotUnroller<Index, 0, Derived1, Derived2>
   static void run(const Derived1&, const Derived2&, typename Derived1::Scalar&) {}
 };
 
+/** \returns the dot product of *this with other.
+  *
+  * \only_for_vectors
+  * 
+  * \note If the scalar type is complex numbers, then this function returns the hermitian
+  * (sesquilinear) dot product, linear in the first variable and anti-linear in the
+  * second variable.
+  *
+  * \sa norm2(), norm()
+  */
 template<typename Scalar, typename Derived>
 template<typename OtherDerived>
 Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const
@@ -76,18 +86,36 @@ Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const
   return res;
 }
 
+/** \returns the squared norm of *this, i.e. the dot product of *this with itself.
+  *
+  * \only_for_vectors
+  *
+  * \sa dot(), norm()
+  */
 template<typename Scalar, typename Derived>
 typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm2() const
 {
   return real(dot(*this));
 }
 
+/** \returns the norm of *this, i.e. the square root of the dot product of *this with itself.
+  *
+  * \only_for_vectors
+  *
+  * \sa dot(), norm2()
+  */
 template<typename Scalar, typename Derived>
 typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm() const
 {
   return sqrt(norm2());
 }
 
+/** \returns an expression of the quotient of *this by its own norm.
+  *
+  * \only_for_vectors
+  *
+  * \sa norm()
+  */
 template<typename Scalar, typename Derived>
 const ScalarMultiple<typename NumTraits<Scalar>::Real, Derived>
 MatrixBase<Scalar, Derived>::normalized() const
diff --git a/Eigen/Core/Eval.h b/Eigen/Core/Eval.h
index b90018e05..a4bf582f1 100644
--- a/Eigen/Core/Eval.h
+++ b/Eigen/Core/Eval.h
@@ -26,6 +26,25 @@
 #ifndef EIGEN_EVAL_H
 #define EIGEN_EVAL_H
 
+/** \class Eval
+  *
+  * \brief Evaluation of an expression
+  *
+  * The template parameter Expression is the type of the expression that we are evaluating.
+  *
+  * This class is the return
+  * type of MatrixBase::eval() and most of the time this is the only way it
+  * is used.
+  *
+  * However, if you want to write a function returning an evaluation of an expression, you
+  * will need to use this class.
+  *
+  * Here is an example illustrating this:
+  * \include class_Eval.cpp
+  * Output: \verbinclude class_Eval.out
+  *
+  * \sa MatrixBase::eval()
+  */
 template<typename Expression> class Eval : NoOperatorEquals,
   public Matrix< typename Expression::Scalar,
                  Expression::RowsAtCompileTime,
@@ -40,8 +59,21 @@ template<typename Expression> class Eval : NoOperatorEquals,
     Eval(const Expression& expression) : MatrixType(expression) {}
 };
 
+/** Evaluates *this, which can be any expression, and returns the obtained matrix.
+  *
+  * A common use case for this is the following. In an expression-templates library
+  * like Eigen, the coefficients of an expression are only computed as they are
+  * accessed, they are not computed when the expression itself is constructed. This is
+  * usually a good thing, as this "lazy evaluation" improves performance, but can also
+  * in certain cases lead to wrong results and/or to redundant computations. In such
+  * cases, one can restore the classical immediate-evaluation behavior by calling eval().
+  *
+  * Example: \include MatrixBase_eval.cpp
+  * Output: \verbinclude MatrixBase_eval.out
+  *
+  * \sa class Eval */
 template<typename Scalar, typename Derived>
-Eval<Derived> MatrixBase<Scalar, Derived>::eval() const
+const Eval<Derived> MatrixBase<Scalar, Derived>::eval() const
 {
   return Eval<Derived>(*static_cast<const Derived*>(this));
 }
diff --git a/Eigen/Core/Map.h b/Eigen/Core/Map.h
index c11ffcfce..2fcda3e3d 100644
--- a/Eigen/Core/Map.h
+++ b/Eigen/Core/Map.h
@@ -44,18 +44,26 @@ template<typename MatrixType> class Map
     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; }
     int _cols() const { return m_cols; }
     
     const Scalar& _coeff(int row, int col) const
     {
-      return m_data[row + col * m_rows];
+      if(_StorageOrder == ColumnDominant)
+        return m_data[row + col * m_rows];
+      else // RowDominant
+        return m_data[col + row * m_cols];
     }
     
     Scalar& _coeffRef(int row, int col)
     {
-      return const_cast<Scalar*>(m_data)[row + col * m_rows];
+      if(_StorageOrder == ColumnDominant)
+        return const_cast<Scalar*>(m_data)[row + col * m_rows];
+      else // RowDominant
+        return const_cast<Scalar*>(m_data)[col + row * m_cols];
     }
     
   protected:
@@ -63,68 +71,74 @@ template<typename MatrixType> class Map
     int m_rows, m_cols;
 };
 
-template<typename Scalar, typename Derived>
-const Map<Derived> MatrixBase<Scalar, Derived>::map(const Scalar* data, int rows, int cols)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int rows, int cols)
 {
-  return Map<Derived>(data, rows, cols);
+  return Map<Matrix>(data, rows, cols);
 }
 
-template<typename Scalar, typename Derived>
-const Map<Derived> MatrixBase<Scalar, Derived>::map(const Scalar* data, int size)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data, int size)
 {
-  assert(IsVectorAtCompileTime);
-  if(ColsAtCompileTime == 1)
-    return Map<Derived>(data, size, 1);
+  assert(_Cols == 1 || _Rows ==1);
+  if(_Cols == 1)
+    return Map<Matrix>(data, size, 1);
   else
-    return Map<Derived>(data, 1, size);
+    return Map<Matrix>(data, 1, size);
 }
 
-template<typename Scalar, typename Derived>
-const Map<Derived> MatrixBase<Scalar, Derived>::map(const Scalar* data)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+const Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(const Scalar* data)
 {
-  return Map<Derived>(data, RowsAtCompileTime, ColsAtCompileTime);
+  return Map<Matrix>(data, _Rows, _Cols);
 }
 
-template<typename Scalar, typename Derived>
-Map<Derived> MatrixBase<Scalar, Derived>::map(Scalar* data, int rows, int cols)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int rows, int cols)
 {
-  return Map<Derived>(data, rows, cols);
+  return Map<Matrix>(data, rows, cols);
 }
 
-template<typename Scalar, typename Derived>
-Map<Derived> MatrixBase<Scalar, Derived>::map(Scalar* data, int size)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data, int size)
 {
-  assert(IsVectorAtCompileTime);
-  if(ColsAtCompileTime == 1)
-    return Map<Derived>(data, size, 1);
+  assert(_Cols == 1 || _Rows ==1);
+  if(_Cols == 1)
+    return Map<Matrix>(data, size, 1);
   else
-    return Map<Derived>(data, 1, size);
+    return Map<Matrix>(data, 1, size);
 }
 
-template<typename Scalar, typename Derived>
-Map<Derived> MatrixBase<Scalar, Derived>::map(Scalar* data)
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Map<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>::map(Scalar* data)
 {
-  return Map<Derived>(data, RowsAtCompileTime, ColsAtCompileTime);
+  return Map<Matrix>(data, _Rows, _Cols);
 }
 
-template<typename _Scalar, int _Rows, int _Cols>
-Matrix<_Scalar, _Rows, _Cols>
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data, int rows, int cols)
   : Storage(rows, cols)
 {
   *this = map(data, rows, cols);
 }
 
-template<typename _Scalar, int _Rows, int _Cols>
-Matrix<_Scalar, _Rows, _Cols>
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data, int size)
   : Storage(size)
 {
   *this = map(data, size);
 }
 
-template<typename _Scalar, int _Rows, int _Cols>
-Matrix<_Scalar, _Rows, _Cols>
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+Matrix<_Scalar, _Rows, _Cols, _StorageOrder>
   ::Matrix(const Scalar *data)
   : Storage()
 {
diff --git a/Eigen/Core/Matrix.h b/Eigen/Core/Matrix.h
index 73d61a978..812691156 100644
--- a/Eigen/Core/Matrix.h
+++ b/Eigen/Core/Matrix.h
@@ -26,13 +26,14 @@
 #ifndef EIGEN_MATRIX_H
 #define EIGEN_MATRIX_H
 
-template<typename _Scalar, int _Rows, int _Cols>
-class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols> >,
+template<typename _Scalar, int _Rows, int _Cols,
+         MatrixStorageOrder _StorageOrder = EIGEN_DEFAULT_MATRIX_STORAGE_ORDER>
+class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >,
                public MatrixStorage<_Scalar, _Rows, _Cols>
 {
   public:
     friend class MatrixBase<_Scalar, Matrix>;
-    typedef      MatrixBase<_Scalar, Matrix>                Base;
+    typedef      MatrixBase<_Scalar, Matrix>            Base;
     typedef      MatrixStorage<_Scalar, _Rows, _Cols>   Storage;
     typedef      _Scalar                                Scalar;
     typedef      MatrixRef<Matrix>                      Ref;
@@ -44,6 +45,8 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols> >,
     Scalar* data()
     { return Storage::m_data; }
     
+    static const MatrixStorageOrder StorageOrder = _StorageOrder;
+    
   private:
     static const int _RowsAtCompileTime = _Rows, _ColsAtCompileTime = _Cols;
     
@@ -51,12 +54,18 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols> >,
     
     const Scalar& _coeff(int row, int col) const
     {
-      return data()[row + col * Storage::_rows()];
+      if(_StorageOrder == ColumnDominant)
+        return (Storage::m_data)[row + col * Storage::_rows()];
+      else // RowDominant
+        return (Storage::m_data)[col + row * Storage::_cols()];
     }
     
     Scalar& _coeffRef(int row, int col)
     {
-      return data()[row + col * Storage::_rows()];
+      if(_StorageOrder == ColumnDominant)
+        return (Storage::m_data)[row + col * Storage::_rows()];
+      else // RowDominant
+        return (Storage::m_data)[col + row * Storage::_cols()];
     }
     
   public:
@@ -87,6 +96,13 @@ class Matrix : public MatrixBase<_Scalar, Matrix<_Scalar, _Rows, _Cols> >,
     EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, *=)
     EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Matrix, /=)
     
+    static const Map<Matrix> map(const Scalar* array, int rows, int cols);
+    static const Map<Matrix> map(const Scalar* array, int size);
+    static const Map<Matrix> map(const Scalar* array);
+    static Map<Matrix> map(Scalar* array, int rows, int cols);
+    static Map<Matrix> map(Scalar* array, int size);
+    static Map<Matrix> map(Scalar* array);
+    
     explicit Matrix() : Storage()
     {
       assert(_RowsAtCompileTime > 0 && _ColsAtCompileTime > 0);
diff --git a/Eigen/Core/MatrixBase.h b/Eigen/Core/MatrixBase.h
index c15167291..6f9ebe651 100644
--- a/Eigen/Core/MatrixBase.h
+++ b/Eigen/Core/MatrixBase.h
@@ -34,8 +34,10 @@
   * types. Most of the Eigen API is contained in this class.
   *
   * This class takes two template parameters:
-  * \param Scalar the type of the coefficients, e.g. float, double, etc.
-  * \param Derived the derived type, e.g. a matrix type, or an expression, etc.
+  *
+  * \a Scalar is the type of the coefficients, e.g. float, double, etc.
+  *
+  * \a Derived is the derived type, e.g. a matrix type, or an expression, etc.
   * Indeed, a separate MatrixBase type is generated for each derived type
   * so one knows from inside MatrixBase, at compile-time, what the derived type is.
   *
@@ -101,6 +103,8 @@ template<typename Scalar, typename Derived> class MatrixBase
       * \sa rows(), cols(), SizeAtCompileTime. */
     int size() const { return rows() * cols(); }
     /** \returns true if either the number of rows or the number of columns is equal to 1.
+      * In other words, this function returns
+      * \code rows()==1 || cols()==1 \endcode
       * \sa rows(), cols(), IsVectorAtCompileTime. */
     bool isVector() const { return rows()==1 || cols()==1; }
     /** \returns a Ref to *this. \sa Ref */
@@ -151,9 +155,9 @@ template<typename Scalar, typename Derived> class MatrixBase
     RealScalar norm()  const;
     const ScalarMultiple<RealScalar, Derived> normalized() const;
     
-    static Eval<Random<Derived> > random(int rows, int cols);
-    static Eval<Random<Derived> > random(int size);
-    static Eval<Random<Derived> > random();
+    static const Eval<Random<Derived> > random(int rows, int cols);
+    static const Eval<Random<Derived> > random(int size);
+    static const Eval<Random<Derived> > random();
     static const Zero<Derived> zero(int rows, int cols);
     static const Zero<Derived> zero(int size);
     static const Zero<Derived> zero();
@@ -169,13 +173,6 @@ template<typename Scalar, typename Derived> class MatrixBase
     DiagonalCoeffs<Derived> diagonal();
     const DiagonalCoeffs<Derived> diagonal() const;
     
-    static const Map<Derived> map(const Scalar* array, int rows, int cols);
-    static const Map<Derived> map(const Scalar* array, int size);
-    static const Map<Derived> map(const Scalar* array);
-    static Map<Derived> map(Scalar* array, int rows, int cols);
-    static Map<Derived> map(Scalar* array, int size);
-    static Map<Derived> map(Scalar* array);
-    
     template<typename OtherDerived>
     bool isApprox(
       const OtherDerived& other,
@@ -237,7 +234,7 @@ template<typename Scalar, typename Derived> class MatrixBase
     Scalar& z();
     Scalar& w();
 
-    Eval<Derived> eval() const EIGEN_ALWAYS_INLINE;
+    const Eval<Derived> eval() const EIGEN_ALWAYS_INLINE;
 };
 
 #endif // EIGEN_MATRIXBASE_H
diff --git a/Eigen/Core/Minor.h b/Eigen/Core/Minor.h
index 1d2bf178c..ab2e76963 100644
--- a/Eigen/Core/Minor.h
+++ b/Eigen/Core/Minor.h
@@ -26,6 +26,18 @@
 #ifndef EIGEN_MINOR_H
 #define EIGEN_MINOR_H
 
+/** \class Minor
+  *
+  * \brief Expression of a minor
+  *
+  * \param MatrixType the type of the object in which we are taking a minor
+  *
+  * This class represents an expression of a minor. It is the return
+  * type of MatrixBase::minor() and most of the time this is the only way it
+  * is used.
+  *
+  * \sa MatrixBase::minor()
+  */
 template<typename MatrixType> class Minor
   : public MatrixBase<typename MatrixType::Scalar, Minor<MatrixType> >
 {
@@ -75,7 +87,13 @@ template<typename MatrixType> class Minor
 
 /** \return an expression of the (\a row, \a col)-minor of *this,
   * i.e. an expression constructed from *this by removing the specified
-  * row and column. */
+  * row and column.
+  *
+  * Example: \include MatrixBase_minor.cpp
+  * Output: \verbinclude MatrixBase_minor.out
+  *
+  * \sa class Minor
+  */
 template<typename Scalar, typename Derived>
 Minor<Derived>
 MatrixBase<Scalar, Derived>::minor(int row, int col)
diff --git a/Eigen/Core/Product.h b/Eigen/Core/Product.h
index 506df2160..4a89e3a99 100644
--- a/Eigen/Core/Product.h
+++ b/Eigen/Core/Product.h
@@ -115,6 +115,15 @@ MatrixBase<Scalar, Derived>::lazyProduct(const MatrixBase<Scalar, OtherDerived>
   return Product<Derived, OtherDerived>(ref(), other.ref());
 }
 
+/** \relates MatrixBase
+  *
+  * \returns the matrix product of \a mat1 and \a mat2.
+  *
+  * \note This function causes an immediate evaluation. If you want to perform a matrix product
+  * without immediate evaluation, use MatrixBase::lazyProduct() instead.
+  *
+  * \sa MatrixBase::lazyProduct(), MatrixBase::operator*=(const MatrixBase&)
+  */
 template<typename Scalar, typename Derived1, typename Derived2>
 Eval<Product<Derived1, Derived2> >
 operator*(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
diff --git a/Eigen/Core/Random.h b/Eigen/Core/Random.h
index 24ef90655..a8c00c858 100644
--- a/Eigen/Core/Random.h
+++ b/Eigen/Core/Random.h
@@ -56,13 +56,13 @@ template<typename MatrixType> class Random : NoOperatorEquals,
 };
 
 template<typename Scalar, typename Derived>
-Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random(int rows, int cols)
+const Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random(int rows, int cols)
 {
   return Random<Derived>(rows, cols).eval();
 }
 
 template<typename Scalar, typename Derived>
-Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random(int size)
+const Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random(int size)
 {
   assert(IsVectorAtCompileTime);
   if(RowsAtCompileTime == 1) return Random<Derived>(1, size).eval();
@@ -70,7 +70,7 @@ Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random(int size)
 }
 
 template<typename Scalar, typename Derived>
-Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random()
+const Eval<Random<Derived> > MatrixBase<Scalar, Derived>::random()
 {
   return Random<Derived>(RowsAtCompileTime, ColsAtCompileTime).eval();
 }
diff --git a/Eigen/Core/ScalarMultiple.h b/Eigen/Core/ScalarMultiple.h
index ef424f79a..bcaabb383 100644
--- a/Eigen/Core/ScalarMultiple.h
+++ b/Eigen/Core/ScalarMultiple.h
@@ -59,6 +59,7 @@ template<typename FactorType, typename MatrixType> class ScalarMultiple : NoOper
 };
 
 #define EIGEN_MAKE_SCALAR_OPS(FactorType)                              \
+/** \relates MatrixBase */                                             \
 template<typename Scalar, typename Derived>                            \
 const ScalarMultiple<FactorType, Derived>                              \
 operator*(const MatrixBase<Scalar, Derived>& matrix,                   \
@@ -67,6 +68,7 @@ operator*(const MatrixBase<Scalar, Derived>& matrix,                   \
   return ScalarMultiple<FactorType, Derived>(matrix.ref(), scalar);    \
 }                                                                      \
                                                                        \
+/** \relates MatrixBase */                                             \
 template<typename Scalar, typename Derived>                            \
 const ScalarMultiple<FactorType, Derived>                              \
 operator*(FactorType scalar,                                           \
@@ -75,6 +77,7 @@ operator*(FactorType scalar,                                           \
   return ScalarMultiple<FactorType, Derived>(matrix.ref(), scalar);    \
 }                                                                      \
                                                                        \
+/** \relates MatrixBase */                                             \
 template<typename Scalar, typename Derived>                            \
 const ScalarMultiple<typename NumTraits<FactorType>::FloatingPoint, Derived> \
 operator/(const MatrixBase<Scalar, Derived>& matrix,                   \
diff --git a/Eigen/Core/Sum.h b/Eigen/Core/Sum.h
index 3f5f45c7d..6385ced2c 100644
--- a/Eigen/Core/Sum.h
+++ b/Eigen/Core/Sum.h
@@ -61,6 +61,7 @@ template<typename Lhs, typename Rhs> class Sum : NoOperatorEquals,
     const RhsRef m_rhs;
 };
 
+/** \relates MatrixBase */
 template<typename Scalar, typename Derived1, typename Derived2>
 const Sum<Derived1, Derived2>
 operator+(const MatrixBase<Scalar, Derived1> &mat1, const MatrixBase<Scalar, Derived2> &mat2)
diff --git a/Eigen/Core/Trace.h b/Eigen/Core/Trace.h
index 214953ac5..cfe7c7937 100644
--- a/Eigen/Core/Trace.h
+++ b/Eigen/Core/Trace.h
@@ -54,6 +54,9 @@ template<int Index, typename Derived> struct TraceUnroller<Index, 0, Derived>
   static void run(const Derived&, typename Derived::Scalar&) {}
 };
 
+/** \returns the trace of *this, which must be a square matrix.
+  *
+  * \sa diagonal() */
 template<typename Scalar, typename Derived>
 Scalar MatrixBase<Scalar, Derived>::trace() const
 {
diff --git a/Eigen/Core/Transpose.h b/Eigen/Core/Transpose.h
index c3723ad82..fb79ea3cd 100644
--- a/Eigen/Core/Transpose.h
+++ b/Eigen/Core/Transpose.h
@@ -26,6 +26,18 @@
 #ifndef EIGEN_TRANSPOSE_H
 #define EIGEN_TRANSPOSE_H
 
+/** \class Transpose
+  *
+  * \brief Expression of the transpose of a matrix
+  *
+  * \param MatrixType the type of the object of which we are taking the transpose
+  *
+  * This class represents an expression of the transpose of a matrix.
+  * It is the return type of MatrixBase::transpose() and MatrixBase::adjoint()
+  * and most of the time this is the only way it is used.
+  *
+  * \sa MatrixBase::transpose(), MatrixBase::adjoint()
+  */
 template<typename MatrixType> class Transpose
   : public MatrixBase<typename MatrixType::Scalar, Transpose<MatrixType> >
 {
@@ -63,6 +75,12 @@ template<typename MatrixType> class Transpose
     MatRef m_matrix;
 };
 
+/** \returns an expression of the transpose of *this.
+  *
+  * Example: \include MatrixBase_transpose.cpp
+  * Output: \verbinclude MatrixBase_transpose.out
+  *
+  * \sa adjoint(), class DiagonalCoeffs */
 template<typename Scalar, typename Derived>
 Transpose<Derived>
 MatrixBase<Scalar, Derived>::transpose()
@@ -70,7 +88,7 @@ MatrixBase<Scalar, Derived>::transpose()
   return Transpose<Derived>(ref());
 }
 
-/** This is the const version of transpose(). */
+/** This is the const version of transpose(). \sa adjoint() */
 template<typename Scalar, typename Derived>
 const Transpose<Derived>
 MatrixBase<Scalar, Derived>::transpose() const
diff --git a/Eigen/Core/Util.h b/Eigen/Core/Util.h
index 0bccc1f75..a38feaf9a 100644
--- a/Eigen/Core/Util.h
+++ b/Eigen/Core/Util.h
@@ -32,6 +32,10 @@
 #define EIGEN_UNROLLED_LOOPS (true)
 #endif
 
+#ifndef EIGEN_DEFAULT_MATRIX_STORAGE_ORDER
+#define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER ColumnDominant
+#endif
+
 #undef minor
 
 #define USING_PART_OF_NAMESPACE_EIGEN \
@@ -82,8 +86,17 @@ EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \
 EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
 EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
 
+const int Dynamic = -1;
+
+enum MatrixStorageOrder
+{
+  ColumnDominant,
+  RowDominant
+};
+
 //forward declarations
-template<typename _Scalar, int _Rows, int _Cols> class Matrix;
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+  class Matrix;
 template<typename MatrixType> class MatrixRef;
 template<typename NewScalar, typename MatrixType> class Cast;
 template<typename MatrixType> class Row;
@@ -112,14 +125,12 @@ template<typename T> struct ForwardDecl
   typedef T Ref;
 };
 
-template<typename _Scalar, int _Rows, int _Cols>
-struct ForwardDecl<Matrix<_Scalar, _Rows, _Cols> >
+template<typename _Scalar, int _Rows, int _Cols, MatrixStorageOrder _StorageOrder>
+struct ForwardDecl<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> >
 {
-  typedef MatrixRef<Matrix<_Scalar, _Rows, _Cols> > Ref;
+  typedef MatrixRef<Matrix<_Scalar, _Rows, _Cols, _StorageOrder> > Ref;
 };
 
-const int Dynamic = -1;
-
 //classes inheriting NoOperatorEquals don't generate a default operator=.
 class NoOperatorEquals
 {
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 721c3de81..8b9aa62af 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -237,13 +237,13 @@ PERLMOD_MAKEVAR_PREFIX =
 # Configuration options related to the preprocessor   
 #---------------------------------------------------------------------------
 ENABLE_PREPROCESSING   = YES
-MACRO_EXPANSION        = NO
-EXPAND_ONLY_PREDEF     = NO
+MACRO_EXPANSION        = YES
+EXPAND_ONLY_PREDEF     = YES
 SEARCH_INCLUDES        = YES
 INCLUDE_PATH           = 
 INCLUDE_FILE_PATTERNS  = 
 PREDEFINED             = 
-EXPAND_AS_DEFINED      = 
+EXPAND_AS_DEFINED      = EIGEN_MAKE_SCALAR_OPS
 SKIP_FUNCTION_MACROS   = YES
 #---------------------------------------------------------------------------
 # Configuration::additions related to external references   
diff --git a/doc/examples/class_Eval.cpp b/doc/examples/class_Eval.cpp
new file mode 100644
index 000000000..bcf66c316
--- /dev/null
+++ b/doc/examples/class_Eval.cpp
@@ -0,0 +1,28 @@
+#include <Eigen/Core.h>
+USING_PART_OF_NAMESPACE_EIGEN
+using namespace std;
+
+template<typename Scalar, typename Derived>
+const Eigen::Eval<Eigen::Transpose<Derived> >
+evaluatedTranspose(const MatrixBase<Scalar, Derived>& m)
+{
+  return m.transpose().eval();
+}
+
+int main(int, char**)
+{
+  Matrix2f M = Matrix2f::random();
+  Matrix2f m;
+  m = M;
+  cout << "Here is the matrix m:" << endl << m << endl;
+  cout << "Now we want to replace m by its own transpose." << endl;
+  cout << "If we do m = m.transpose(), then m becomes:" << endl;
+  m = m.transpose();
+  cout << m << endl << "which is wrong!" << endl;
+  cout << "Now let us instead do m = evaluatedTranspose(m). Then m becomes" << endl;
+  m = M;
+  m = evaluatedTranspose(m);
+  cout << m << endl << "which is right." << endl;
+
+  return 0;
+}
diff --git a/doc/snippets/MatrixBase_adjoint.cpp b/doc/snippets/MatrixBase_adjoint.cpp
new file mode 100644
index 000000000..661803ca1
--- /dev/null
+++ b/doc/snippets/MatrixBase_adjoint.cpp
@@ -0,0 +1,3 @@
+Matrix2cf m = Matrix2cf::random();
+cout << "Here is the 2x2 complex matrix m:" << endl << m << endl;
+cout << "Here is the adjoint of m:" << endl << m.adjoint() << endl;
diff --git a/doc/snippets/MatrixBase_block.cpp b/doc/snippets/MatrixBase_block.cpp
index 0a1b6d6d3..d8a5ecc71 100644
--- a/doc/snippets/MatrixBase_block.cpp
+++ b/doc/snippets/MatrixBase_block.cpp
@@ -1,3 +1,5 @@
 Matrix4d m = Matrix4d::diagonal(Vector4d(1,2,3,4));
+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);
-cout << m << endl;
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_diagonal.cpp b/doc/snippets/MatrixBase_diagonal.cpp
new file mode 100644
index 000000000..005fba7b6
--- /dev/null
+++ b/doc/snippets/MatrixBase_diagonal.cpp
@@ -0,0 +1,4 @@
+Matrix3i m = Matrix3i::random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here are the coefficients on the main diagonal of m:" << endl
+     << m.diagonal() << endl;
diff --git a/doc/snippets/MatrixBase_dynBlock.cpp b/doc/snippets/MatrixBase_dynBlock.cpp
index 4fc6aceb3..edea29da6 100644
--- a/doc/snippets/MatrixBase_dynBlock.cpp
+++ b/doc/snippets/MatrixBase_dynBlock.cpp
@@ -1,3 +1,5 @@
 Matrix3d m = Matrix3d::diagonal(Vector3d(1,2,3));
+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);
-cout << m << endl;
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_eval.cpp b/doc/snippets/MatrixBase_eval.cpp
new file mode 100644
index 000000000..376957359
--- /dev/null
+++ b/doc/snippets/MatrixBase_eval.cpp
@@ -0,0 +1,12 @@
+Matrix2f M = Matrix2f::random();
+Matrix2f m;
+m = M;
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Now we want to replace m by its own transpose." << endl;
+cout << "If we do m = m.transpose(), then m becomes:" << endl;
+m = m.transpose();
+cout << m << endl << "which is wrong!" << endl;
+cout << "Now let us instead do m = m.transpose().eval(). Then m becomes" << endl;
+m = M;
+m = m.transpose().eval();
+cout << m << endl << "which is right." << endl;
diff --git a/doc/snippets/MatrixBase_minor.cpp b/doc/snippets/MatrixBase_minor.cpp
new file mode 100644
index 000000000..5f9dc5c96
--- /dev/null
+++ b/doc/snippets/MatrixBase_minor.cpp
@@ -0,0 +1,3 @@
+Matrix3i m = Matrix3i::random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.minor(1,1):" << endl << m.minor(1,1) << endl;
diff --git a/doc/snippets/MatrixBase_transpose.cpp b/doc/snippets/MatrixBase_transpose.cpp
new file mode 100644
index 000000000..0a0586d53
--- /dev/null
+++ b/doc/snippets/MatrixBase_transpose.cpp
@@ -0,0 +1,8 @@
+Matrix2f m = Matrix2f::random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the transpose of m:" << endl << m.transpose() << endl;
+cout << "Here is the coefficient (1,0) in the transpose of m:" << endl
+     << m.transpose()(1,0) << endl;
+cout << "Let us overwrite this coefficient with the value 0." << endl;
+m.transpose()(1,0) = 0;
+cout << "Now the matrix m is:" << endl << m << endl;