diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h
index 1c66cde8e..f2bde0e71 100644
--- a/Eigen/src/Core/PermutationMatrix.h
+++ b/Eigen/src/Core/PermutationMatrix.h
@@ -25,18 +25,24 @@
 #ifndef EIGEN_PERMUTATIONMATRIX_H
 #define EIGEN_PERMUTATIONMATRIX_H
 
-/** \nonstableyet
-  * \class PermutationMatrix
+/** \class PermutationMatrix
   *
   * \brief Permutation matrix
   *
   * \param SizeAtCompileTime the number of rows/cols, or Dynamic
-  * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime.
+  * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
   *
   * This class represents a permutation matrix, internally stored as a vector of integers.
-  * The convention followed here is the same as on <a href="http://en.wikipedia.org/wiki/Permutation_matrix">Wikipedia</a>,
-  * namely: the matrix of permutation \a p is the matrix such that on each row \a i, the only nonzero coefficient is
-  * in column p(i).
+  * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
+  * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
+  *  \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
+  * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
+  *  \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
+  *
+  * Permutation matrices are square and invertible.
+  *
+  * Notice that in addition to the member functions and operators listed here, there also are non-member
+  * operator* to multiply a PermutationMatrix with any kind of matrix expression (MatrixBase) on either side.
   *
   * \sa class DiagonalMatrix
   */
@@ -53,6 +59,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
 {
   public:
 
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
     typedef ei_traits<PermutationMatrix> Traits;
     typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime>
             DenseMatrixType;
@@ -65,25 +72,37 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
     };
     typedef typename Traits::Scalar Scalar;
+    #endif
 
-    typedef Matrix<int, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> IndicesType;
+    typedef Matrix<int, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
 
     inline PermutationMatrix()
     {
     }
 
+    /** Copy constructor. */
     template<int OtherSize, int OtherMaxSize>
     inline PermutationMatrix(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
       : m_indices(other.indices()) {}
 
-    /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** Standard copy constructor. Defined only to prevent a default copy constructor
+      * from hiding the other templated constructor */
     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
+    #endif
 
-    /** generic constructor from expression of the indices */
+    /** Generic constructor from expression of the indices. The indices
+      * array has the meaning that the permutations sends each integer i to indices[i].
+      *
+      * \warning It is your responsibility to check that the indices array that you passes actually
+      * describes a permutation, \ie each value between 0 and n-1 occurs exactly once, where n is the
+      * array's size.
+      */
     template<typename Other>
-    explicit inline PermutationMatrix(const MatrixBase<Other>& other) : m_indices(other)
+    explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
     {}
 
+    /** Copies the other permutation into *this */
     template<int OtherSize, int OtherMaxSize>
     PermutationMatrix& operator=(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
     {
@@ -91,6 +110,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
       return *this;
     }
 
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** This is a special case of the templated operator=. Its purpose is to
       * prevent a default operator= from hiding the templated operator=.
       */
@@ -99,8 +119,12 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
       m_indices = other.m_indices();
       return *this;
     }
+    #endif
 
-    inline PermutationMatrix(int rows, int cols) : m_indices(rows)
+    /** Constructs an uninitialized permutation matrix of given size. Note that it is required
+      * that rows == cols, since permutation matrices are square. The \a cols parameter may be omitted.
+      */
+    inline PermutationMatrix(int rows, int cols = rows) : m_indices(rows)
     {
       ei_assert(rows == cols);
     }
@@ -111,24 +135,33 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
     /** \returns the number of columns */
     inline int cols() const { return m_indices.size(); }
 
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename DenseDerived>
     void evalTo(MatrixBase<DenseDerived>& other) const
     {
       other.setZero();
       for (int i=0; i<rows();++i)
-        other.coeffRef(i,m_indices.coeff(i)) = typename DenseDerived::Scalar(1);
+        other.coeffRef(m_indices.coeff(i),i) = typename DenseDerived::Scalar(1);
     }
+    #endif
 
+    /** \returns a Matrix object initialized from this permutation matrix. Notice that it
+      * is inefficient to return this Matrix object by value. For efficiency, favor using
+      * the Matrix constructor taking AnyMatrixBase objects.
+      */
     DenseMatrixType toDenseMatrix() const
     {
       return *this;
     }
 
+    /** const version of indices(). */
     const IndicesType& indices() const { return m_indices; }
+    /** \returns a reference to the stored array representing the permutation. */
     IndicesType& indices() { return m_indices; }
 
     /**** inversion and multiplication helpers to hopefully get RVO ****/
 
+#ifndef EIGEN_PARSED_BY_DOXYGEN
   protected:
     enum Inverse_t {Inverse};
     PermutationMatrix(Inverse_t, const PermutationMatrix& other)
@@ -143,10 +176,19 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
       ei_assert(lhs.cols() == rhs.rows());
       for (int i=0; i<rows();++i) m_indices.coeffRef(i) = lhs.m_indices.coeff(rhs.m_indices.coeff(i));
     }
+#endif
 
   public:
+    /** \returns the inverse permutation matrix.
+      *
+      * \note \note_try_to_help_rvo
+      */
     inline PermutationMatrix inverse() const
     { return PermutationMatrix(Inverse, *this); }
+    /** \returns the product permutation matrix.
+      *
+      * \note \note_try_to_help_rvo
+      */
     template<int OtherSize, int OtherMaxSize>
     inline PermutationMatrix operator*(const PermutationMatrix<OtherSize, OtherMaxSize>& other) const
     { return PermutationMatrix(Product, *this, other); }
@@ -209,7 +251,7 @@ struct ei_permut_matrix_product_retval
           Dest,
           Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime,
           Side==OnTheRight ? 1 : Dest::ColsAtCompileTime
-        >(dst, Side==OnTheRight ? m_permutation.indices().coeff(i) : i)
+        >(dst, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i)
 
         =
 
@@ -217,7 +259,7 @@ struct ei_permut_matrix_product_retval
           MatrixTypeNestedCleaned,
           Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,
           Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime
-        >(m_matrix, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i);
+        >(m_matrix, Side==OnTheRight ? m_permutation.indices().coeff(i) : i);
       }
     }
 
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index da99d3592..82ebf6aa0 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -218,7 +218,8 @@ ALIASES                = "only_for_vectors=This is only for vectors (either row-
                          "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \
                          "note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \
                          "note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \
-                         "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values."
+                         "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values." \
+                         "note_try_to_help_rvo=This function returns the result by value. In order to make that efficient, it is implemented as just a return statement using a special constructor, hopefully allowing the compiler to perform a RVO (return value optimization)."
 
 # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
 # sources only. Doxygen will then generate output that is more tailored for C.
diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp
index ec3a8541c..c4affc795 100644
--- a/test/permutationmatrices.cpp
+++ b/test/permutationmatrices.cpp
@@ -66,7 +66,7 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
 
   for (int i=0; i<rows; i++)
     for (int j=0; j<cols; j++)
-        VERIFY_IS_APPROX(m_original(lv(i),j), m_permuted(i,rv(j)));
+        VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j)));
 
   Matrix<Scalar,Rows,Rows> lm(lp);
   Matrix<Scalar,Cols,Cols> rm(rp);
@@ -80,7 +80,7 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
   randomPermutationVector(lv2, rows);
   LeftPermutationType lp2(lv2);
   Matrix<Scalar,Rows,Rows> lm2(lp2);
-  VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm2*lm);
+  VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
 }
 
 void test_permutationmatrices()