diff options
| -rw-r--r-- | Eigen/src/Core/PermutationMatrix.h | 69 | ||||
| -rw-r--r-- | Eigen/src/LU/FullPivLU.h | 10 | ||||
| -rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 5 | ||||
| -rw-r--r-- | test/permutationmatrices.cpp | 24 |
4 files changed, 90 insertions, 18 deletions
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; diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index c4affc795..0ef0a371a 100644 --- a/test/permutationmatrices.cpp +++ b/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() |