From 07412b2075a0d2924ae302f4af7ca54c923c9115 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Mon, 16 Nov 2009 21:28:26 -0500
Subject: PermutationMatrix: add setIdentity and transpositions methods LU:
 make use of that

---
 Eigen/src/Core/PermutationMatrix.h | 69 +++++++++++++++++++++++++++++++++-----
 Eigen/src/LU/FullPivLU.h           | 10 +++---
 Eigen/src/LU/PartialPivLU.h        |  5 ++-
 3 files changed, 66 insertions(+), 18 deletions(-)

(limited to 'Eigen')

diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h
index b2acdb50d..147c48734 100644
--- a/Eigen/src/Core/PermutationMatrix.h
+++ b/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 ****/
 
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 62bdfba56..607fe261f 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/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;
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 4eb1162c1..975f79287 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/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;
-- 
cgit v1.2.3