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 (at) 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 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 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
