diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 0a444bd02..40e44c9b3 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -277,6 +277,7 @@ template class MatrixBase static const BasisReturnType UnitW(); const DiagonalWrapper asDiagonal() const; + const PermutationWrapper asPermutation() const; Derived& setIdentity(); Derived& setIdentity(Index rows, Index cols); diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index f53dd8772..06787a66d 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Benoit Jacob -// Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2009-2011 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -26,15 +26,17 @@ #ifndef EIGEN_PERMUTATIONMATRIX_H #define EIGEN_PERMUTATIONMATRIX_H -/** \class PermutationMatrix +template class PermutedImpl; + +/** \class PermutationBase * \ingroup Core_Module * - * \brief Permutation matrix + * \brief Base class for permutations * - * \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. Most of the time, you should not have to specify it. + * \param Derived the derived class * - * This class represents a permutation matrix, internally stored as a vector of integers. + * This class is the base class for all expressions representing a permutation matrix, + * internally stored as a vector of integers. * 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] @@ -44,31 +46,29 @@ * 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. + * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase) + * on either side. * - * \sa class DiagonalMatrix + * \sa class PermutationMatrix, class PermutationWrapper */ namespace internal { -template struct permut_matrix_product_retval; - -template -struct traits > - : traits > -{}; +template +struct permut_matrix_product_retval; +enum PermPermProduct_t {PermPermProduct}; } // end namespace internal -template -class PermutationMatrix : public EigenBase > +template +class PermutationBase : public EigenBase { + typedef internal::traits Traits; + typedef EigenBase Base; public: #ifndef EIGEN_PARSED_BY_DOXYGEN - typedef internal::traits Traits; - typedef Matrix - DenseMatrixType; + typedef typename Traits::IndicesType IndicesType; enum { Flags = Traits::Flags, CoeffReadCost = Traits::CoeffReadCost, @@ -79,9 +79,231 @@ class PermutationMatrix : public EigenBase + DenseMatrixType; + typedef PermutationMatrix + PlainPermutationType; + using Base::derived; #endif - typedef Matrix IndicesType; + + + inline PermutationBase() {} + + /** Copies the other permutation into *this */ + template + Derived& operator=(const PermutationBase& other) + { + indices() = other.indices(); + return derived(); + } + + /** Assignment from the Transpositions \a tr */ + template + Derived& operator=(const TranspositionsBase& tr) + { + setIdentity(tr.size()); + for(Index k=size()-1; k>=0; --k) + applyTranspositionOnTheRight(k,tr.coeff(k)); + return derived(); + } + + #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=. + */ + Derived& operator=(const PermutationBase& other) + { + indices() = other.indices(); + return derived(); + } + #endif + + /** \returns the number of rows */ + inline Index rows() const { return indices().size(); } + + /** \returns the number of columns */ + inline Index cols() const { return indices().size(); } + + /** \returns the size of a side of the respective square matrix, i.e., the number of indices */ + inline Index size() const { return indices().size(); } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template + void evalTo(MatrixBase& other) const + { + other.setZero(); + for (int i=0; i=0 && j>=0 && i=0 && j>=0 && i inverse() const + { return derived(); } + /** \returns the tranpose permutation matrix. + * + * \note \note_try_to_help_rvo + */ + inline Transpose transpose() const + { return derived(); } + + /**** multiplication helpers to hopefully get RVO ****/ + + +#ifndef EIGEN_PARSED_BY_DOXYGEN + protected: + template + void assignTranspose(const PermutationBase& other) + { + for (int i=0; i + void assignProduct(const Lhs& lhs, const Rhs& rhs) + { + eigen_assert(lhs.cols() == rhs.rows()); + for (int i=0; i + inline PlainPermutationType operator*(const PermutationBase& other) const + { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); } + + /** \returns the product of a permutation with another inverse permutation. + * + * \note \note_try_to_help_rvo + */ + template + inline PlainPermutationType operator*(const Transpose >& other) const + { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); } + + /** \returns the product of an inverse permutation with another permutation. + * + * \note \note_try_to_help_rvo + */ + template friend + inline PlainPermutationType operator*(const Transpose >& other, const PermutationBase& perm) + { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); } + + protected: + +}; + +/** \class PermutationMatrix + * \ingroup Core_Module + * + * \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. Most of the time, you should not have to specify it. + * \param IndexType the interger type of the indices + * + * This class represents a permutation matrix, internally stored as a vector of integers. + * + * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix + */ + +namespace internal { +template +struct traits > + : traits > +{ + typedef IndexType Index; + typedef Matrix IndicesType; +}; +} + +template +class PermutationMatrix : public PermutationBase > +{ + typedef PermutationBase Base; + typedef internal::traits Traits; + public: + + #ifndef EIGEN_PARSED_BY_DOXYGEN + typedef typename Traits::IndicesType IndicesType; + #endif inline PermutationMatrix() {} @@ -92,8 +314,8 @@ class PermutationMatrix : public EigenBase - inline PermutationMatrix(const PermutationMatrix& other) + template + inline PermutationMatrix(const PermutationBase& other) : m_indices(other.indices()) {} #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -114,29 +336,26 @@ class PermutationMatrix : public EigenBase - explicit PermutationMatrix(const Transpositions& tr) + template + explicit PermutationMatrix(const TranspositionsBase& tr) : m_indices(tr.size()) { *this = tr; } /** Copies the other permutation into *this */ - template - PermutationMatrix& operator=(const PermutationMatrix& other) + template + PermutationMatrix& operator=(const PermutationBase& other) { m_indices = other.indices(); return *this; } /** Assignment from the Transpositions \a tr */ - template - PermutationMatrix& operator=(const Transpositions& tr) + template + PermutationMatrix& operator=(const TranspositionsBase& tr) { - setIdentity(tr.size()); - for(Index k=size()-1; k>=0; --k) - applyTranspositionOnTheRight(k,tr.coeff(k)); - return *this; + return Base::operator=(tr.derived()); } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -150,182 +369,178 @@ class PermutationMatrix : public EigenBase - void evalTo(MatrixBase& other) const - { - other.setZero(); - for (int i=0; i=0 && j>=0 && i=0 && j>=0 && i inverse() const - { return *this; } - /** \returns the tranpose permutation matrix. - * - * \note \note_try_to_help_rvo - */ - inline Transpose transpose() const - { return *this; } /**** multiplication helpers to hopefully get RVO ****/ #ifndef EIGEN_PARSED_BY_DOXYGEN - template - PermutationMatrix(const Transpose >& other) + template + PermutationMatrix(const Transpose >& other) : m_indices(other.nestedPermutation().size()) { - for (int i=0; i + PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs) + : m_indices(lhs.indices().size()) { - eigen_assert(lhs.cols() == rhs.rows()); - for (int i=0; i - inline PermutationMatrix operator*(const PermutationMatrix& other) const - { return PermutationMatrix(Product, *this, other); } - - /** \returns the product of a permutation with another inverse permutation. - * - * \note \note_try_to_help_rvo - */ - template - inline PermutationMatrix operator*(const Transpose >& other) const - { return PermutationMatrix(Product, *this, other.eval()); } - - /** \returns the product of an inverse permutation with another permutation. - * - * \note \note_try_to_help_rvo - */ - template friend - inline PermutationMatrix operator*(const Transpose >& other, const PermutationMatrix& perm) - { return PermutationMatrix(Product, other.eval(), perm); } - protected: IndicesType m_indices; }; + +namespace internal { +template +struct traits,_PacketAccess> > + : traits > +{ + typedef IndexType Index; + typedef Map, _PacketAccess> IndicesType; +}; +} + +template +class Map,_PacketAccess> + : public PermutationBase,_PacketAccess> > +{ + typedef PermutationBase Base; + typedef internal::traits Traits; + public: + + #ifndef EIGEN_PARSED_BY_DOXYGEN + typedef typename Traits::IndicesType IndicesType; + typedef typename IndicesType::Scalar Index; + #endif + + inline Map(const Index* indices) + : m_indices(indices) + {} + + inline Map(const Index* indices, Index size) + : m_indices(indices,size) + {} + + /** Copies the other permutation into *this */ + template + Map& operator=(const PermutationBase& other) + { return Base::operator=(other.derived()); } + + /** Assignment from the Transpositions \a tr */ + template + Map& operator=(const TranspositionsBase& tr) + { return Base::operator=(tr.derived()); } + + #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=. + */ + Map& operator=(const Map& other) + { + m_indices = other.m_indices; + return *this; + } + #endif + + /** 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; } + + protected: + + IndicesType m_indices; +}; + +/** \class PermutationWrapper + * \ingroup Core_Module + * + * \brief Class to view a vector of integers as a permutation matrix + * + * \param _IndicesType the type of the vector of integer (can be any compatible expression) + * + * This class allows to view any vector expression of integers as a permutation matrix. + * + * \sa class PermutationBase, class PermutationMatrix + */ + +struct PermutationStorage {}; + +template class TranspositionsWrapper; +namespace internal { +template +struct traits > +{ + typedef PermutationStorage StorageKind; + typedef typename _IndicesType::Scalar Scalar; + typedef typename _IndicesType::Scalar Index; + typedef _IndicesType IndicesType; + enum { + RowsAtCompileTime = _IndicesType::SizeAtCompileTime, + ColsAtCompileTime = _IndicesType::SizeAtCompileTime, + MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime, + Flags = 0, + CoeffReadCost = _IndicesType::CoeffReadCost + }; +}; +} + +template +class PermutationWrapper : public PermutationBase > +{ + typedef PermutationBase Base; + typedef internal::traits Traits; + public: + + #ifndef EIGEN_PARSED_BY_DOXYGEN + typedef typename Traits::IndicesType IndicesType; + #endif + + inline PermutationWrapper(const IndicesType& indices) + : m_indices(indices) + {} + + /** const version of indices(). */ + const typename internal::remove_all::type& + indices() const { return m_indices; } + + protected: + + const typename IndicesType::Nested m_indices; +}; + /** \returns the matrix with the permutation applied to the columns. */ -template -inline const internal::permut_matrix_product_retval, Derived, OnTheRight> +template +inline const internal::permut_matrix_product_retval operator*(const MatrixBase& matrix, - const PermutationMatrix &permutation) + const PermutationBase &permutation) { return internal::permut_matrix_product_retval - , Derived, OnTheRight> - (permutation, matrix.derived()); + + (permutation.derived(), matrix.derived()); } /** \returns the matrix with the permutation applied to the rows. */ -template +template inline const internal::permut_matrix_product_retval - , Derived, OnTheLeft> -operator*(const PermutationMatrix &permutation, + +operator*(const PermutationBase &permutation, const MatrixBase& matrix) { return internal::permut_matrix_product_retval - , Derived, OnTheLeft> - (permutation, matrix.derived()); + + (permutation.derived(), matrix.derived()); } namespace internal { @@ -402,25 +617,25 @@ struct permut_matrix_product_retval /* Template partial specialization for transposed/inverse permutations */ -template -struct traits > > - : traits > +template +struct traits > > + : traits {}; } // end namespace internal -template -class Transpose > - : public EigenBase > > +template +class Transpose > + : public EigenBase > > { - typedef PermutationMatrix PermutationType; + typedef Derived PermutationType; typedef typename PermutationType::IndicesType IndicesType; + typedef typename PermutationType::PlainPermutationType PlainPermutationType; public: #ifndef EIGEN_PARSED_BY_DOXYGEN typedef internal::traits Traits; - typedef Matrix - DenseMatrixType; + typedef typename Derived::DenseMatrixType DenseMatrixType; enum { Flags = Traits::Flags, CoeffReadCost = Traits::CoeffReadCost, @@ -448,26 +663,26 @@ class Transpose > #endif /** \return the equivalent permutation matrix */ - PermutationType eval() const { return *this; } + PlainPermutationType eval() const { return *this; } DenseMatrixType toDenseMatrix() const { return *this; } /** \returns the matrix with the inverse permutation applied to the columns. */ - template friend - inline const internal::permut_matrix_product_retval - operator*(const MatrixBase& matrix, const Transpose& trPerm) + template friend + inline const internal::permut_matrix_product_retval + operator*(const MatrixBase& matrix, const Transpose& trPerm) { - return internal::permut_matrix_product_retval(trPerm.m_permutation, matrix.derived()); + return internal::permut_matrix_product_retval(trPerm.m_permutation, matrix.derived()); } /** \returns the matrix with the inverse permutation applied to the rows. */ - template - inline const internal::permut_matrix_product_retval - operator*(const MatrixBase& matrix) const + template + inline const internal::permut_matrix_product_retval + operator*(const MatrixBase& matrix) const { - return internal::permut_matrix_product_retval(m_permutation, matrix.derived()); + return internal::permut_matrix_product_retval(m_permutation, matrix.derived()); } const PermutationType& nestedPermutation() const { return m_permutation; } @@ -476,4 +691,10 @@ class Transpose > const PermutationType& m_permutation; }; +template +const PermutationWrapper MatrixBase::asPermutation() const +{ + return derived(); +} + #endif // EIGEN_PERMUTATIONMATRIX_H diff --git a/Eigen/src/Core/Transpositions.h b/Eigen/src/Core/Transpositions.h index 9735b4eea..928ccb02b 100644 --- a/Eigen/src/Core/Transpositions.h +++ b/Eigen/src/Core/Transpositions.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2010 Gael Guennebaud +// Copyright (C) 2010-2011 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -58,88 +58,72 @@ namespace internal { template struct transposition_matrix_product_retval; } -template -class Transpositions +template +class TranspositionsBase { + typedef internal::traits Traits; + public: - typedef Matrix IndicesType; - typedef typename IndicesType::Index Index; + typedef typename Traits::IndicesType IndicesType; + typedef typename IndicesType::Scalar Index; - inline Transpositions() {} + Derived& derived() { return *static_cast(this); } + const Derived& derived() const { return *static_cast(this); } - /** Copy constructor. */ - template - inline Transpositions(const Transpositions& other) - : m_indices(other.indices()) {} - - #ifndef EIGEN_PARSED_BY_DOXYGEN - /** Standard copy constructor. Defined only to prevent a default copy constructor - * from hiding the other templated constructor */ - inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {} - #endif - - /** Generic constructor from expression of the transposition indices. */ - template - explicit inline Transpositions(const MatrixBase& indices) : m_indices(indices) - {} + inline TranspositionsBase() {} /** Copies the \a other transpositions into \c *this */ - template - Transpositions& operator=(const Transpositions& other) + template + Derived& operator=(const TranspositionsBase& other) { - m_indices = other.indices(); - return *this; + indices() = other.indices(); + return derived(); } #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=. */ - Transpositions& operator=(const Transpositions& other) + Derived& operator=(const TranspositionsBase& other) { - m_indices = other.m_indices; - return *this; + indices() = other.indices(); + return derived(); } #endif - /** Constructs an uninitialized permutation matrix of given size. - */ - inline Transpositions(Index size) : m_indices(size) - {} - /** \returns the number of transpositions */ - inline Index size() const { return m_indices.size(); } + inline Index size() const { return indices().size(); } /** Direct access to the underlying index vector */ - inline const Index& coeff(Index i) const { return m_indices.coeff(i); } + inline const Index& coeff(Index i) const { return indices().coeff(i); } /** Direct access to the underlying index vector */ - inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); } + inline Index& coeffRef(Index i) { return indices().coeffRef(i); } /** Direct access to the underlying index vector */ - inline const Index& operator()(Index i) const { return m_indices(i); } + inline const Index& operator()(Index i) const { return indices()(i); } /** Direct access to the underlying index vector */ - inline Index& operator()(Index i) { return m_indices(i); } + inline Index& operator()(Index i) { return indices()(i); } /** Direct access to the underlying index vector */ - inline const Index& operator[](Index i) const { return m_indices(i); } + inline const Index& operator[](Index i) const { return indices()(i); } /** Direct access to the underlying index vector */ - inline Index& operator[](Index i) { return m_indices(i); } + inline Index& operator[](Index i) { return indices()(i); } /** const version of indices(). */ - const IndicesType& indices() const { return m_indices; } + const IndicesType& indices() const { return derived().indices(); } /** \returns a reference to the stored array representing the transpositions. */ - IndicesType& indices() { return m_indices; } + IndicesType& indices() { return derived().indices(); } /** Resizes to given size. */ inline void resize(int size) { - m_indices.resize(size); + indices().resize(size); } /** Sets \c *this to represents an identity transformation */ void setIdentity() { - for(int i = 0; i < m_indices.size(); ++i) - m_indices.coeffRef(i) = i; + for(int i = 0; i < indices().size(); ++i) + coeffRef(i) = i; } // FIXME: do we want such methods ? @@ -164,53 +148,220 @@ class Transpositions */ /** \returns the inverse transformation */ - inline Transpose inverse() const - { return *this; } + inline Transpose inverse() const + { return Transpose(derived()); } /** \returns the tranpose transformation */ - inline Transpose transpose() const - { return *this; } + inline Transpose transpose() const + { return Transpose(derived()); } -#ifndef EIGEN_PARSED_BY_DOXYGEN - template - Transpositions(const Transpose >& other) - : m_indices(other.size()) + protected: +}; + +namespace internal { +template +struct traits > +{ + typedef IndexType Index; + typedef Matrix IndicesType; +}; +} + +template +class Transpositions : public TranspositionsBase > +{ + typedef internal::traits Traits; + public: + + typedef TranspositionsBase Base; + typedef typename Traits::IndicesType IndicesType; + typedef typename IndicesType::Scalar Index; + + inline Transpositions() {} + + /** Copy constructor. */ + template + inline Transpositions(const TranspositionsBase& other) + : m_indices(other.indices()) {} + + #ifndef EIGEN_PARSED_BY_DOXYGEN + /** Standard copy constructor. Defined only to prevent a default copy constructor + * from hiding the other templated constructor */ + inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {} + #endif + + /** Generic constructor from expression of the transposition indices. */ + template + explicit inline Transpositions(const MatrixBase& indices) : m_indices(indices) + {} + + /** Copies the \a other transpositions into \c *this */ + template + Transpositions& operator=(const TranspositionsBase& other) { - Index n = size(); - Index j = size-1; - for(Index i=0; i +struct traits,_PacketAccess> > +{ + typedef IndexType Index; + typedef Map, _PacketAccess> IndicesType; +}; +} + +template +class Map,PacketAccess> + : public TranspositionsBase,PacketAccess> > +{ + typedef internal::traits Traits; + public: + + typedef TranspositionsBase Base; + typedef typename Traits::IndicesType IndicesType; + typedef typename IndicesType::Scalar Index; + + inline Map(const Index* indices) + : m_indices(indices) + {} + + inline Map(const Index* indices, Index size) + : m_indices(indices,size) + {} + + /** Copies the \a other transpositions into \c *this */ + template + Map& operator=(const TranspositionsBase& other) + { + return Base::operator=(other); + } + + #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=. + */ + Map& operator=(const Map& other) + { + m_indices = other.m_indices; + return *this; + } + #endif + + /** const version of indices(). */ + const IndicesType& indices() const { return m_indices; } + + /** \returns a reference to the stored array representing the transpositions. */ + IndicesType& indices() { return m_indices; } + + protected: + + IndicesType m_indices; +}; + +namespace internal { +template +struct traits > +{ + typedef typename _IndicesType::Scalar Index; + typedef _IndicesType IndicesType; +}; +} + +template +class TranspositionsWrapper + : public TranspositionsBase > +{ + typedef internal::traits Traits; + public: + + typedef TranspositionsBase Base; + typedef typename Traits::IndicesType IndicesType; + typedef typename IndicesType::Scalar Index; + + inline TranspositionsWrapper(IndicesType& indices) + : m_indices(indices) + {} + + /** Copies the \a other transpositions into \c *this */ + template + TranspositionsWrapper& operator=(const TranspositionsBase& other) + { + return Base::operator=(other); + } + + #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=. + */ + TranspositionsWrapper& operator=(const TranspositionsWrapper& other) + { + m_indices = other.m_indices; + return *this; + } + #endif + + /** const version of indices(). */ + const IndicesType& indices() const { return m_indices; } + + /** \returns a reference to the stored array representing the transpositions. */ + IndicesType& indices() { return m_indices; } + + protected: + + const typename IndicesType::Nested m_indices; +}; + /** \returns the \a matrix with the \a transpositions applied to the columns. */ -template -inline const internal::transposition_matrix_product_retval, Derived, OnTheRight> +template +inline const internal::transposition_matrix_product_retval operator*(const MatrixBase& matrix, - const Transpositions &transpositions) + const TranspositionsBase &transpositions) { return internal::transposition_matrix_product_retval - , Derived, OnTheRight> - (transpositions, matrix.derived()); + + (transpositions.derived(), matrix.derived()); } /** \returns the \a matrix with the \a transpositions applied to the rows. */ -template +template inline const internal::transposition_matrix_product_retval - , Derived, OnTheLeft> -operator*(const Transpositions &transpositions, + +operator*(const TranspositionsBase &transpositions, const MatrixBase& matrix) { return internal::transposition_matrix_product_retval - , Derived, OnTheLeft> - (transpositions, matrix.derived()); + + (transpositions.derived(), matrix.derived()); } namespace internal { @@ -262,10 +413,10 @@ struct transposition_matrix_product_retval /* Template partial specialization for transposed/inverse transpositions */ -template -class Transpose > +template +class Transpose > { - typedef Transpositions TranspositionType; + typedef TranspositionsDerived TranspositionType; typedef typename TranspositionType::IndicesType IndicesType; public: @@ -291,8 +442,6 @@ class Transpose > return internal::transposition_matrix_product_retval(m_transpositions, matrix.derived()); } - const TranspositionType& nestedTranspositions() const { return m_transpositions; } - protected: const TranspositionType& m_transpositions; }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 2a65c666d..d4eb2c31c 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -101,8 +101,12 @@ template class DiagonalWrapper; template class DiagonalMatrix; template class DiagonalProduct; template class Diagonal; -template class PermutationMatrix; -template class Transpositions; +template class PermutationMatrix; +template class Transpositions; +template class PermutationBase; +template class TranspositionsBase; +template class PermutationWrapper; +template class TranspositionsWrapper; template::has_write_access ? WriteAccessors : ReadOnlyAccessors diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 8694abae7..f9e645ec1 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -225,7 +225,7 @@ PartialPivLU::PartialPivLU(const MatrixType& matrix) namespace internal { /** \internal This is the blocked version of fullpivlu_unblocked() */ -template +template struct partial_lu_impl { // FIXME add a stride to Map, so that the following mapping becomes easier, @@ -384,13 +384,13 @@ struct partial_lu_impl /** \internal performs the LU decomposition with partial pivoting in-place. */ template -void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename MatrixType::Index& nb_transpositions) +void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions) { eigen_assert(lu.cols() == row_transpositions.size()); eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); partial_lu_impl - + ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); } @@ -406,7 +406,7 @@ PartialPivLU& PartialPivLU::compute(const MatrixType& ma m_rowsTranspositions.resize(size); - Index nb_transpositions; + typename TranspositionType::Index nb_transpositions; internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index b9b3bbbca..d0fa01310 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp @@ -51,8 +51,10 @@ template void permutationmatrices(const MatrixType& m) Options = MatrixType::Options }; typedef PermutationMatrix LeftPermutationType; typedef Matrix LeftPermutationVectorType; + typedef Map MapLeftPerm; typedef PermutationMatrix RightPermutationType; typedef Matrix RightPermutationVectorType; + typedef Map MapRightPerm; Index rows = m.rows(); Index cols = m.cols(); @@ -76,13 +78,20 @@ template void permutationmatrices(const MatrixType& m) VERIFY_IS_APPROX(m_permuted, lm*m_original*rm); VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original); + VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original); + VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original); + VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity()); + VERIFY((lv.asPermutation()*lv.asPermutation().inverse()).toDenseMatrix().isIdentity()); + VERIFY((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv.data(),lv.size()).inverse()).toDenseMatrix().isIdentity()); LeftPermutationVectorType lv2; randomPermutationVector(lv2, rows); LeftPermutationType lp2(lv2); Matrix lm2(lp2); VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast(), lm*lm2); + VERIFY_IS_APPROX((lv.asPermutation()*lv2.asPermutation()).toDenseMatrix().template cast(), lm*lm2); + VERIFY_IS_APPROX((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv2.data(),lv2.size())).toDenseMatrix().template cast(), lm*lm2); LeftPermutationType identityp; identityp.setIdentity(rows); @@ -123,7 +132,7 @@ template void permutationmatrices(const MatrixType& m) rm = rp; rm.col(i).swap(rm.col(j)); VERIFY_IS_APPROX(rm, rp2.toDenseMatrix().template cast()); - } + } } void test_permutationmatrices()