diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index d3e36c73a..6ce357091 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -47,7 +47,6 @@ * * \sa class DiagonalMatrix */ -template class PermutationMatrix; template struct ei_permut_matrix_product_retval; template @@ -78,8 +77,12 @@ class PermutationMatrix : public EigenBase IndicesType; inline PermutationMatrix() - { - } + {} + + /** Constructs an uninitialized permutation matrix of given size. + */ + inline PermutationMatrix(int size) : m_indices(size) + {} /** Copy constructor. */ template @@ -103,6 +106,14 @@ class PermutationMatrix : public EigenBase& indices) : m_indices(indices) {} + /** Convert the Transpositions \a tr to a permutation matrix */ + template + explicit PermutationMatrix(const Transpositions& tr) + : m_indices(tr.size()) + { + *this = tr; + } + /** Copies the other permutation into *this */ template PermutationMatrix& operator=(const PermutationMatrix& other) @@ -111,6 +122,15 @@ class PermutationMatrix : public EigenBase + PermutationMatrix& operator=(const Transpositions& tr) + { + setIdentity(tr.size()); + for(int k=size()-1; k>=0; --k) + applyTranspositionOnTheRight(k,tr.coeff(k)); + } + #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=. @@ -122,11 +142,6 @@ class PermutationMatrix : public EigenBase class Transpositions; template struct ei_transposition_matrix_product_retval; template @@ -108,10 +107,18 @@ class Transpositions /** \returns the number of transpositions */ inline Index size() const { return m_indices.size(); } + /** Direct access to the underlying index vector */ inline const Index& coeff(Index i) const { return m_indices.coeff(i); } + /** Direct access to the underlying index vector */ inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); } + /** Direct access to the underlying index vector */ inline const Index& operator()(Index i) const { return m_indices(i); } + /** Direct access to the underlying index vector */ inline Index& operator()(Index i) { return m_indices(i); } + /** Direct access to the underlying index vector */ + inline const Index& operator[](Index i) const { return m_indices(i); } + /** Direct access to the underlying index vector */ + inline Index& operator[](Index i) { return m_indices(i); } /** const version of indices(). */ const IndicesType& indices() const { return m_indices; } diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 5cf62f4c6..6a9a7941c 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -77,6 +77,8 @@ template class DiagonalWrapper; template class DiagonalMatrix; template class DiagonalProduct; template class Diagonal; +template class PermutationMatrix; +template class Transpositions; template class Stride; template > class Map; diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 0bf1ac3ce..a9172289c 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -73,8 +73,8 @@ template class PartialPivLU typedef typename NumTraits::Real RealScalar; typedef typename ei_traits::StorageKind StorageKind; typedef typename MatrixType::Index Index; - typedef typename ei_plain_col_type::type PermutationVectorType; typedef PermutationMatrix PermutationType; + typedef Transpositions TranspositionType; /** @@ -186,7 +186,7 @@ template class PartialPivLU protected: MatrixType m_lu; PermutationType m_p; - PermutationVectorType m_rowsTranspositions; + TranspositionType m_rowsTranspositions; Index m_det_p; bool m_isInitialized; }; @@ -389,8 +389,8 @@ struct ei_partial_lu_impl /** \internal performs the LU decomposition with partial pivoting in-place. */ -template -void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, typename MatrixType::Index& nb_transpositions) +template +void ei_partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename MatrixType::Index& nb_transpositions) { ei_assert(lu.cols() == row_transpositions.size()); ei_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); @@ -414,9 +414,7 @@ PartialPivLU& PartialPivLU::compute(const MatrixType& ma ei_partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; - m_p.setIdentity(size); - for(Index k = size-1; k >= 0; --k) - m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k)); + m_p = m_rowsTranspositions; m_isInitialized = true; return *this;