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