From bfeba41174638c1a19df74436a1572b6f8a6da33 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 4 Jun 2010 23:17:57 +0200 Subject: Add a Transpositions class to ease the representation and manipulation of permutations as a sequence of transpositions. Make LDLT use it. --- Eigen/src/Core/Transpositions.h | 285 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 285 insertions(+) create mode 100644 Eigen/src/Core/Transpositions.h (limited to 'Eigen/src/Core/Transpositions.h') diff --git a/Eigen/src/Core/Transpositions.h b/Eigen/src/Core/Transpositions.h new file mode 100644 index 000000000..39cb24fd7 --- /dev/null +++ b/Eigen/src/Core/Transpositions.h @@ -0,0 +1,285 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_TRANSPOSITIONS_H +#define EIGEN_TRANSPOSITIONS_H + +/** \class Transpositions + * + * \brief Represents a sequence of transpositions (row/column interchange) + * + * \param SizeAtCompileTime the number of transpositions, or Dynamic + * \param MaxSizeAtCompileTime the maximum number of transpositions, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it. + * + * This class represents a permutation transformation as a sequence of \em n transpositions + * \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices. + * Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges + * the rows \c i and \c indices[i] of the matrix \c M. + * A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange. + * + * Compared to the class PermutationMatrix, such a sequence of transpositions is what is + * computed during a decomposition with pivoting, and it is faster when applying the permutation in-place. + * + * To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example: + * \code + * Transpositions tr; + * MatrixXf mat; + * mat = tr * mat; + * \endcode + * In this example, we detect that the matrix appears on both side, and so the transpositions + * are applied in-place without any temporary or extra copy. + * + * \sa class PermutationMatrix + */ +template class Transpositions; +template struct ei_transposition_matrix_product_retval; + +template +class Transpositions +{ + public: + + typedef Matrix IndicesType; + typedef typename IndicesType::Index Index; + + inline Transpositions() {} + + /** 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) + {} + + /** Copies the \a other transpositions into \c *this */ + template + Transpositions& operator=(const Transpositions& other) + { + m_indices = other.indices(); + return *this; + } + + #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) + { + m_indices = other.m_indices; + return *this; + } + #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 const Index& coeff(Index i) const { return m_indices.coeff(i); } + inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); } + inline const Index& operator()(Index i) const { return m_indices(i); } + inline Index& operator()(Index i) { return m_indices(i); } + + /** 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; } + + /** Resizes to given size. */ + inline void resize(int size) + { + m_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; + } + + // FIXME: do we want such methods ? + // might be usefull when the target matrix expression is complex, e.g.: + // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..); + /* + template + void applyForwardToRows(MatrixType& mat) const + { + for(Index k=0 ; k + void applyBackwardToRows(MatrixType& mat) const + { + for(Index k=size()-1 ; k>=0 ; --k) + if(m_indices(k)!=k) + mat.row(k).swap(mat.row(m_indices(k))); + } + */ + + /** \returns the inverse transformation */ + inline Transpose inverse() const + { return *this; } + + /** \returns the tranpose transformation */ + inline Transpose transpose() const + { return *this; } + +#ifndef EIGEN_PARSED_BY_DOXYGEN + template + Transpositions(const Transpose >& other) + : m_indices(other.size()) + { + Index n = size(); + Index j = size-1; + for(Index i=0; i +inline const ei_transposition_matrix_product_retval, Derived, OnTheRight> +operator*(const MatrixBase& matrix, + const Transpositions &transpositions) +{ + return ei_transposition_matrix_product_retval + , Derived, OnTheRight> + (transpositions, matrix.derived()); +} + +/** \returns the \a matrix with the \a transpositions applied to the rows. + */ +template +inline const ei_transposition_matrix_product_retval + , Derived, OnTheLeft> +operator*(const Transpositions &transpositions, + const MatrixBase& matrix) +{ + return ei_transposition_matrix_product_retval + , Derived, OnTheLeft> + (transpositions, matrix.derived()); +} + +template +struct ei_traits > +{ + typedef typename MatrixType::PlainObject ReturnType; +}; + +template +struct ei_transposition_matrix_product_retval + : public ReturnByValue > +{ + typedef typename ei_cleantype::type MatrixTypeNestedCleaned; + typedef typename TranspositionType::Index Index; + + ei_transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix) + : m_transpositions(tr), m_matrix(matrix) + {} + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + template inline void evalTo(Dest& dst) const + { + const int size = m_transpositions.size(); + Index j = 0; + + if(!(ei_is_same_type::ret && ei_extract_data(dst) == ei_extract_data(m_matrix))) + dst = m_matrix; + + for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k +class Transpose > +{ + typedef Transpositions TranspositionType; + typedef typename TranspositionType::IndicesType IndicesType; + public: + + Transpose(const TranspositionType& t) : m_transpositions(t) {} + + inline int size() const { return m_transpositions.size(); } + + /** \returns the \a matrix with the inverse transpositions applied to the columns. + */ + template friend + inline const ei_transposition_matrix_product_retval + operator*(const MatrixBase& matrix, const Transpose& trt) + { + return ei_transposition_matrix_product_retval(trt.m_transpositions, matrix.derived()); + } + + /** \returns the \a matrix with the inverse transpositions applied to the rows. + */ + template + inline const ei_transposition_matrix_product_retval + operator*(const MatrixBase& matrix) const + { + return ei_transposition_matrix_product_retval(m_transpositions, matrix.derived()); + } + + const TranspositionType& nestedTranspositions() const { return m_transpositions; } + + protected: + const TranspositionType& m_transpositions; +}; + +#endif // EIGEN_TRANSPOSITIONS_H -- cgit v1.2.3