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 template<typename Derived> 16 class 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 118 namespace internal { 119 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 120 struct 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 157 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 158 class 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 219 namespace internal { 220 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess> 221 struct 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 230 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int PacketAccess> 231 class 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 278 namespace internal { 279 template<typename _IndicesType> 280 struct traits<TranspositionsWrapper<_IndicesType> > 281 : traits<PermutationWrapper<_IndicesType> > 282 { 283 typedef TranspositionsStorage StorageKind; 284 }; 285 } 286 287 template<typename _IndicesType> 288 class 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 */ 335 template<typename MatrixDerived, typename TranspositionsDerived> 336 EIGEN_DEVICE_FUNC 337 const Product<MatrixDerived, TranspositionsDerived, AliasFreeProduct> 338 operator*(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 */ 347 template<typename TranspositionsDerived, typename MatrixDerived> 348 EIGEN_DEVICE_FUNC 349 const Product<TranspositionsDerived, MatrixDerived, AliasFreeProduct> 350 operator*(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 359 namespace internal { 360 361 template<typename Derived> 362 struct traits<Transpose<TranspositionsBase<Derived> > > 363 : traits<Derived> 364 {}; 365 366 } // end namespace internal 367 368 template<typename TranspositionsDerived> 369 class 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
