Transpositions.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_TRANSPOSITIONS_H
00011 #define EIGEN_TRANSPOSITIONS_H
00012 
00013 namespace Eigen { 
00014 
00015 /** \class Transpositions
00016   * \ingroup Core_Module
00017   *
00018   * \brief Represents a sequence of transpositions (row/column interchange)
00019   *
00020   * \param SizeAtCompileTime the number of transpositions, or Dynamic
00021   * \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.
00022   *
00023   * This class represents a permutation transformation as a sequence of \em n transpositions
00024   * \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices.
00025   * Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges
00026   * the rows \c i and \c indices[i] of the matrix \c M.
00027   * A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange.
00028   *
00029   * Compared to the class PermutationMatrix, such a sequence of transpositions is what is
00030   * computed during a decomposition with pivoting, and it is faster when applying the permutation in-place.
00031   * 
00032   * To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example:
00033   * \code
00034   * Transpositions tr;
00035   * MatrixXf mat;
00036   * mat = tr * mat;
00037   * \endcode
00038   * In this example, we detect that the matrix appears on both side, and so the transpositions
00039   * are applied in-place without any temporary or extra copy.
00040   *
00041   * \sa class PermutationMatrix
00042   */
00043 
00044 namespace internal {
00045 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct transposition_matrix_product_retval;
00046 }
00047 
00048 template<typename Derived>
00049 class TranspositionsBase
00050 {
00051     typedef internal::traits<Derived> Traits;
00052     
00053   public:
00054 
00055     typedef typename Traits::IndicesType IndicesType;
00056     typedef typename IndicesType::Scalar Index;
00057 
00058     Derived& derived() { return *static_cast<Derived*>(this); }
00059     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00060 
00061     /** Copies the \a other transpositions into \c *this */
00062     template<typename OtherDerived>
00063     Derived& operator=(const TranspositionsBase<OtherDerived>& other)
00064     {
00065       indices() = other.indices();
00066       return derived();
00067     }
00068 
00069     #ifndef EIGEN_PARSED_BY_DOXYGEN
00070     /** This is a special case of the templated operator=. Its purpose is to
00071       * prevent a default operator= from hiding the templated operator=.
00072       */
00073     Derived& operator=(const TranspositionsBase& other)
00074     {
00075       indices() = other.indices();
00076       return derived();
00077     }
00078     #endif
00079 
00080     /** \returns the number of transpositions */
00081     inline Index size() const { return indices().size(); }
00082 
00083     /** Direct access to the underlying index vector */
00084     inline const Index& coeff(Index i) const { return indices().coeff(i); }
00085     /** Direct access to the underlying index vector */
00086     inline Index& coeffRef(Index i) { return indices().coeffRef(i); }
00087     /** Direct access to the underlying index vector */
00088     inline const Index& operator()(Index i) const { return indices()(i); }
00089     /** Direct access to the underlying index vector */
00090     inline Index& operator()(Index i) { return indices()(i); }
00091     /** Direct access to the underlying index vector */
00092     inline const Index& operator[](Index i) const { return indices()(i); }
00093     /** Direct access to the underlying index vector */
00094     inline Index& operator[](Index i) { return indices()(i); }
00095 
00096     /** const version of indices(). */
00097     const IndicesType& indices() const { return derived().indices(); }
00098     /** \returns a reference to the stored array representing the transpositions. */
00099     IndicesType& indices() { return derived().indices(); }
00100 
00101     /** Resizes to given size. */
00102     inline void resize(int newSize)
00103     {
00104       indices().resize(newSize);
00105     }
00106 
00107     /** Sets \c *this to represents an identity transformation */
00108     void setIdentity()
00109     {
00110       for(int i = 0; i < indices().size(); ++i)
00111         coeffRef(i) = i;
00112     }
00113 
00114     // FIXME: do we want such methods ?
00115     // might be usefull when the target matrix expression is complex, e.g.:
00116     // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
00117     /*
00118     template<typename MatrixType>
00119     void applyForwardToRows(MatrixType& mat) const
00120     {
00121       for(Index k=0 ; k<size() ; ++k)
00122         if(m_indices(k)!=k)
00123           mat.row(k).swap(mat.row(m_indices(k)));
00124     }
00125 
00126     template<typename MatrixType>
00127     void applyBackwardToRows(MatrixType& mat) const
00128     {
00129       for(Index k=size()-1 ; k>=0 ; --k)
00130         if(m_indices(k)!=k)
00131           mat.row(k).swap(mat.row(m_indices(k)));
00132     }
00133     */
00134 
00135     /** \returns the inverse transformation */
00136     inline Transpose<TranspositionsBase> inverse() const
00137     { return Transpose<TranspositionsBase>(derived()); }
00138 
00139     /** \returns the tranpose transformation */
00140     inline Transpose<TranspositionsBase> transpose() const
00141     { return Transpose<TranspositionsBase>(derived()); }
00142 
00143   protected:
00144 };
00145 
00146 namespace internal {
00147 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00148 struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
00149 {
00150   typedef IndexType Index;
00151   typedef Matrix<Index, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00152 };
00153 }
00154 
00155 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00156 class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
00157 {
00158     typedef internal::traits<Transpositions> Traits;
00159   public:
00160 
00161     typedef TranspositionsBase<Transpositions> Base;
00162     typedef typename Traits::IndicesType IndicesType;
00163     typedef typename IndicesType::Scalar Index;
00164 
00165     inline Transpositions() {}
00166 
00167     /** Copy constructor. */
00168     template<typename OtherDerived>
00169     inline Transpositions(const TranspositionsBase<OtherDerived>& other)
00170       : m_indices(other.indices()) {}
00171 
00172     #ifndef EIGEN_PARSED_BY_DOXYGEN
00173     /** Standard copy constructor. Defined only to prevent a default copy constructor
00174       * from hiding the other templated constructor */
00175     inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {}
00176     #endif
00177 
00178     /** Generic constructor from expression of the transposition indices. */
00179     template<typename Other>
00180     explicit inline Transpositions(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
00181     {}
00182 
00183     /** Copies the \a other transpositions into \c *this */
00184     template<typename OtherDerived>
00185     Transpositions& operator=(const TranspositionsBase<OtherDerived>& other)
00186     {
00187       return Base::operator=(other);
00188     }
00189 
00190     #ifndef EIGEN_PARSED_BY_DOXYGEN
00191     /** This is a special case of the templated operator=. Its purpose is to
00192       * prevent a default operator= from hiding the templated operator=.
00193       */
00194     Transpositions& operator=(const Transpositions& other)
00195     {
00196       m_indices = other.m_indices;
00197       return *this;
00198     }
00199     #endif
00200 
00201     /** Constructs an uninitialized permutation matrix of given size.
00202       */
00203     inline Transpositions(Index size) : m_indices(size)
00204     {}
00205 
00206     /** const version of indices(). */
00207     const IndicesType& indices() const { return m_indices; }
00208     /** \returns a reference to the stored array representing the transpositions. */
00209     IndicesType& indices () { return m_indices; }
00210 
00211   protected:
00212 
00213     IndicesType m_indices;
00214 };
00215 
00216 
00217 namespace internal {
00218 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00219 struct traits<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,_PacketAccess> >
00220 {
00221   typedef IndexType Index;
00222   typedef Map<const Matrix<Index,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1>, _PacketAccess> IndicesType;
00223 };
00224 }
00225 
00226 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int PacketAccess>
00227 class Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess>
00228  : public TranspositionsBase<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess> >
00229 {
00230     typedef internal::traits<Map> Traits;
00231   public:
00232 
00233     typedef TranspositionsBase<Map> Base;
00234     typedef typename Traits::IndicesType IndicesType;
00235     typedef typename IndicesType::Scalar Index;
00236 
00237     inline Map(const Index* indicesPtr)
00238       : m_indices(indicesPtr)
00239     {}
00240 
00241     inline Map(const Index* indicesPtr, Index size)
00242       : m_indices(indicesPtr,size)
00243     {}
00244 
00245     /** Copies the \a other transpositions into \c *this */
00246     template<typename OtherDerived>
00247     Map& operator=(const TranspositionsBase<OtherDerived>& other)
00248     {
00249       return Base::operator=(other);
00250     }
00251 
00252     #ifndef EIGEN_PARSED_BY_DOXYGEN
00253     /** This is a special case of the templated operator=. Its purpose is to
00254       * prevent a default operator= from hiding the templated operator=.
00255       */
00256     Map& operator=(const Map& other)
00257     {
00258       m_indices = other.m_indices;
00259       return *this;
00260     }
00261     #endif
00262 
00263     /** const version of indices(). */
00264     const IndicesType& indices() const { return m_indices; }
00265     
00266     /** \returns a reference to the stored array representing the transpositions. */
00267     IndicesType& indices() { return m_indices; }
00268 
00269   protected:
00270 
00271     IndicesType m_indices;
00272 };
00273 
00274 namespace internal {
00275 template<typename _IndicesType>
00276 struct traits<TranspositionsWrapper<_IndicesType> >
00277 {
00278   typedef typename _IndicesType::Scalar Index;
00279   typedef _IndicesType IndicesType;
00280 };
00281 }
00282 
00283 template<typename _IndicesType>
00284 class TranspositionsWrapper
00285  : public TranspositionsBase<TranspositionsWrapper<_IndicesType> >
00286 {
00287     typedef internal::traits<TranspositionsWrapper> Traits;
00288   public:
00289 
00290     typedef TranspositionsBase<TranspositionsWrapper> Base;
00291     typedef typename Traits::IndicesType IndicesType;
00292     typedef typename IndicesType::Scalar Index;
00293 
00294     inline TranspositionsWrapper(IndicesType& a_indices)
00295       : m_indices(a_indices)
00296     {}
00297 
00298     /** Copies the \a other transpositions into \c *this */
00299     template<typename OtherDerived>
00300     TranspositionsWrapper& operator=(const TranspositionsBase<OtherDerived>& other)
00301     {
00302       return Base::operator=(other);
00303     }
00304 
00305     #ifndef EIGEN_PARSED_BY_DOXYGEN
00306     /** This is a special case of the templated operator=. Its purpose is to
00307       * prevent a default operator= from hiding the templated operator=.
00308       */
00309     TranspositionsWrapper& operator=(const TranspositionsWrapper& other)
00310     {
00311       m_indices = other.m_indices;
00312       return *this;
00313     }
00314     #endif
00315 
00316     /** const version of indices(). */
00317     const IndicesType& indices() const { return m_indices; }
00318 
00319     /** \returns a reference to the stored array representing the transpositions. */
00320     IndicesType& indices() { return m_indices; }
00321 
00322   protected:
00323 
00324     const typename IndicesType::Nested m_indices;
00325 };
00326 
00327 /** \returns the \a matrix with the \a transpositions applied to the columns.
00328   */
00329 template<typename Derived, typename TranspositionsDerived>
00330 inline const internal::transposition_matrix_product_retval<TranspositionsDerived, Derived, OnTheRight>
00331 operator* (const MatrixBase<Derived>& matrix,
00332           const TranspositionsBase<TranspositionsDerived> &transpositions)
00333 {
00334   return internal::transposition_matrix_product_retval
00335            <TranspositionsDerived, Derived, OnTheRight>
00336            (transpositions.derived(), matrix.derived());
00337 }
00338 
00339 /** \returns the \a matrix with the \a transpositions applied to the rows.
00340   */
00341 template<typename Derived, typename TranspositionDerived>
00342 inline const internal::transposition_matrix_product_retval
00343                <TranspositionDerived, Derived, OnTheLeft>
00344 operator* (const TranspositionsBase<TranspositionDerived> &transpositions,
00345           const MatrixBase<Derived>& matrix)
00346 {
00347   return internal::transposition_matrix_product_retval
00348            <TranspositionDerived, Derived, OnTheLeft>
00349            (transpositions.derived(), matrix.derived());
00350 }
00351 
00352 namespace internal {
00353 
00354 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
00355 struct traits<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
00356 {
00357   typedef typename MatrixType::PlainObject ReturnType;
00358 };
00359 
00360 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
00361 struct transposition_matrix_product_retval
00362  : public ReturnByValue<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
00363 {
00364     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00365     typedef typename TranspositionType::Index Index;
00366 
00367     transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix)
00368       : m_transpositions(tr), m_matrix(matrix)
00369     {}
00370 
00371     inline int rows() const { return m_matrix.rows(); }
00372     inline int cols() const { return m_matrix.cols(); }
00373 
00374     template<typename Dest> inline void evalTo(Dest& dst) const
00375     {
00376       const int size = m_transpositions.size();
00377       Index j = 0;
00378 
00379       if(!(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix)))
00380         dst = m_matrix;
00381 
00382       for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
00383         if((j=m_transpositions.coeff(k))!=k)
00384         {
00385           if(Side==OnTheLeft)
00386             dst.row(k).swap(dst.row(j));
00387           else if(Side==OnTheRight)
00388             dst.col(k).swap(dst.col(j));
00389         }
00390     }
00391 
00392   protected:
00393     const TranspositionType& m_transpositions;
00394     typename MatrixType::Nested m_matrix;
00395 };
00396 
00397 } // end namespace internal
00398 
00399 /* Template partial specialization for transposed/inverse transpositions */
00400 
00401 template<typename TranspositionsDerived>
00402 class Transpose<TranspositionsBase<TranspositionsDerived> >
00403 {
00404     typedef TranspositionsDerived TranspositionType;
00405     typedef typename TranspositionType::IndicesType IndicesType;
00406   public:
00407 
00408     Transpose(const TranspositionType& t) : m_transpositions(t) {}
00409 
00410     inline int size() const { return m_transpositions.size(); }
00411 
00412     /** \returns the \a matrix with the inverse transpositions applied to the columns.
00413       */
00414     template<typename Derived> friend
00415     inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>
00416     operator* (const MatrixBase<Derived>& matrix, const Transpose& trt)
00417     {
00418       return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>(trt.m_transpositions, matrix.derived());
00419     }
00420 
00421     /** \returns the \a matrix with the inverse transpositions applied to the rows.
00422       */
00423     template<typename Derived>
00424     inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>
00425     operator* (const MatrixBase<Derived>& matrix) const
00426     {
00427       return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>(m_transpositions, matrix.derived());
00428     }
00429 
00430   protected:
00431     const TranspositionType& m_transpositions;
00432 };
00433 
00434 } // end namespace Eigen
00435 
00436 #endif // EIGEN_TRANSPOSITIONS_H