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);
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