• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
13 namespace 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 
44 namespace internal {
45 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct transposition_matrix_product_retval;
46 }
47 
48 template<typename Derived>
49 class 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 
derived()58     Derived& derived() { return *static_cast<Derived*>(this); }
derived()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 */
size()81     inline Index size() const { return indices().size(); }
82 
83     /** Direct access to the underlying index vector */
coeff(Index i)84     inline const Index& coeff(Index i) const { return indices().coeff(i); }
85     /** Direct access to the underlying index vector */
coeffRef(Index i)86     inline Index& coeffRef(Index i) { return indices().coeffRef(i); }
87     /** Direct access to the underlying index vector */
operator()88     inline const Index& operator()(Index i) const { return indices()(i); }
89     /** Direct access to the underlying index vector */
operator()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(). */
indices()97     const IndicesType& indices() const { return derived().indices(); }
98     /** \returns a reference to the stored array representing the transpositions. */
indices()99     IndicesType& indices() { return derived().indices(); }
100 
101     /** Resizes to given size. */
resize(int newSize)102     inline void resize(int newSize)
103     {
104       indices().resize(newSize);
105     }
106 
107     /** Sets \c *this to represents an identity transformation */
setIdentity()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 */
inverse()136     inline Transpose<TranspositionsBase> inverse() const
137     { return Transpose<TranspositionsBase>(derived()); }
138 
139     /** \returns the tranpose transformation */
transpose()140     inline Transpose<TranspositionsBase> transpose() const
141     { return Transpose<TranspositionsBase>(derived()); }
142 
143   protected:
144 };
145 
146 namespace internal {
147 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
148 struct traits<Transpositions<SizeAtCompileTime,MaxSizeAtCompileTime,IndexType> >
149 {
150   typedef IndexType Index;
151   typedef Matrix<Index, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
152 };
153 }
154 
155 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
156 class 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 
217 namespace internal {
218 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
219 struct 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 
226 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int PacketAccess>
227 class 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 
274 namespace internal {
275 template<typename _IndicesType>
276 struct traits<TranspositionsWrapper<_IndicesType> >
277 {
278   typedef typename _IndicesType::Scalar Index;
279   typedef _IndicesType IndicesType;
280 };
281 }
282 
283 template<typename _IndicesType>
284 class 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   */
329 template<typename Derived, typename TranspositionsDerived>
330 inline const internal::transposition_matrix_product_retval<TranspositionsDerived, Derived, OnTheRight>
331 operator*(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   */
341 template<typename Derived, typename TranspositionDerived>
342 inline const internal::transposition_matrix_product_retval
343                <TranspositionDerived, Derived, OnTheLeft>
344 operator*(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 
352 namespace internal {
353 
354 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
355 struct traits<transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
356 {
357   typedef typename MatrixType::PlainObject ReturnType;
358 };
359 
360 template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
361 struct 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 
401 template<typename TranspositionsDerived>
402 class 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