Yoji KURODA / Eigen

Eigne Matrix Class Library

Dependents:   Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more

Eigen Matrix Class Library for mbed.

Finally, you can use Eigen on your mbed!!!

src/Core/PermutationMatrix.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
ykuroda 0:13a5d365ba16 5 // Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 6 //
ykuroda 0:13a5d365ba16 7 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 8 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 10
ykuroda 0:13a5d365ba16 11 #ifndef EIGEN_PERMUTATIONMATRIX_H
ykuroda 0:13a5d365ba16 12 #define EIGEN_PERMUTATIONMATRIX_H
ykuroda 0:13a5d365ba16 13
ykuroda 0:13a5d365ba16 14 namespace Eigen {
ykuroda 0:13a5d365ba16 15
ykuroda 0:13a5d365ba16 16 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
ykuroda 0:13a5d365ba16 17
ykuroda 0:13a5d365ba16 18 /** \class PermutationBase
ykuroda 0:13a5d365ba16 19 * \ingroup Core_Module
ykuroda 0:13a5d365ba16 20 *
ykuroda 0:13a5d365ba16 21 * \brief Base class for permutations
ykuroda 0:13a5d365ba16 22 *
ykuroda 0:13a5d365ba16 23 * \param Derived the derived class
ykuroda 0:13a5d365ba16 24 *
ykuroda 0:13a5d365ba16 25 * This class is the base class for all expressions representing a permutation matrix,
ykuroda 0:13a5d365ba16 26 * internally stored as a vector of integers.
ykuroda 0:13a5d365ba16 27 * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix
ykuroda 0:13a5d365ba16 28 * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have:
ykuroda 0:13a5d365ba16 29 * \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f]
ykuroda 0:13a5d365ba16 30 * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have:
ykuroda 0:13a5d365ba16 31 * \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f]
ykuroda 0:13a5d365ba16 32 *
ykuroda 0:13a5d365ba16 33 * Permutation matrices are square and invertible.
ykuroda 0:13a5d365ba16 34 *
ykuroda 0:13a5d365ba16 35 * Notice that in addition to the member functions and operators listed here, there also are non-member
ykuroda 0:13a5d365ba16 36 * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase)
ykuroda 0:13a5d365ba16 37 * on either side.
ykuroda 0:13a5d365ba16 38 *
ykuroda 0:13a5d365ba16 39 * \sa class PermutationMatrix, class PermutationWrapper
ykuroda 0:13a5d365ba16 40 */
ykuroda 0:13a5d365ba16 41
ykuroda 0:13a5d365ba16 42 namespace internal {
ykuroda 0:13a5d365ba16 43
ykuroda 0:13a5d365ba16 44 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
ykuroda 0:13a5d365ba16 45 struct permut_matrix_product_retval;
ykuroda 0:13a5d365ba16 46 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
ykuroda 0:13a5d365ba16 47 struct permut_sparsematrix_product_retval;
ykuroda 0:13a5d365ba16 48 enum PermPermProduct_t {PermPermProduct};
ykuroda 0:13a5d365ba16 49
ykuroda 0:13a5d365ba16 50 } // end namespace internal
ykuroda 0:13a5d365ba16 51
ykuroda 0:13a5d365ba16 52 template<typename Derived>
ykuroda 0:13a5d365ba16 53 class PermutationBase : public EigenBase<Derived>
ykuroda 0:13a5d365ba16 54 {
ykuroda 0:13a5d365ba16 55 typedef internal::traits<Derived> Traits;
ykuroda 0:13a5d365ba16 56 typedef EigenBase<Derived> Base;
ykuroda 0:13a5d365ba16 57 public:
ykuroda 0:13a5d365ba16 58
ykuroda 0:13a5d365ba16 59 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 60 typedef typename Traits::IndicesType IndicesType;
ykuroda 0:13a5d365ba16 61 enum {
ykuroda 0:13a5d365ba16 62 Flags = Traits::Flags,
ykuroda 0:13a5d365ba16 63 CoeffReadCost = Traits::CoeffReadCost,
ykuroda 0:13a5d365ba16 64 RowsAtCompileTime = Traits::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 65 ColsAtCompileTime = Traits::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 66 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 67 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 68 };
ykuroda 0:13a5d365ba16 69 typedef typename Traits::Scalar Scalar;
ykuroda 0:13a5d365ba16 70 typedef typename Traits::Index Index;
ykuroda 0:13a5d365ba16 71 typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
ykuroda 0:13a5d365ba16 72 DenseMatrixType;
ykuroda 0:13a5d365ba16 73 typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
ykuroda 0:13a5d365ba16 74 PlainPermutationType;
ykuroda 0:13a5d365ba16 75 using Base::derived;
ykuroda 0:13a5d365ba16 76 #endif
ykuroda 0:13a5d365ba16 77
ykuroda 0:13a5d365ba16 78 /** Copies the other permutation into *this */
ykuroda 0:13a5d365ba16 79 template<typename OtherDerived>
ykuroda 0:13a5d365ba16 80 Derived& operator=(const PermutationBase<OtherDerived>& other)
ykuroda 0:13a5d365ba16 81 {
ykuroda 0:13a5d365ba16 82 indices() = other.indices();
ykuroda 0:13a5d365ba16 83 return derived();
ykuroda 0:13a5d365ba16 84 }
ykuroda 0:13a5d365ba16 85
ykuroda 0:13a5d365ba16 86 /** Assignment from the Transpositions \a tr */
ykuroda 0:13a5d365ba16 87 template<typename OtherDerived>
ykuroda 0:13a5d365ba16 88 Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
ykuroda 0:13a5d365ba16 89 {
ykuroda 0:13a5d365ba16 90 setIdentity(tr.size());
ykuroda 0:13a5d365ba16 91 for(Index k=size()-1; k>=0; --k)
ykuroda 0:13a5d365ba16 92 applyTranspositionOnTheRight(k,tr.coeff(k));
ykuroda 0:13a5d365ba16 93 return derived();
ykuroda 0:13a5d365ba16 94 }
ykuroda 0:13a5d365ba16 95
ykuroda 0:13a5d365ba16 96 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 97 /** This is a special case of the templated operator=. Its purpose is to
ykuroda 0:13a5d365ba16 98 * prevent a default operator= from hiding the templated operator=.
ykuroda 0:13a5d365ba16 99 */
ykuroda 0:13a5d365ba16 100 Derived& operator=(const PermutationBase& other)
ykuroda 0:13a5d365ba16 101 {
ykuroda 0:13a5d365ba16 102 indices() = other.indices();
ykuroda 0:13a5d365ba16 103 return derived();
ykuroda 0:13a5d365ba16 104 }
ykuroda 0:13a5d365ba16 105 #endif
ykuroda 0:13a5d365ba16 106
ykuroda 0:13a5d365ba16 107 /** \returns the number of rows */
ykuroda 0:13a5d365ba16 108 inline Index rows() const { return Index(indices().size()); }
ykuroda 0:13a5d365ba16 109
ykuroda 0:13a5d365ba16 110 /** \returns the number of columns */
ykuroda 0:13a5d365ba16 111 inline Index cols() const { return Index(indices().size()); }
ykuroda 0:13a5d365ba16 112
ykuroda 0:13a5d365ba16 113 /** \returns the size of a side of the respective square matrix, i.e., the number of indices */
ykuroda 0:13a5d365ba16 114 inline Index size() const { return Index(indices().size()); }
ykuroda 0:13a5d365ba16 115
ykuroda 0:13a5d365ba16 116 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 117 template<typename DenseDerived>
ykuroda 0:13a5d365ba16 118 void evalTo(MatrixBase<DenseDerived>& other) const
ykuroda 0:13a5d365ba16 119 {
ykuroda 0:13a5d365ba16 120 other.setZero();
ykuroda 0:13a5d365ba16 121 for (int i=0; i<rows();++i)
ykuroda 0:13a5d365ba16 122 other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
ykuroda 0:13a5d365ba16 123 }
ykuroda 0:13a5d365ba16 124 #endif
ykuroda 0:13a5d365ba16 125
ykuroda 0:13a5d365ba16 126 /** \returns a Matrix object initialized from this permutation matrix. Notice that it
ykuroda 0:13a5d365ba16 127 * is inefficient to return this Matrix object by value. For efficiency, favor using
ykuroda 0:13a5d365ba16 128 * the Matrix constructor taking EigenBase objects.
ykuroda 0:13a5d365ba16 129 */
ykuroda 0:13a5d365ba16 130 DenseMatrixType toDenseMatrix() const
ykuroda 0:13a5d365ba16 131 {
ykuroda 0:13a5d365ba16 132 return derived();
ykuroda 0:13a5d365ba16 133 }
ykuroda 0:13a5d365ba16 134
ykuroda 0:13a5d365ba16 135 /** const version of indices(). */
ykuroda 0:13a5d365ba16 136 const IndicesType& indices() const { return derived().indices(); }
ykuroda 0:13a5d365ba16 137 /** \returns a reference to the stored array representing the permutation. */
ykuroda 0:13a5d365ba16 138 IndicesType& indices() { return derived().indices(); }
ykuroda 0:13a5d365ba16 139
ykuroda 0:13a5d365ba16 140 /** Resizes to given size.
ykuroda 0:13a5d365ba16 141 */
ykuroda 0:13a5d365ba16 142 inline void resize(Index newSize)
ykuroda 0:13a5d365ba16 143 {
ykuroda 0:13a5d365ba16 144 indices().resize(newSize);
ykuroda 0:13a5d365ba16 145 }
ykuroda 0:13a5d365ba16 146
ykuroda 0:13a5d365ba16 147 /** Sets *this to be the identity permutation matrix */
ykuroda 0:13a5d365ba16 148 void setIdentity()
ykuroda 0:13a5d365ba16 149 {
ykuroda 0:13a5d365ba16 150 for(Index i = 0; i < size(); ++i)
ykuroda 0:13a5d365ba16 151 indices().coeffRef(i) = i;
ykuroda 0:13a5d365ba16 152 }
ykuroda 0:13a5d365ba16 153
ykuroda 0:13a5d365ba16 154 /** Sets *this to be the identity permutation matrix of given size.
ykuroda 0:13a5d365ba16 155 */
ykuroda 0:13a5d365ba16 156 void setIdentity(Index newSize)
ykuroda 0:13a5d365ba16 157 {
ykuroda 0:13a5d365ba16 158 resize(newSize);
ykuroda 0:13a5d365ba16 159 setIdentity();
ykuroda 0:13a5d365ba16 160 }
ykuroda 0:13a5d365ba16 161
ykuroda 0:13a5d365ba16 162 /** Multiplies *this by the transposition \f$(ij)\f$ on the left.
ykuroda 0:13a5d365ba16 163 *
ykuroda 0:13a5d365ba16 164 * \returns a reference to *this.
ykuroda 0:13a5d365ba16 165 *
ykuroda 0:13a5d365ba16 166 * \warning This is much slower than applyTranspositionOnTheRight(int,int):
ykuroda 0:13a5d365ba16 167 * this has linear complexity and requires a lot of branching.
ykuroda 0:13a5d365ba16 168 *
ykuroda 0:13a5d365ba16 169 * \sa applyTranspositionOnTheRight(int,int)
ykuroda 0:13a5d365ba16 170 */
ykuroda 0:13a5d365ba16 171 Derived& applyTranspositionOnTheLeft(Index i, Index j)
ykuroda 0:13a5d365ba16 172 {
ykuroda 0:13a5d365ba16 173 eigen_assert(i>=0 && j>=0 && i<size() && j<size());
ykuroda 0:13a5d365ba16 174 for(Index k = 0; k < size(); ++k)
ykuroda 0:13a5d365ba16 175 {
ykuroda 0:13a5d365ba16 176 if(indices().coeff(k) == i) indices().coeffRef(k) = j;
ykuroda 0:13a5d365ba16 177 else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
ykuroda 0:13a5d365ba16 178 }
ykuroda 0:13a5d365ba16 179 return derived();
ykuroda 0:13a5d365ba16 180 }
ykuroda 0:13a5d365ba16 181
ykuroda 0:13a5d365ba16 182 /** Multiplies *this by the transposition \f$(ij)\f$ on the right.
ykuroda 0:13a5d365ba16 183 *
ykuroda 0:13a5d365ba16 184 * \returns a reference to *this.
ykuroda 0:13a5d365ba16 185 *
ykuroda 0:13a5d365ba16 186 * This is a fast operation, it only consists in swapping two indices.
ykuroda 0:13a5d365ba16 187 *
ykuroda 0:13a5d365ba16 188 * \sa applyTranspositionOnTheLeft(int,int)
ykuroda 0:13a5d365ba16 189 */
ykuroda 0:13a5d365ba16 190 Derived& applyTranspositionOnTheRight(Index i, Index j)
ykuroda 0:13a5d365ba16 191 {
ykuroda 0:13a5d365ba16 192 eigen_assert(i>=0 && j>=0 && i<size() && j<size());
ykuroda 0:13a5d365ba16 193 std::swap(indices().coeffRef(i), indices().coeffRef(j));
ykuroda 0:13a5d365ba16 194 return derived();
ykuroda 0:13a5d365ba16 195 }
ykuroda 0:13a5d365ba16 196
ykuroda 0:13a5d365ba16 197 /** \returns the inverse permutation matrix.
ykuroda 0:13a5d365ba16 198 *
ykuroda 0:13a5d365ba16 199 * \note \note_try_to_help_rvo
ykuroda 0:13a5d365ba16 200 */
ykuroda 0:13a5d365ba16 201 inline Transpose<PermutationBase> inverse() const
ykuroda 0:13a5d365ba16 202 { return derived(); }
ykuroda 0:13a5d365ba16 203 /** \returns the tranpose permutation matrix.
ykuroda 0:13a5d365ba16 204 *
ykuroda 0:13a5d365ba16 205 * \note \note_try_to_help_rvo
ykuroda 0:13a5d365ba16 206 */
ykuroda 0:13a5d365ba16 207 inline Transpose<PermutationBase> transpose() const
ykuroda 0:13a5d365ba16 208 { return derived(); }
ykuroda 0:13a5d365ba16 209
ykuroda 0:13a5d365ba16 210 /**** multiplication helpers to hopefully get RVO ****/
ykuroda 0:13a5d365ba16 211
ykuroda 0:13a5d365ba16 212
ykuroda 0:13a5d365ba16 213 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 214 protected:
ykuroda 0:13a5d365ba16 215 template<typename OtherDerived>
ykuroda 0:13a5d365ba16 216 void assignTranspose(const PermutationBase<OtherDerived>& other)
ykuroda 0:13a5d365ba16 217 {
ykuroda 0:13a5d365ba16 218 for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
ykuroda 0:13a5d365ba16 219 }
ykuroda 0:13a5d365ba16 220 template<typename Lhs,typename Rhs>
ykuroda 0:13a5d365ba16 221 void assignProduct(const Lhs& lhs, const Rhs& rhs)
ykuroda 0:13a5d365ba16 222 {
ykuroda 0:13a5d365ba16 223 eigen_assert(lhs.cols() == rhs.rows());
ykuroda 0:13a5d365ba16 224 for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
ykuroda 0:13a5d365ba16 225 }
ykuroda 0:13a5d365ba16 226 #endif
ykuroda 0:13a5d365ba16 227
ykuroda 0:13a5d365ba16 228 public:
ykuroda 0:13a5d365ba16 229
ykuroda 0:13a5d365ba16 230 /** \returns the product permutation matrix.
ykuroda 0:13a5d365ba16 231 *
ykuroda 0:13a5d365ba16 232 * \note \note_try_to_help_rvo
ykuroda 0:13a5d365ba16 233 */
ykuroda 0:13a5d365ba16 234 template<typename Other>
ykuroda 0:13a5d365ba16 235 inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
ykuroda 0:13a5d365ba16 236 { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
ykuroda 0:13a5d365ba16 237
ykuroda 0:13a5d365ba16 238 /** \returns the product of a permutation with another inverse permutation.
ykuroda 0:13a5d365ba16 239 *
ykuroda 0:13a5d365ba16 240 * \note \note_try_to_help_rvo
ykuroda 0:13a5d365ba16 241 */
ykuroda 0:13a5d365ba16 242 template<typename Other>
ykuroda 0:13a5d365ba16 243 inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
ykuroda 0:13a5d365ba16 244 { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
ykuroda 0:13a5d365ba16 245
ykuroda 0:13a5d365ba16 246 /** \returns the product of an inverse permutation with another permutation.
ykuroda 0:13a5d365ba16 247 *
ykuroda 0:13a5d365ba16 248 * \note \note_try_to_help_rvo
ykuroda 0:13a5d365ba16 249 */
ykuroda 0:13a5d365ba16 250 template<typename Other> friend
ykuroda 0:13a5d365ba16 251 inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
ykuroda 0:13a5d365ba16 252 { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
ykuroda 0:13a5d365ba16 253
ykuroda 0:13a5d365ba16 254 /** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation.
ykuroda 0:13a5d365ba16 255 *
ykuroda 0:13a5d365ba16 256 * This function is O(\c n) procedure allocating a buffer of \c n booleans.
ykuroda 0:13a5d365ba16 257 */
ykuroda 0:13a5d365ba16 258 Index determinant() const
ykuroda 0:13a5d365ba16 259 {
ykuroda 0:13a5d365ba16 260 Index res = 1;
ykuroda 0:13a5d365ba16 261 Index n = size();
ykuroda 0:13a5d365ba16 262 Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n);
ykuroda 0:13a5d365ba16 263 mask.fill(false);
ykuroda 0:13a5d365ba16 264 Index r = 0;
ykuroda 0:13a5d365ba16 265 while(r < n)
ykuroda 0:13a5d365ba16 266 {
ykuroda 0:13a5d365ba16 267 // search for the next seed
ykuroda 0:13a5d365ba16 268 while(r<n && mask[r]) r++;
ykuroda 0:13a5d365ba16 269 if(r>=n)
ykuroda 0:13a5d365ba16 270 break;
ykuroda 0:13a5d365ba16 271 // we got one, let's follow it until we are back to the seed
ykuroda 0:13a5d365ba16 272 Index k0 = r++;
ykuroda 0:13a5d365ba16 273 mask.coeffRef(k0) = true;
ykuroda 0:13a5d365ba16 274 for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
ykuroda 0:13a5d365ba16 275 {
ykuroda 0:13a5d365ba16 276 mask.coeffRef(k) = true;
ykuroda 0:13a5d365ba16 277 res = -res;
ykuroda 0:13a5d365ba16 278 }
ykuroda 0:13a5d365ba16 279 }
ykuroda 0:13a5d365ba16 280 return res;
ykuroda 0:13a5d365ba16 281 }
ykuroda 0:13a5d365ba16 282
ykuroda 0:13a5d365ba16 283 protected:
ykuroda 0:13a5d365ba16 284
ykuroda 0:13a5d365ba16 285 };
ykuroda 0:13a5d365ba16 286
ykuroda 0:13a5d365ba16 287 /** \class PermutationMatrix
ykuroda 0:13a5d365ba16 288 * \ingroup Core_Module
ykuroda 0:13a5d365ba16 289 *
ykuroda 0:13a5d365ba16 290 * \brief Permutation matrix
ykuroda 0:13a5d365ba16 291 *
ykuroda 0:13a5d365ba16 292 * \param SizeAtCompileTime the number of rows/cols, or Dynamic
ykuroda 0:13a5d365ba16 293 * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
ykuroda 0:13a5d365ba16 294 * \param IndexType the interger type of the indices
ykuroda 0:13a5d365ba16 295 *
ykuroda 0:13a5d365ba16 296 * This class represents a permutation matrix, internally stored as a vector of integers.
ykuroda 0:13a5d365ba16 297 *
ykuroda 0:13a5d365ba16 298 * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix
ykuroda 0:13a5d365ba16 299 */
ykuroda 0:13a5d365ba16 300
ykuroda 0:13a5d365ba16 301 namespace internal {
ykuroda 0:13a5d365ba16 302 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
ykuroda 0:13a5d365ba16 303 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
ykuroda 0:13a5d365ba16 304 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
ykuroda 0:13a5d365ba16 305 {
ykuroda 0:13a5d365ba16 306 typedef IndexType Index;
ykuroda 0:13a5d365ba16 307 typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
ykuroda 0:13a5d365ba16 308 };
ykuroda 0:13a5d365ba16 309 }
ykuroda 0:13a5d365ba16 310
ykuroda 0:13a5d365ba16 311 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
ykuroda 0:13a5d365ba16 312 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
ykuroda 0:13a5d365ba16 313 {
ykuroda 0:13a5d365ba16 314 typedef PermutationBase<PermutationMatrix> Base;
ykuroda 0:13a5d365ba16 315 typedef internal::traits<PermutationMatrix> Traits;
ykuroda 0:13a5d365ba16 316 public:
ykuroda 0:13a5d365ba16 317
ykuroda 0:13a5d365ba16 318 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 319 typedef typename Traits::IndicesType IndicesType;
ykuroda 0:13a5d365ba16 320 #endif
ykuroda 0:13a5d365ba16 321
ykuroda 0:13a5d365ba16 322 inline PermutationMatrix()
ykuroda 0:13a5d365ba16 323 {}
ykuroda 0:13a5d365ba16 324
ykuroda 0:13a5d365ba16 325 /** Constructs an uninitialized permutation matrix of given size.
ykuroda 0:13a5d365ba16 326 */
ykuroda 0:13a5d365ba16 327 inline PermutationMatrix(int size) : m_indices(size)
ykuroda 0:13a5d365ba16 328 {}
ykuroda 0:13a5d365ba16 329
ykuroda 0:13a5d365ba16 330 /** Copy constructor. */
ykuroda 0:13a5d365ba16 331 template<typename OtherDerived>
ykuroda 0:13a5d365ba16 332 inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
ykuroda 0:13a5d365ba16 333 : m_indices(other.indices()) {}
ykuroda 0:13a5d365ba16 334
ykuroda 0:13a5d365ba16 335 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 336 /** Standard copy constructor. Defined only to prevent a default copy constructor
ykuroda 0:13a5d365ba16 337 * from hiding the other templated constructor */
ykuroda 0:13a5d365ba16 338 inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
ykuroda 0:13a5d365ba16 339 #endif
ykuroda 0:13a5d365ba16 340
ykuroda 0:13a5d365ba16 341 /** Generic constructor from expression of the indices. The indices
ykuroda 0:13a5d365ba16 342 * array has the meaning that the permutations sends each integer i to indices[i].
ykuroda 0:13a5d365ba16 343 *
ykuroda 0:13a5d365ba16 344 * \warning It is your responsibility to check that the indices array that you passes actually
ykuroda 0:13a5d365ba16 345 * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the
ykuroda 0:13a5d365ba16 346 * array's size.
ykuroda 0:13a5d365ba16 347 */
ykuroda 0:13a5d365ba16 348 template<typename Other>
ykuroda 0:13a5d365ba16 349 explicit inline PermutationMatrix(const MatrixBase<Other>& a_indices) : m_indices(a_indices)
ykuroda 0:13a5d365ba16 350 {}
ykuroda 0:13a5d365ba16 351
ykuroda 0:13a5d365ba16 352 /** Convert the Transpositions \a tr to a permutation matrix */
ykuroda 0:13a5d365ba16 353 template<typename Other>
ykuroda 0:13a5d365ba16 354 explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
ykuroda 0:13a5d365ba16 355 : m_indices(tr.size())
ykuroda 0:13a5d365ba16 356 {
ykuroda 0:13a5d365ba16 357 *this = tr;
ykuroda 0:13a5d365ba16 358 }
ykuroda 0:13a5d365ba16 359
ykuroda 0:13a5d365ba16 360 /** Copies the other permutation into *this */
ykuroda 0:13a5d365ba16 361 template<typename Other>
ykuroda 0:13a5d365ba16 362 PermutationMatrix& operator=(const PermutationBase<Other>& other)
ykuroda 0:13a5d365ba16 363 {
ykuroda 0:13a5d365ba16 364 m_indices = other.indices();
ykuroda 0:13a5d365ba16 365 return *this;
ykuroda 0:13a5d365ba16 366 }
ykuroda 0:13a5d365ba16 367
ykuroda 0:13a5d365ba16 368 /** Assignment from the Transpositions \a tr */
ykuroda 0:13a5d365ba16 369 template<typename Other>
ykuroda 0:13a5d365ba16 370 PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
ykuroda 0:13a5d365ba16 371 {
ykuroda 0:13a5d365ba16 372 return Base::operator=(tr.derived());
ykuroda 0:13a5d365ba16 373 }
ykuroda 0:13a5d365ba16 374
ykuroda 0:13a5d365ba16 375 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 376 /** This is a special case of the templated operator=. Its purpose is to
ykuroda 0:13a5d365ba16 377 * prevent a default operator= from hiding the templated operator=.
ykuroda 0:13a5d365ba16 378 */
ykuroda 0:13a5d365ba16 379 PermutationMatrix& operator=(const PermutationMatrix& other)
ykuroda 0:13a5d365ba16 380 {
ykuroda 0:13a5d365ba16 381 m_indices = other.m_indices;
ykuroda 0:13a5d365ba16 382 return *this;
ykuroda 0:13a5d365ba16 383 }
ykuroda 0:13a5d365ba16 384 #endif
ykuroda 0:13a5d365ba16 385
ykuroda 0:13a5d365ba16 386 /** const version of indices(). */
ykuroda 0:13a5d365ba16 387 const IndicesType& indices() const { return m_indices; }
ykuroda 0:13a5d365ba16 388 /** \returns a reference to the stored array representing the permutation. */
ykuroda 0:13a5d365ba16 389 IndicesType& indices() { return m_indices; }
ykuroda 0:13a5d365ba16 390
ykuroda 0:13a5d365ba16 391
ykuroda 0:13a5d365ba16 392 /**** multiplication helpers to hopefully get RVO ****/
ykuroda 0:13a5d365ba16 393
ykuroda 0:13a5d365ba16 394 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 395 template<typename Other>
ykuroda 0:13a5d365ba16 396 PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
ykuroda 0:13a5d365ba16 397 : m_indices(other.nestedPermutation().size())
ykuroda 0:13a5d365ba16 398 {
ykuroda 0:13a5d365ba16 399 for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
ykuroda 0:13a5d365ba16 400 }
ykuroda 0:13a5d365ba16 401 template<typename Lhs,typename Rhs>
ykuroda 0:13a5d365ba16 402 PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
ykuroda 0:13a5d365ba16 403 : m_indices(lhs.indices().size())
ykuroda 0:13a5d365ba16 404 {
ykuroda 0:13a5d365ba16 405 Base::assignProduct(lhs,rhs);
ykuroda 0:13a5d365ba16 406 }
ykuroda 0:13a5d365ba16 407 #endif
ykuroda 0:13a5d365ba16 408
ykuroda 0:13a5d365ba16 409 protected:
ykuroda 0:13a5d365ba16 410
ykuroda 0:13a5d365ba16 411 IndicesType m_indices;
ykuroda 0:13a5d365ba16 412 };
ykuroda 0:13a5d365ba16 413
ykuroda 0:13a5d365ba16 414
ykuroda 0:13a5d365ba16 415 namespace internal {
ykuroda 0:13a5d365ba16 416 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
ykuroda 0:13a5d365ba16 417 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
ykuroda 0:13a5d365ba16 418 : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
ykuroda 0:13a5d365ba16 419 {
ykuroda 0:13a5d365ba16 420 typedef IndexType Index;
ykuroda 0:13a5d365ba16 421 typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
ykuroda 0:13a5d365ba16 422 };
ykuroda 0:13a5d365ba16 423 }
ykuroda 0:13a5d365ba16 424
ykuroda 0:13a5d365ba16 425 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
ykuroda 0:13a5d365ba16 426 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
ykuroda 0:13a5d365ba16 427 : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
ykuroda 0:13a5d365ba16 428 {
ykuroda 0:13a5d365ba16 429 typedef PermutationBase<Map> Base;
ykuroda 0:13a5d365ba16 430 typedef internal::traits<Map> Traits;
ykuroda 0:13a5d365ba16 431 public:
ykuroda 0:13a5d365ba16 432
ykuroda 0:13a5d365ba16 433 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 434 typedef typename Traits::IndicesType IndicesType;
ykuroda 0:13a5d365ba16 435 typedef typename IndicesType::Scalar Index;
ykuroda 0:13a5d365ba16 436 #endif
ykuroda 0:13a5d365ba16 437
ykuroda 0:13a5d365ba16 438 inline Map(const Index* indicesPtr)
ykuroda 0:13a5d365ba16 439 : m_indices(indicesPtr)
ykuroda 0:13a5d365ba16 440 {}
ykuroda 0:13a5d365ba16 441
ykuroda 0:13a5d365ba16 442 inline Map(const Index* indicesPtr, Index size)
ykuroda 0:13a5d365ba16 443 : m_indices(indicesPtr,size)
ykuroda 0:13a5d365ba16 444 {}
ykuroda 0:13a5d365ba16 445
ykuroda 0:13a5d365ba16 446 /** Copies the other permutation into *this */
ykuroda 0:13a5d365ba16 447 template<typename Other>
ykuroda 0:13a5d365ba16 448 Map& operator=(const PermutationBase<Other>& other)
ykuroda 0:13a5d365ba16 449 { return Base::operator=(other.derived()); }
ykuroda 0:13a5d365ba16 450
ykuroda 0:13a5d365ba16 451 /** Assignment from the Transpositions \a tr */
ykuroda 0:13a5d365ba16 452 template<typename Other>
ykuroda 0:13a5d365ba16 453 Map& operator=(const TranspositionsBase<Other>& tr)
ykuroda 0:13a5d365ba16 454 { return Base::operator=(tr.derived()); }
ykuroda 0:13a5d365ba16 455
ykuroda 0:13a5d365ba16 456 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 457 /** This is a special case of the templated operator=. Its purpose is to
ykuroda 0:13a5d365ba16 458 * prevent a default operator= from hiding the templated operator=.
ykuroda 0:13a5d365ba16 459 */
ykuroda 0:13a5d365ba16 460 Map& operator=(const Map& other)
ykuroda 0:13a5d365ba16 461 {
ykuroda 0:13a5d365ba16 462 m_indices = other.m_indices;
ykuroda 0:13a5d365ba16 463 return *this;
ykuroda 0:13a5d365ba16 464 }
ykuroda 0:13a5d365ba16 465 #endif
ykuroda 0:13a5d365ba16 466
ykuroda 0:13a5d365ba16 467 /** const version of indices(). */
ykuroda 0:13a5d365ba16 468 const IndicesType& indices() const { return m_indices; }
ykuroda 0:13a5d365ba16 469 /** \returns a reference to the stored array representing the permutation. */
ykuroda 0:13a5d365ba16 470 IndicesType& indices() { return m_indices; }
ykuroda 0:13a5d365ba16 471
ykuroda 0:13a5d365ba16 472 protected:
ykuroda 0:13a5d365ba16 473
ykuroda 0:13a5d365ba16 474 IndicesType m_indices;
ykuroda 0:13a5d365ba16 475 };
ykuroda 0:13a5d365ba16 476
ykuroda 0:13a5d365ba16 477 /** \class PermutationWrapper
ykuroda 0:13a5d365ba16 478 * \ingroup Core_Module
ykuroda 0:13a5d365ba16 479 *
ykuroda 0:13a5d365ba16 480 * \brief Class to view a vector of integers as a permutation matrix
ykuroda 0:13a5d365ba16 481 *
ykuroda 0:13a5d365ba16 482 * \param _IndicesType the type of the vector of integer (can be any compatible expression)
ykuroda 0:13a5d365ba16 483 *
ykuroda 0:13a5d365ba16 484 * This class allows to view any vector expression of integers as a permutation matrix.
ykuroda 0:13a5d365ba16 485 *
ykuroda 0:13a5d365ba16 486 * \sa class PermutationBase, class PermutationMatrix
ykuroda 0:13a5d365ba16 487 */
ykuroda 0:13a5d365ba16 488
ykuroda 0:13a5d365ba16 489 struct PermutationStorage {};
ykuroda 0:13a5d365ba16 490
ykuroda 0:13a5d365ba16 491 template<typename _IndicesType> class TranspositionsWrapper;
ykuroda 0:13a5d365ba16 492 namespace internal {
ykuroda 0:13a5d365ba16 493 template<typename _IndicesType>
ykuroda 0:13a5d365ba16 494 struct traits<PermutationWrapper<_IndicesType> >
ykuroda 0:13a5d365ba16 495 {
ykuroda 0:13a5d365ba16 496 typedef PermutationStorage StorageKind;
ykuroda 0:13a5d365ba16 497 typedef typename _IndicesType::Scalar Scalar;
ykuroda 0:13a5d365ba16 498 typedef typename _IndicesType::Scalar Index;
ykuroda 0:13a5d365ba16 499 typedef _IndicesType IndicesType;
ykuroda 0:13a5d365ba16 500 enum {
ykuroda 0:13a5d365ba16 501 RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
ykuroda 0:13a5d365ba16 502 ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
ykuroda 0:13a5d365ba16 503 MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 504 MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
ykuroda 0:13a5d365ba16 505 Flags = 0,
ykuroda 0:13a5d365ba16 506 CoeffReadCost = _IndicesType::CoeffReadCost
ykuroda 0:13a5d365ba16 507 };
ykuroda 0:13a5d365ba16 508 };
ykuroda 0:13a5d365ba16 509 }
ykuroda 0:13a5d365ba16 510
ykuroda 0:13a5d365ba16 511 template<typename _IndicesType>
ykuroda 0:13a5d365ba16 512 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
ykuroda 0:13a5d365ba16 513 {
ykuroda 0:13a5d365ba16 514 typedef PermutationBase<PermutationWrapper> Base;
ykuroda 0:13a5d365ba16 515 typedef internal::traits<PermutationWrapper> Traits;
ykuroda 0:13a5d365ba16 516 public:
ykuroda 0:13a5d365ba16 517
ykuroda 0:13a5d365ba16 518 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 519 typedef typename Traits::IndicesType IndicesType;
ykuroda 0:13a5d365ba16 520 #endif
ykuroda 0:13a5d365ba16 521
ykuroda 0:13a5d365ba16 522 inline PermutationWrapper(const IndicesType& a_indices)
ykuroda 0:13a5d365ba16 523 : m_indices(a_indices)
ykuroda 0:13a5d365ba16 524 {}
ykuroda 0:13a5d365ba16 525
ykuroda 0:13a5d365ba16 526 /** const version of indices(). */
ykuroda 0:13a5d365ba16 527 const typename internal::remove_all<typename IndicesType::Nested>::type&
ykuroda 0:13a5d365ba16 528 indices() const { return m_indices; }
ykuroda 0:13a5d365ba16 529
ykuroda 0:13a5d365ba16 530 protected:
ykuroda 0:13a5d365ba16 531
ykuroda 0:13a5d365ba16 532 typename IndicesType::Nested m_indices;
ykuroda 0:13a5d365ba16 533 };
ykuroda 0:13a5d365ba16 534
ykuroda 0:13a5d365ba16 535 /** \returns the matrix with the permutation applied to the columns.
ykuroda 0:13a5d365ba16 536 */
ykuroda 0:13a5d365ba16 537 template<typename Derived, typename PermutationDerived>
ykuroda 0:13a5d365ba16 538 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
ykuroda 0:13a5d365ba16 539 operator*(const MatrixBase<Derived>& matrix,
ykuroda 0:13a5d365ba16 540 const PermutationBase<PermutationDerived> &permutation)
ykuroda 0:13a5d365ba16 541 {
ykuroda 0:13a5d365ba16 542 return internal::permut_matrix_product_retval
ykuroda 0:13a5d365ba16 543 <PermutationDerived, Derived, OnTheRight>
ykuroda 0:13a5d365ba16 544 (permutation.derived(), matrix.derived());
ykuroda 0:13a5d365ba16 545 }
ykuroda 0:13a5d365ba16 546
ykuroda 0:13a5d365ba16 547 /** \returns the matrix with the permutation applied to the rows.
ykuroda 0:13a5d365ba16 548 */
ykuroda 0:13a5d365ba16 549 template<typename Derived, typename PermutationDerived>
ykuroda 0:13a5d365ba16 550 inline const internal::permut_matrix_product_retval
ykuroda 0:13a5d365ba16 551 <PermutationDerived, Derived, OnTheLeft>
ykuroda 0:13a5d365ba16 552 operator*(const PermutationBase<PermutationDerived> &permutation,
ykuroda 0:13a5d365ba16 553 const MatrixBase<Derived>& matrix)
ykuroda 0:13a5d365ba16 554 {
ykuroda 0:13a5d365ba16 555 return internal::permut_matrix_product_retval
ykuroda 0:13a5d365ba16 556 <PermutationDerived, Derived, OnTheLeft>
ykuroda 0:13a5d365ba16 557 (permutation.derived(), matrix.derived());
ykuroda 0:13a5d365ba16 558 }
ykuroda 0:13a5d365ba16 559
ykuroda 0:13a5d365ba16 560 namespace internal {
ykuroda 0:13a5d365ba16 561
ykuroda 0:13a5d365ba16 562 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
ykuroda 0:13a5d365ba16 563 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
ykuroda 0:13a5d365ba16 564 {
ykuroda 0:13a5d365ba16 565 typedef typename MatrixType::PlainObject ReturnType;
ykuroda 0:13a5d365ba16 566 };
ykuroda 0:13a5d365ba16 567
ykuroda 0:13a5d365ba16 568 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
ykuroda 0:13a5d365ba16 569 struct permut_matrix_product_retval
ykuroda 0:13a5d365ba16 570 : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
ykuroda 0:13a5d365ba16 571 {
ykuroda 0:13a5d365ba16 572 typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
ykuroda 0:13a5d365ba16 573 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 574
ykuroda 0:13a5d365ba16 575 permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
ykuroda 0:13a5d365ba16 576 : m_permutation(perm), m_matrix(matrix)
ykuroda 0:13a5d365ba16 577 {}
ykuroda 0:13a5d365ba16 578
ykuroda 0:13a5d365ba16 579 inline Index rows() const { return m_matrix.rows(); }
ykuroda 0:13a5d365ba16 580 inline Index cols() const { return m_matrix.cols(); }
ykuroda 0:13a5d365ba16 581
ykuroda 0:13a5d365ba16 582 template<typename Dest> inline void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 583 {
ykuroda 0:13a5d365ba16 584 const Index n = Side==OnTheLeft ? rows() : cols();
ykuroda 0:13a5d365ba16 585 // FIXME we need an is_same for expression that is not sensitive to constness. For instance
ykuroda 0:13a5d365ba16 586 // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
ykuroda 0:13a5d365ba16 587 if( is_same<MatrixTypeNestedCleaned,Dest>::value
ykuroda 0:13a5d365ba16 588 && blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
ykuroda 0:13a5d365ba16 589 && blas_traits<Dest>::HasUsableDirectAccess
ykuroda 0:13a5d365ba16 590 && extract_data(dst) == extract_data(m_matrix))
ykuroda 0:13a5d365ba16 591 {
ykuroda 0:13a5d365ba16 592 // apply the permutation inplace
ykuroda 0:13a5d365ba16 593 Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
ykuroda 0:13a5d365ba16 594 mask.fill(false);
ykuroda 0:13a5d365ba16 595 Index r = 0;
ykuroda 0:13a5d365ba16 596 while(r < m_permutation.size())
ykuroda 0:13a5d365ba16 597 {
ykuroda 0:13a5d365ba16 598 // search for the next seed
ykuroda 0:13a5d365ba16 599 while(r<m_permutation.size() && mask[r]) r++;
ykuroda 0:13a5d365ba16 600 if(r>=m_permutation.size())
ykuroda 0:13a5d365ba16 601 break;
ykuroda 0:13a5d365ba16 602 // we got one, let's follow it until we are back to the seed
ykuroda 0:13a5d365ba16 603 Index k0 = r++;
ykuroda 0:13a5d365ba16 604 Index kPrev = k0;
ykuroda 0:13a5d365ba16 605 mask.coeffRef(k0) = true;
ykuroda 0:13a5d365ba16 606 for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
ykuroda 0:13a5d365ba16 607 {
ykuroda 0:13a5d365ba16 608 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
ykuroda 0:13a5d365ba16 609 .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
ykuroda 0:13a5d365ba16 610 (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
ykuroda 0:13a5d365ba16 611
ykuroda 0:13a5d365ba16 612 mask.coeffRef(k) = true;
ykuroda 0:13a5d365ba16 613 kPrev = k;
ykuroda 0:13a5d365ba16 614 }
ykuroda 0:13a5d365ba16 615 }
ykuroda 0:13a5d365ba16 616 }
ykuroda 0:13a5d365ba16 617 else
ykuroda 0:13a5d365ba16 618 {
ykuroda 0:13a5d365ba16 619 for(int i = 0; i < n; ++i)
ykuroda 0:13a5d365ba16 620 {
ykuroda 0:13a5d365ba16 621 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
ykuroda 0:13a5d365ba16 622 (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
ykuroda 0:13a5d365ba16 623
ykuroda 0:13a5d365ba16 624 =
ykuroda 0:13a5d365ba16 625
ykuroda 0:13a5d365ba16 626 Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
ykuroda 0:13a5d365ba16 627 (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
ykuroda 0:13a5d365ba16 628 }
ykuroda 0:13a5d365ba16 629 }
ykuroda 0:13a5d365ba16 630 }
ykuroda 0:13a5d365ba16 631
ykuroda 0:13a5d365ba16 632 protected:
ykuroda 0:13a5d365ba16 633 const PermutationType& m_permutation;
ykuroda 0:13a5d365ba16 634 typename MatrixType::Nested m_matrix;
ykuroda 0:13a5d365ba16 635 };
ykuroda 0:13a5d365ba16 636
ykuroda 0:13a5d365ba16 637 /* Template partial specialization for transposed/inverse permutations */
ykuroda 0:13a5d365ba16 638
ykuroda 0:13a5d365ba16 639 template<typename Derived>
ykuroda 0:13a5d365ba16 640 struct traits<Transpose<PermutationBase<Derived> > >
ykuroda 0:13a5d365ba16 641 : traits<Derived>
ykuroda 0:13a5d365ba16 642 {};
ykuroda 0:13a5d365ba16 643
ykuroda 0:13a5d365ba16 644 } // end namespace internal
ykuroda 0:13a5d365ba16 645
ykuroda 0:13a5d365ba16 646 template<typename Derived>
ykuroda 0:13a5d365ba16 647 class Transpose<PermutationBase<Derived> >
ykuroda 0:13a5d365ba16 648 : public EigenBase<Transpose<PermutationBase<Derived> > >
ykuroda 0:13a5d365ba16 649 {
ykuroda 0:13a5d365ba16 650 typedef Derived PermutationType;
ykuroda 0:13a5d365ba16 651 typedef typename PermutationType::IndicesType IndicesType;
ykuroda 0:13a5d365ba16 652 typedef typename PermutationType::PlainPermutationType PlainPermutationType;
ykuroda 0:13a5d365ba16 653 public:
ykuroda 0:13a5d365ba16 654
ykuroda 0:13a5d365ba16 655 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 656 typedef internal::traits<PermutationType> Traits;
ykuroda 0:13a5d365ba16 657 typedef typename Derived::DenseMatrixType DenseMatrixType;
ykuroda 0:13a5d365ba16 658 enum {
ykuroda 0:13a5d365ba16 659 Flags = Traits::Flags,
ykuroda 0:13a5d365ba16 660 CoeffReadCost = Traits::CoeffReadCost,
ykuroda 0:13a5d365ba16 661 RowsAtCompileTime = Traits::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 662 ColsAtCompileTime = Traits::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 663 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 664 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 665 };
ykuroda 0:13a5d365ba16 666 typedef typename Traits::Scalar Scalar;
ykuroda 0:13a5d365ba16 667 #endif
ykuroda 0:13a5d365ba16 668
ykuroda 0:13a5d365ba16 669 Transpose(const PermutationType& p) : m_permutation(p) {}
ykuroda 0:13a5d365ba16 670
ykuroda 0:13a5d365ba16 671 inline int rows() const { return m_permutation.rows(); }
ykuroda 0:13a5d365ba16 672 inline int cols() const { return m_permutation.cols(); }
ykuroda 0:13a5d365ba16 673
ykuroda 0:13a5d365ba16 674 #ifndef EIGEN_PARSED_BY_DOXYGEN
ykuroda 0:13a5d365ba16 675 template<typename DenseDerived>
ykuroda 0:13a5d365ba16 676 void evalTo(MatrixBase<DenseDerived>& other) const
ykuroda 0:13a5d365ba16 677 {
ykuroda 0:13a5d365ba16 678 other.setZero();
ykuroda 0:13a5d365ba16 679 for (int i=0; i<rows();++i)
ykuroda 0:13a5d365ba16 680 other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
ykuroda 0:13a5d365ba16 681 }
ykuroda 0:13a5d365ba16 682 #endif
ykuroda 0:13a5d365ba16 683
ykuroda 0:13a5d365ba16 684 /** \return the equivalent permutation matrix */
ykuroda 0:13a5d365ba16 685 PlainPermutationType eval() const { return *this; }
ykuroda 0:13a5d365ba16 686
ykuroda 0:13a5d365ba16 687 DenseMatrixType toDenseMatrix() const { return *this; }
ykuroda 0:13a5d365ba16 688
ykuroda 0:13a5d365ba16 689 /** \returns the matrix with the inverse permutation applied to the columns.
ykuroda 0:13a5d365ba16 690 */
ykuroda 0:13a5d365ba16 691 template<typename OtherDerived> friend
ykuroda 0:13a5d365ba16 692 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
ykuroda 0:13a5d365ba16 693 operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
ykuroda 0:13a5d365ba16 694 {
ykuroda 0:13a5d365ba16 695 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
ykuroda 0:13a5d365ba16 696 }
ykuroda 0:13a5d365ba16 697
ykuroda 0:13a5d365ba16 698 /** \returns the matrix with the inverse permutation applied to the rows.
ykuroda 0:13a5d365ba16 699 */
ykuroda 0:13a5d365ba16 700 template<typename OtherDerived>
ykuroda 0:13a5d365ba16 701 inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
ykuroda 0:13a5d365ba16 702 operator*(const MatrixBase<OtherDerived>& matrix) const
ykuroda 0:13a5d365ba16 703 {
ykuroda 0:13a5d365ba16 704 return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
ykuroda 0:13a5d365ba16 705 }
ykuroda 0:13a5d365ba16 706
ykuroda 0:13a5d365ba16 707 const PermutationType& nestedPermutation() const { return m_permutation; }
ykuroda 0:13a5d365ba16 708
ykuroda 0:13a5d365ba16 709 protected:
ykuroda 0:13a5d365ba16 710 const PermutationType& m_permutation;
ykuroda 0:13a5d365ba16 711 };
ykuroda 0:13a5d365ba16 712
ykuroda 0:13a5d365ba16 713 template<typename Derived>
ykuroda 0:13a5d365ba16 714 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
ykuroda 0:13a5d365ba16 715 {
ykuroda 0:13a5d365ba16 716 return derived();
ykuroda 0:13a5d365ba16 717 }
ykuroda 0:13a5d365ba16 718
ykuroda 0:13a5d365ba16 719 } // end namespace Eigen
ykuroda 0:13a5d365ba16 720
ykuroda 0:13a5d365ba16 721 #endif // EIGEN_PERMUTATIONMATRIX_H