PermutationMatrix: add setIdentity and transpositions methods

LU: make use of that

Eigen/src/Core/PermutationMatrix.h

@@ -116,18 +116,15 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
       */
     PermutationMatrix& operator=(const PermutationMatrix& other)
     {
-      m_indices = other.m_indices();
+      m_indices = other.m_indices;
       return *this;
     }
     #endif
 
-    /** 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.
+    /** Constructs an uninitialized permutation matrix of given size.
       */
-    inline PermutationMatrix(int rows, int cols = rows) : m_indices(rows)
-    {
-      ei_assert(rows == cols);
-    }
+    inline PermutationMatrix(int size) : m_indices(size)
+    {}
 
     /** \returns the number of rows */
     inline int rows() const { return m_indices.size(); }
@@ -159,8 +156,62 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi
     /** \returns a reference to the stored array representing the permutation. */
     IndicesType& indices() { return m_indices; }
 
-    /** Resizes to given size. */
-    inline void resize(int size) { m_indices.resize(size); }
+    /** Resizes to given size.
+      */
+    inline void resize(int size)
+    {
+      m_indices.resize(size);
+    }
+
+    /** Sets *this to be the identity permutation matrix */
+    void setIdentity()
+    {
+      for(int i = 0; i < m_indices.size(); ++i)
+        m_indices.coeffRef(i) = i;
+    }
+
+    /** Sets *this to be the identity permutation matrix of given size.
+      */
+    void setIdentity(int size)
+    {
+      resize(size);
+      setIdentity();
+    }
+
+    /** Multiplies *this by the transposition \f$(ij)\f$ on the left.
+      *
+      * \returns a reference to *this.
+      *
+      * \warning This is much slower than applyTranspositionOnTheRight(int,int):
+      * this has linear complexity and requires a lot of branching.
+      *
+      * \sa applyTranspositionOnTheRight(int,int)
+      */
+    PermutationMatrix& applyTranspositionOnTheLeft(int i, int j)
+    {
+      ei_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size());
+      for(int k = 0; k < m_indices.size(); ++k)
+      {
+        if(m_indices.coeff(k) == i) m_indices.coeffRef(k) = j;
+        else if(m_indices.coeff(k) == j) m_indices.coeffRef(k) = i;
+      }
+      return *this;
+    }
+
+    /** Multiplies *this by the transposition \f$(ij)\f$ on the right.
+      *
+      * \returns a reference to *this.
+      *
+      * This is a fast operation, it only consists in swapping two indices.
+      *
+      * \sa applyTranspositionOnTheLeft(int,int)
+      */
+    PermutationMatrix& applyTranspositionOnTheRight(int i, int j)
+    {
+      ei_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size());
+      std::swap(m_indices.coeffRef(i), m_indices.coeffRef(j));
+      return *this;
+    }
 
     /**** inversion and multiplication helpers to hopefully get RVO ****/
 

Eigen/src/LU/FullPivLU.h

@@ -389,8 +389,6 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
 {
   m_isInitialized = true;
   m_lu = matrix;
-  m_p.resize(matrix.rows());
-  m_q.resize(matrix.cols());
 
   const int size = matrix.diagonalSize();
   const int rows = matrix.rows();
@@ -459,13 +457,13 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
   // the main loop is over, we still have to accumulate the transpositions to find the
   // permutations P and Q
 
-  for(int k = 0; k < matrix.rows(); ++k) m_p.indices().coeffRef(k) = k;
+  m_p.setIdentity(rows);
   for(int k = size-1; k >= 0; --k)
-    std::swap(m_p.indices().coeffRef(k), m_p.indices().coeffRef(rows_transpositions.coeff(k)));
+    m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
 
-  for(int k = 0; k < matrix.cols(); ++k) m_q.indices().coeffRef(k) = k;
+  m_q.setIdentity(cols);
   for(int k = 0; k < size; ++k)
-    std::swap(m_q.indices().coeffRef(k), m_q.indices().coeffRef(cols_transpositions.coeff(k)));
+    m_q.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   return *this;

Eigen/src/LU/PartialPivLU.h

@@ -375,7 +375,6 @@ template<typename MatrixType>
 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
 {
   m_lu = matrix;
-  m_p.resize(matrix.rows());
 
   ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
   const int size = matrix.rows();
@@ -386,9 +385,9 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma
   ei_partial_lu_inplace(m_lu, rows_transpositions, nb_transpositions);
   m_det_p = (nb_transpositions%2) ? -1 : 1;
 
-  for(int k = 0; k < size; ++k) m_p.indices().coeffRef(k) = k;
+  m_p.setIdentity(size);
   for(int k = size-1; k >= 0; --k)
-    std::swap(m_p.indices().coeffRef(k), m_p.indices().coeffRef(rows_transpositions.coeff(k)));
+    m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
 
   m_isInitialized = true;
   return *this;

test/permutationmatrices.cpp

@@ -81,6 +81,30 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
   LeftPermutationType lp2(lv2);
   Matrix<Scalar,Rows,Rows> lm2(lp2);
   VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
+
+  LeftPermutationType identityp;
+  identityp.setIdentity(rows);
+  VERIFY_IS_APPROX(m_original, identityp*m_original);
+
+  if(rows>1 && cols>1)
+  {
+    lp2 = lp;
+    int i = ei_random<int>(0, rows-1);
+    int j;
+    do j = ei_random<int>(0, rows-1); while(j==i);
+    lp2.applyTranspositionOnTheLeft(i, j);
+    lm = lp;
+    lm.row(i).swap(lm.row(j));
+    VERIFY_IS_APPROX(lm, lp2.toDenseMatrix().template cast<Scalar>());
+
+    RightPermutationType rp2 = rp;
+    i = ei_random<int>(0, cols-1);
+    do j = ei_random<int>(0, cols-1); while(j==i);
+    rp2.applyTranspositionOnTheRight(i, j);
+    rm = rp;
+    rm.col(i).swap(rm.col(j));
+    VERIFY_IS_APPROX(rm, rp2.toDenseMatrix().template cast<Scalar>());
+  }
 }
 
 void test_permutationmatrices()