00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_PERMUTATIONMATRIX_H
00027 #define EIGEN_PERMUTATIONMATRIX_H
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 namespace internal {
00053
00054 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false> struct permut_matrix_product_retval;
00055
00056 template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
00057 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> >
00058 : traits<Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00059 {};
00060
00061 }
00062
00063 template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
00064 class PermutationMatrix : public EigenBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> >
00065 {
00066 public:
00067
00068 #ifndef EIGEN_PARSED_BY_DOXYGEN
00069 typedef internal::traits<PermutationMatrix> Traits;
00070 typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime>
00071 DenseMatrixType;
00072 enum {
00073 Flags = Traits::Flags,
00074 CoeffReadCost = Traits::CoeffReadCost,
00075 RowsAtCompileTime = Traits::RowsAtCompileTime,
00076 ColsAtCompileTime = Traits::ColsAtCompileTime,
00077 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00078 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00079 };
00080 typedef typename Traits::Scalar Scalar;
00081 typedef typename Traits::Index Index;
00082 #endif
00083
00084 typedef Matrix<int, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00085
00086 inline PermutationMatrix()
00087 {}
00088
00089
00090
00091 inline PermutationMatrix(int size) : m_indices(size)
00092 {}
00093
00094
00095 template<int OtherSize, int OtherMaxSize>
00096 inline PermutationMatrix(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
00097 : m_indices(other.indices()) {}
00098
00099 #ifndef EIGEN_PARSED_BY_DOXYGEN
00100
00101
00102 inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00103 #endif
00104
00105
00106
00107
00108
00109
00110
00111
00112 template<typename Other>
00113 explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00114 {}
00115
00116
00117 template<int OtherSize, int OtherMaxSize>
00118 explicit PermutationMatrix(const Transpositions<OtherSize,OtherMaxSize>& tr)
00119 : m_indices(tr.size())
00120 {
00121 *this = tr;
00122 }
00123
00124
00125 template<int OtherSize, int OtherMaxSize>
00126 PermutationMatrix& operator=(const PermutationMatrix<OtherSize, OtherMaxSize>& other)
00127 {
00128 m_indices = other.indices();
00129 return *this;
00130 }
00131
00132
00133 template<int OtherSize, int OtherMaxSize>
00134 PermutationMatrix& operator=(const Transpositions<OtherSize,OtherMaxSize>& tr)
00135 {
00136 setIdentity(tr.size());
00137 for(Index k=size()-1; k>=0; --k)
00138 applyTranspositionOnTheRight(k,tr.coeff(k));
00139 return *this;
00140 }
00141
00142 #ifndef EIGEN_PARSED_BY_DOXYGEN
00143
00144
00145
00146 PermutationMatrix& operator=(const PermutationMatrix& other)
00147 {
00148 m_indices = other.m_indices;
00149 return *this;
00150 }
00151 #endif
00152
00153
00154 inline Index rows() const { return m_indices.size(); }
00155
00156
00157 inline Index cols() const { return m_indices.size(); }
00158
00159
00160 inline Index size() const { return m_indices.size(); }
00161
00162 #ifndef EIGEN_PARSED_BY_DOXYGEN
00163 template<typename DenseDerived>
00164 void evalTo(MatrixBase<DenseDerived>& other) const
00165 {
00166 other.setZero();
00167 for (int i=0; i<rows();++i)
00168 other.coeffRef(m_indices.coeff(i),i) = typename DenseDerived::Scalar(1);
00169 }
00170 #endif
00171
00172
00173
00174
00175
00176 DenseMatrixType toDenseMatrix() const
00177 {
00178 return *this;
00179 }
00180
00181
00182 const IndicesType& indices() const { return m_indices; }
00183
00184 IndicesType& indices() { return m_indices; }
00185
00186
00187
00188 inline void resize(Index size)
00189 {
00190 m_indices.resize(size);
00191 }
00192
00193
00194 void setIdentity()
00195 {
00196 for(Index i = 0; i < m_indices.size(); ++i)
00197 m_indices.coeffRef(i) = i;
00198 }
00199
00200
00201
00202 void setIdentity(Index size)
00203 {
00204 resize(size);
00205 setIdentity();
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 PermutationMatrix& applyTranspositionOnTheLeft(Index i, Index j)
00218 {
00219 eigen_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size());
00220 for(Index k = 0; k < m_indices.size(); ++k)
00221 {
00222 if(m_indices.coeff(k) == i) m_indices.coeffRef(k) = j;
00223 else if(m_indices.coeff(k) == j) m_indices.coeffRef(k) = i;
00224 }
00225 return *this;
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 PermutationMatrix& applyTranspositionOnTheRight(Index i, Index j)
00237 {
00238 eigen_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size());
00239 std::swap(m_indices.coeffRef(i), m_indices.coeffRef(j));
00240 return *this;
00241 }
00242
00243
00244
00245
00246
00247 inline Transpose<PermutationMatrix> inverse() const
00248 { return *this; }
00249
00250
00251
00252
00253 inline Transpose<PermutationMatrix> transpose() const
00254 { return *this; }
00255
00256
00257
00258 #ifndef EIGEN_PARSED_BY_DOXYGEN
00259 template<int OtherSize, int OtherMaxSize>
00260 PermutationMatrix(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other)
00261 : m_indices(other.nestedPermutation().size())
00262 {
00263 for (int i=0; i<rows();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00264 }
00265 protected:
00266 enum Product_t {Product};
00267 PermutationMatrix(Product_t, const PermutationMatrix& lhs, const PermutationMatrix& rhs)
00268 : m_indices(lhs.m_indices.size())
00269 {
00270 eigen_assert(lhs.cols() == rhs.rows());
00271 for (int i=0; i<rows();++i) m_indices.coeffRef(i) = lhs.m_indices.coeff(rhs.m_indices.coeff(i));
00272 }
00273 #endif
00274
00275 public:
00276
00277
00278
00279
00280
00281 template<int OtherSize, int OtherMaxSize>
00282 inline PermutationMatrix operator*(const PermutationMatrix<OtherSize, OtherMaxSize>& other) const
00283 { return PermutationMatrix(Product, *this, other); }
00284
00285
00286
00287
00288
00289 template<int OtherSize, int OtherMaxSize>
00290 inline PermutationMatrix operator*(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other) const
00291 { return PermutationMatrix(Product, *this, other.eval()); }
00292
00293
00294
00295
00296
00297 template<int OtherSize, int OtherMaxSize> friend
00298 inline PermutationMatrix operator*(const Transpose<PermutationMatrix<OtherSize,OtherMaxSize> >& other, const PermutationMatrix& perm)
00299 { return PermutationMatrix(Product, other.eval(), perm); }
00300
00301 protected:
00302
00303 IndicesType m_indices;
00304 };
00305
00306
00307
00308 template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
00309 inline const internal::permut_matrix_product_retval<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
00310 operator*(const MatrixBase<Derived>& matrix,
00311 const PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> &permutation)
00312 {
00313 return internal::permut_matrix_product_retval
00314 <PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
00315 (permutation, matrix.derived());
00316 }
00317
00318
00319
00320 template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
00321 inline const internal::permut_matrix_product_retval
00322 <PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
00323 operator*(const PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> &permutation,
00324 const MatrixBase<Derived>& matrix)
00325 {
00326 return internal::permut_matrix_product_retval
00327 <PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
00328 (permutation, matrix.derived());
00329 }
00330
00331 namespace internal {
00332
00333 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00334 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00335 {
00336 typedef typename MatrixType::PlainObject ReturnType;
00337 };
00338
00339 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00340 struct permut_matrix_product_retval
00341 : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00342 {
00343 typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00344
00345 permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00346 : m_permutation(perm), m_matrix(matrix)
00347 {}
00348
00349 inline int rows() const { return m_matrix.rows(); }
00350 inline int cols() const { return m_matrix.cols(); }
00351
00352 template<typename Dest> inline void evalTo(Dest& dst) const
00353 {
00354 const int n = Side==OnTheLeft ? rows() : cols();
00355
00356 if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
00357 {
00358
00359 Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00360 mask.fill(false);
00361 int r = 0;
00362 while(r < m_permutation.size())
00363 {
00364
00365 while(r<m_permutation.size() && mask[r]) r++;
00366 if(r>=m_permutation.size())
00367 break;
00368
00369 int k0 = r++;
00370 int kPrev = k0;
00371 mask.coeffRef(k0) = true;
00372 for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00373 {
00374 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00375 .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00376 (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00377
00378 mask.coeffRef(k) = true;
00379 kPrev = k;
00380 }
00381 }
00382 }
00383 else
00384 {
00385 for(int i = 0; i < n; ++i)
00386 {
00387 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00388 (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00389
00390 =
00391
00392 Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00393 (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00394 }
00395 }
00396 }
00397
00398 protected:
00399 const PermutationType& m_permutation;
00400 const typename MatrixType::Nested m_matrix;
00401 };
00402
00403
00404
00405 template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
00406 struct traits<Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > >
00407 : traits<Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00408 {};
00409
00410 }
00411
00412 template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
00413 class Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> >
00414 : public EigenBase<Transpose<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> > >
00415 {
00416 typedef PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime> PermutationType;
00417 typedef typename PermutationType::IndicesType IndicesType;
00418 public:
00419
00420 #ifndef EIGEN_PARSED_BY_DOXYGEN
00421 typedef internal::traits<PermutationType> Traits;
00422 typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime>
00423 DenseMatrixType;
00424 enum {
00425 Flags = Traits::Flags,
00426 CoeffReadCost = Traits::CoeffReadCost,
00427 RowsAtCompileTime = Traits::RowsAtCompileTime,
00428 ColsAtCompileTime = Traits::ColsAtCompileTime,
00429 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00430 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00431 };
00432 typedef typename Traits::Scalar Scalar;
00433 #endif
00434
00435 Transpose(const PermutationType& p) : m_permutation(p) {}
00436
00437 inline int rows() const { return m_permutation.rows(); }
00438 inline int cols() const { return m_permutation.cols(); }
00439
00440 #ifndef EIGEN_PARSED_BY_DOXYGEN
00441 template<typename DenseDerived>
00442 void evalTo(MatrixBase<DenseDerived>& other) const
00443 {
00444 other.setZero();
00445 for (int i=0; i<rows();++i)
00446 other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00447 }
00448 #endif
00449
00450
00451 PermutationType eval() const { return *this; }
00452
00453 DenseMatrixType toDenseMatrix() const { return *this; }
00454
00455
00456
00457 template<typename Derived> friend
00458 inline const internal::permut_matrix_product_retval<PermutationType, Derived, OnTheRight, true>
00459 operator*(const MatrixBase<Derived>& matrix, const Transpose& trPerm)
00460 {
00461 return internal::permut_matrix_product_retval<PermutationType, Derived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00462 }
00463
00464
00465
00466 template<typename Derived>
00467 inline const internal::permut_matrix_product_retval<PermutationType, Derived, OnTheLeft, true>
00468 operator*(const MatrixBase<Derived>& matrix) const
00469 {
00470 return internal::permut_matrix_product_retval<PermutationType, Derived, OnTheLeft, true>(m_permutation, matrix.derived());
00471 }
00472
00473 const PermutationType& nestedPermutation() const { return m_permutation; }
00474
00475 protected:
00476 const PermutationType& m_permutation;
00477 };
00478
00479 #endif // EIGEN_PERMUTATIONMATRIX_H