1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_TRANSPOSITIONS_H
11#define EIGEN_TRANSPOSITIONS_H
12
13namespace Eigen {
14
15/** \class Transpositions
16  * \ingroup Core_Module
17  *
18  * \brief Represents a sequence of transpositions (row/column interchange)
19  *
20  * \param SizeAtCompileTime the number of transpositions, or Dynamic
21  * \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.
22  *
23  * This class represents a permutation transformation as a sequence of \em n transpositions
24  * \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices.
25  * Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges
26  * the rows \c i and \c indices[i] of the matrix \c M.
27  * A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange.
28  *
29  * Compared to the class PermutationMatrix, such a sequence of transpositions is what is
30  * computed during a decomposition with pivoting, and it is faster when applying the permutation in-place.
31  *
32  * To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example:
33  * \code
34  * Transpositions tr;
35  * MatrixXf mat;
36  * mat = tr * mat;
37  * \endcode
38  * In this example, we detect that the matrix appears on both side, and so the transpositions
39  * are applied in-place without any temporary or extra copy.
40  *
41  * \sa class PermutationMatrix
42  */
43
44namespace internal {
45template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct transposition_matrix_product_retval;
46}
47
48template<typename Derived>
49class TranspositionsBase
50{
51    typedef internal::traits<Derived> Traits;
52
53  public:
54
55    typedef typename Traits::IndicesType IndicesType;
56    typedef typename IndicesType::Scalar Index;
57
58    Derived& derived() { return *static_cast<Derived*>(this); }
59    const Derived& derived() const { return *static_cast<const Derived*>(this); }
60
61    /** Copies the \a other transpositions into \c *this */
62    template<typename OtherDerived>
63    Derived& operator=(const TranspositionsBase<OtherDerived>& other)
64    {
65      indices() = other.indices();
66      return derived();
67    }
68
69    #ifndef EIGEN_PARSED_BY_DOXYGEN
70    /** This is a special case of the templated operator=. Its purpose is to
71      * prevent a default operator= from hiding the templated operator=.
72      */
73    Derived& operator=(const TranspositionsBase& other)
74    {
75      indices() = other.indices();
76      return derived();
77    }
78    #endif
79
80    /** \returns the number of transpositions */
81    inline Index size() const { return indices().size(); }
82
83    /** Direct access to the underlying index vector */
84    inline const Index& coeff(Index i) const { return indices().coeff(i); }
85    /** Direct access to the underlying index vector */
86    inline Index& coeffRef(Index i) { return indices().coeffRef(i); }
87    /** Direct access to the underlying index vector */
88    inline const Index& operator()(Index i) const { return indices()(i); }
89    /** Direct access to the underlying index vector */
90    inline Index& operator()(Index i) { return indices()(i); }
91    /** Direct access to the underlying index vector */
92    inline const Index& operator[](Index i) const { return indices()(i); }
93    /** Direct access to the underlying index vector */
94    inline Index& operator[](Index i) { return indices()(i); }
95
96    /** const version of indices(). */
97    const IndicesType& indices() const { return derived().indices(); }
98    /** \returns a reference to the stored array representing the transpositions. */
99    IndicesType& indices() { return derived().indices(); }
100
101    /** Resizes to given size. */
102    inline void resize(int newSize)
103    {
104      indices().resize(newSize);
105    }
106
107    /** Sets \c *this to represents an identity transformation */
108    void setIdentity()
109    {
110      for(int i = 0; i < indices().size(); ++i)
111        coeffRef(i) = i;
112    }
113
114    // FIXME: do we want such methods ?
115    // might be usefull when the target matrix expression is complex, e.g.:
116    // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
117    /*
118    template<typename MatrixType>
119    void applyForwardToRows(MatrixType& mat) const
120    {
121      for(Index k=0 ; k<size() ; ++k)
122        if(m_indices(k)!=k)
123          mat.row(k).swap(mat.row(m_indices(k)));
124    }
125
126    template<typename MatrixType>
127    void applyBackwardToRows(MatrixType& mat) const
128    {
129      for(Index k=size()-1 ; k>=0 ; --k)
130        if(m_indices(k)!=k)
131          mat.row(k).swap(mat.row(m_indices(k)));
132    }
133    */
134
135    /** \returns the inverse transformation */
136    inline Transpose<TranspositionsBase> inverse() const
137    { return Transpose<TranspositionsBase>(derived()); }
138
139    /** \returns the tranpose transformation */
140    inline Transpose<TranspositionsBase> transpose() const
141    { return Transpose<TranspositionsBase>(derived()); }
142
143  protected:
144};
145
146namespace internal {
147template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
148struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
149{
150  typedef IndexType Index;
151  typedef Matrix<Index, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
152};
153}
154
155template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
156class Transpositions : public TranspositionsBase<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
157{
158    typedef internal::traits<Transpositions> Traits;
159  public:
160
161    typedef TranspositionsBase<Transpositions> Base;
162    typedef typename Traits::IndicesType IndicesType;
163    typedef typename IndicesType::Scalar Index;
164
165    inline Transpositions() {}
166
167    /** Copy constructor. */
168    template<typename OtherDerived>
169    inline Transpositions(const TranspositionsBase<OtherDerived>& other)
170      : m_indices(other.indices()) {}
171
172    #ifndef EIGEN_PARSED_BY_DOXYGEN
173    /** Standard copy constructor. Defined only to prevent a default copy constructor
174      * from hiding the other templated constructor */
175    inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {}
176    #endif
177
178    /** Generic constructor from expression of the transposition indices. */
179    template<typename Other>
180    explicit inline Transpositions(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
181    {}
182
183    /** Copies the \a other transpositions into \c *this */
184    template<typename OtherDerived>
185    Transpositions& operator=(const TranspositionsBase<OtherDerived>& other)
186    {
187      return Base::operator=(other);
188    }
189
190    #ifndef EIGEN_PARSED_BY_DOXYGEN
191    /** This is a special case of the templated operator=. Its purpose is to
192      * prevent a default operator= from hiding the templated operator=.
193      */
194    Transpositions& operator=(const Transpositions& other)
195    {
196      m_indices = other.m_indices;
197      return *this;
198    }
199    #endif
200
201    /** Constructs an uninitialized permutation matrix of given size.
202      */
203    inline Transpositions(Index size) : m_indices(size)
204    {}
205
206    /** const version of indices(). */
207    const IndicesType& indices() const { return m_indices; }
208    /** \returns a reference to the stored array representing the transpositions. */
209    IndicesType& indices() { return m_indices; }
210
211  protected:
212
213    IndicesType m_indices;
214};
215
216
217namespace internal {
218template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
219struct traits<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,_PacketAccess> >
220{
221  typedef IndexType Index;
222  typedef Map<const Matrix<Index,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1>, _PacketAccess> IndicesType;
223};
224}
225
226template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int PacketAccess>
227class Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess>
228 : public TranspositionsBase<Map<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType>,PacketAccess> >
229{
230    typedef internal::traits<Map> Traits;
231  public:
232
233    typedef TranspositionsBase<Map> Base;
234    typedef typename Traits::IndicesType IndicesType;
235    typedef typename IndicesType::Scalar Index;
236
237    inline Map(const Index* indicesPtr)
238      : m_indices(indicesPtr)
239    {}
240
241    inline Map(const Index* indicesPtr, Index size)
242      : m_indices(indicesPtr,size)
243    {}
244
245    /** Copies the \a other transpositions into \c *this */
246    template<typename OtherDerived>
247    Map& operator=(const TranspositionsBase<OtherDerived>& other)
248    {
249      return Base::operator=(other);
250    }
251
252    #ifndef EIGEN_PARSED_BY_DOXYGEN
253    /** This is a special case of the templated operator=. Its purpose is to
254      * prevent a default operator= from hiding the templated operator=.
255      */
256    Map& operator=(const Map& other)
257    {
258      m_indices = other.m_indices;
259      return *this;
260    }
261    #endif
262
263    /** const version of indices(). */
264    const IndicesType& indices() const { return m_indices; }
265
266    /** \returns a reference to the stored array representing the transpositions. */
267    IndicesType& indices() { return m_indices; }
268
269  protected:
270
271    IndicesType m_indices;
272};
273
274namespace internal {
275template<typename _IndicesType>
276struct traits<TranspositionsWrapper<_IndicesType> >
277{
278  typedef typename _IndicesType::Scalar Index;
279  typedef _IndicesType IndicesType;
280};
281}
282
283template<typename _IndicesType>
284class TranspositionsWrapper
285 : public TranspositionsBase<TranspositionsWrapper<_IndicesType> >
286{
287    typedef internal::traits<TranspositionsWrapper> Traits;
288  public:
289
290    typedef TranspositionsBase<TranspositionsWrapper> Base;
291    typedef typename Traits::IndicesType IndicesType;
292    typedef typename IndicesType::Scalar Index;
293
294    inline TranspositionsWrapper(IndicesType& a_indices)
295      : m_indices(a_indices)
296    {}
297
298    /** Copies the \a other transpositions into \c *this */
299    template<typename OtherDerived>
300    TranspositionsWrapper& operator=(const TranspositionsBase<OtherDerived>& other)
301    {
302      return Base::operator=(other);
303    }
304
305    #ifndef EIGEN_PARSED_BY_DOXYGEN
306    /** This is a special case of the templated operator=. Its purpose is to
307      * prevent a default operator= from hiding the templated operator=.
308      */
309    TranspositionsWrapper& operator=(const TranspositionsWrapper& other)
310    {
311      m_indices = other.m_indices;
312      return *this;
313    }
314    #endif
315
316    /** const version of indices(). */
317    const IndicesType& indices() const { return m_indices; }
318
319    /** \returns a reference to the stored array representing the transpositions. */
320    IndicesType& indices() { return m_indices; }
321
322  protected:
323
324    const typename IndicesType::Nested m_indices;
325};
326
327/** \returns the \a matrix with the \a transpositions applied to the columns.
328  */
329template<typename Derived, typename TranspositionsDerived>
330inline const internal::transposition_matrix_product_retval<TranspositionsDerived, Derived, OnTheRight>
331operator*(const MatrixBase<Derived>& matrix,
332          const TranspositionsBase<TranspositionsDerived> &transpositions)
333{
334  return internal::transposition_matrix_product_retval
335           <TranspositionsDerived, Derived, OnTheRight>
336           (transpositions.derived(), matrix.derived());
337}
338
339/** \returns the \a matrix with the \a transpositions applied to the rows.
340  */
341template<typename Derived, typename TranspositionDerived>
342inline const internal::transposition_matrix_product_retval
343               <TranspositionDerived, Derived, OnTheLeft>
344operator*(const TranspositionsBase<TranspositionDerived> &transpositions,
345          const MatrixBase<Derived>& matrix)
346{
347  return internal::transposition_matrix_product_retval
348           <TranspositionDerived, Derived, OnTheLeft>
349           (transpositions.derived(), matrix.derived());
350}
351
352namespace internal {
353
354template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
355struct traits<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
356{
357  typedef typename MatrixType::PlainObject ReturnType;
358};
359
360template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
361struct transposition_matrix_product_retval
362 : public ReturnByValue<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
363{
364    typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
365    typedef typename TranspositionType::Index Index;
366
367    transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix)
368      : m_transpositions(tr), m_matrix(matrix)
369    {}
370
371    inline int rows() const { return m_matrix.rows(); }
372    inline int cols() const { return m_matrix.cols(); }
373
374    template<typename Dest> inline void evalTo(Dest& dst) const
375    {
376      const int size = m_transpositions.size();
377      Index j = 0;
378
379      if(!(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix)))
380        dst = m_matrix;
381
382      for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
383        if((j=m_transpositions.coeff(k))!=k)
384        {
385          if(Side==OnTheLeft)
386            dst.row(k).swap(dst.row(j));
387          else if(Side==OnTheRight)
388            dst.col(k).swap(dst.col(j));
389        }
390    }
391
392  protected:
393    const TranspositionType& m_transpositions;
394    typename MatrixType::Nested m_matrix;
395};
396
397} // end namespace internal
398
399/* Template partial specialization for transposed/inverse transpositions */
400
401template<typename TranspositionsDerived>
402class Transpose<TranspositionsBase<TranspositionsDerived> >
403{
404    typedef TranspositionsDerived TranspositionType;
405    typedef typename TranspositionType::IndicesType IndicesType;
406  public:
407
408    Transpose(const TranspositionType& t) : m_transpositions(t) {}
409
410    inline int size() const { return m_transpositions.size(); }
411
412    /** \returns the \a matrix with the inverse transpositions applied to the columns.
413      */
414    template<typename Derived> friend
415    inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>
416    operator*(const MatrixBase<Derived>& matrix, const Transpose& trt)
417    {
418      return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>(trt.m_transpositions, matrix.derived());
419    }
420
421    /** \returns the \a matrix with the inverse transpositions applied to the rows.
422      */
423    template<typename Derived>
424    inline const internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>
425    operator*(const MatrixBase<Derived>& matrix) const
426    {
427      return internal::transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>(m_transpositions, matrix.derived());
428    }
429
430  protected:
431    const TranspositionType& m_transpositions;
432};
433
434} // end namespace Eigen
435
436#endif // EIGEN_TRANSPOSITIONS_H
437