GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 19:59:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Clean implementation of permutation * matrix products.

Browse Source
This commit is contained in:
Gael Guennebaud 2015-06-19 10:51:57 +02:00
parent 06036d8bb1
commit fad36cc814
3 changed files with 108 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
Eigen/src/Core/PermutationMatrix.h
View File

@ -42,10 +42,6 @@ namespace Eigen {
namespace internal { namespace internal {
template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
struct permut_matrix_product_retval;
template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
struct permut_sparsematrix_product_retval;
enum PermPermProduct_t {PermPermProduct}; enum PermPermProduct_t {PermPermProduct};
} // end namespace internal } // end namespace internal
@ -570,80 +566,6 @@ operator*(const PermutationBase<PermutationDerived> &permutation,
namespace internal { namespace internal {
template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
: traits<typename MatrixType::PlainObject>
{
typedef typename MatrixType::PlainObject ReturnType;
};
template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
struct permut_matrix_product_retval
: public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
{
typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
typedef typename MatrixType::StorageIndex StorageIndex;
permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
: m_permutation(perm), m_matrix(matrix)
{}
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{
const Index n = Side==OnTheLeft ? rows() : cols();
// FIXME we need an is_same for expression that is not sensitive to constness. For instance
// is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
//if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
if(is_same_dense(dst, m_matrix))
{
// apply the permutation inplace
Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
mask.fill(false);
Index r = 0;
while(r < m_permutation.size())
{
// search for the next seed
while(r<m_permutation.size() && mask[r]) r++;
if(r>=m_permutation.size())
break;
// we got one, let's follow it until we are back to the seed
Index k0 = r++;
Index kPrev = k0;
mask.coeffRef(k0) = true;
for(Index k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
.swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
mask.coeffRef(k) = true;
kPrev = k;
}
}
}
else
{
for(Index i = 0; i < n; ++i)
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
=
Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
(m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
}
}
}
protected:
const PermutationType& m_permutation;
typename MatrixType::Nested m_matrix;
};
/* Template partial specialization for transposed/inverse permutations */ /* Template partial specialization for transposed/inverse permutations */
template<typename Derived> template<typename Derived>

93
Eigen/src/Core/ProductEvaluators.h
View File

@ -825,48 +825,107 @@ struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DenseShape,
/*************************************************************************** /***************************************************************************
* Products with permutation matrices * Products with permutation matrices
***************************************************************************/ ***************************************************************************/
template<typename Lhs, typename Rhs, int ProductTag> /** \internal
struct generic_product_impl<Lhs, Rhs, PermutationShape, DenseShape, ProductTag> * \class permutation_matrix_product
* Internal helper class implementing the product between a permutation matrix and a matrix.
* This class is specialized for DenseShape below and for SparseShape in SparseCore/SparsePermutation.h
*/
template<typename MatrixType, int Side, bool Transposed, typename MatrixShape>
struct permutation_matrix_product;
template<typename MatrixType, int Side, bool Transposed>
struct permutation_matrix_product<MatrixType, Side, Transposed, DenseShape>
{
typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
template<typename Dest, typename PermutationType>
static inline void run(Dest& dst, const PermutationType& perm, const MatrixType& mat)
{
const Index n = Side==OnTheLeft ? mat.rows() : mat.cols();
// FIXME we need an is_same for expression that is not sensitive to constness. For instance
// is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
//if(is_same<MatrixTypeCleaned,Dest>::value && extract_data(dst) == extract_data(mat))
if(is_same_dense(dst, mat))
{
// apply the permutation inplace
Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(perm.size());
mask.fill(false);
Index r = 0;
while(r < perm.size())
{
// search for the next seed
while(r<perm.size() && mask[r]) r++;
if(r>=perm.size())
break;
// we got one, let's follow it until we are back to the seed
Index k0 = r++;
Index kPrev = k0;
mask.coeffRef(k0) = true;
for(Index k=perm.indices().coeff(k0); k!=k0; k=perm.indices().coeff(k))
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
.swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
mask.coeffRef(k) = true;
kPrev = k;
}
}
}
else
{
for(Index i = 0; i < n; ++i)
{
Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
(dst, ((Side==OnTheLeft) ^ Transposed) ? perm.indices().coeff(i) : i)
=
Block<const MatrixTypeCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
(mat, ((Side==OnTheRight) ^ Transposed) ? perm.indices().coeff(i) : i);
}
}
}
};
template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>
{ {
template<typename Dest> template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{ {
permut_matrix_product_retval<Lhs, Rhs, OnTheLeft, false> pmpr(lhs, rhs); permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
pmpr.evalTo(dst);
} }
}; };
template<typename Lhs, typename Rhs, int ProductTag> template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Rhs, DenseShape, PermutationShape, ProductTag> struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>
{ {
template<typename Dest> template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{ {
permut_matrix_product_retval<Rhs, Lhs, OnTheRight, false> pmpr(rhs, lhs); permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
pmpr.evalTo(dst);
} }
}; };
template<typename Lhs, typename Rhs, int ProductTag> template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Transpose<Lhs>, Rhs, PermutationShape, DenseShape, ProductTag> struct generic_product_impl<Transpose<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>
{ {
template<typename Dest> template<typename Dest>
static void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs) static void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
{ {
permut_matrix_product_retval<Lhs, Rhs, OnTheLeft, true> pmpr(lhs.nestedPermutation(), rhs); permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedPermutation(), rhs);
pmpr.evalTo(dst);
} }
}; };
template<typename Lhs, typename Rhs, int ProductTag> template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
struct generic_product_impl<Lhs, Transpose<Rhs>, DenseShape, PermutationShape, ProductTag> struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, PermutationShape, ProductTag>
{ {
template<typename Dest> template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs) static void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
{ {
permut_matrix_product_retval<Rhs, Lhs, OnTheRight, true> pmpr(rhs.nestedPermutation(), lhs); permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedPermutation(), lhs);
pmpr.evalTo(dst);
} }
}; };

133
Eigen/src/SparseCore/SparsePermutation.h
View File

@ -16,25 +16,8 @@ namespace Eigen {
namespace internal { namespace internal {
template<typename PermutationType, typename MatrixType, int Side, bool Transposed> template<typename MatrixType, int Side, bool Transposed>
struct traits<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> > struct permutation_matrix_product<MatrixType, Side, Transposed, SparseShape>
{
typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
typedef typename MatrixTypeNestedCleaned::StorageIndex StorageIndex;
enum {
SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
};
typedef typename internal::conditional<MoveOuter,
SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>,
SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> >::type ReturnType;
};
template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
struct permut_sparsematrix_product_retval
: public ReturnByValue<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
{ {
typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
typedef typename MatrixTypeNestedCleaned::Scalar Scalar; typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
@ -44,61 +27,55 @@ struct permut_sparsematrix_product_retval
SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
}; };
typedef typename internal::conditional<MoveOuter,
SparseMatrix<Scalar,SrcStorageOrder,StorageIndex>,
SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> >::type ReturnType;
permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix) template<typename Dest,typename PermutationType>
: m_permutation(perm), m_matrix(matrix) static inline void run(Dest& dst, const PermutationType& perm, const MatrixType& mat)
{}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{ {
if(MoveOuter) if(MoveOuter)
{ {
SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(m_matrix.rows(), m_matrix.cols()); SparseMatrix<Scalar,SrcStorageOrder,StorageIndex> tmp(mat.rows(), mat.cols());
Matrix<StorageIndex,Dynamic,1> sizes(m_matrix.outerSize()); Matrix<StorageIndex,Dynamic,1> sizes(mat.outerSize());
for(Index j=0; j<m_matrix.outerSize(); ++j) for(Index j=0; j<mat.outerSize(); ++j)
{ {
Index jp = m_permutation.indices().coeff(j); Index jp = perm.indices().coeff(j);
sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros()); sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = StorageIndex(mat.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros());
} }
tmp.reserve(sizes); tmp.reserve(sizes);
for(Index j=0; j<m_matrix.outerSize(); ++j) for(Index j=0; j<mat.outerSize(); ++j)
{ {
Index jp = m_permutation.indices().coeff(j); Index jp = perm.indices().coeff(j);
Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j; Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j; Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,jsrc); it; ++it) for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,jsrc); it; ++it)
tmp.insertByOuterInner(jdst,it.index()) = it.value(); tmp.insertByOuterInner(jdst,it.index()) = it.value();
} }
dst = tmp; dst = tmp;
} }
else else
{ {
SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(m_matrix.rows(), m_matrix.cols()); SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,StorageIndex> tmp(mat.rows(), mat.cols());
Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize()); Matrix<StorageIndex,Dynamic,1> sizes(tmp.outerSize());
sizes.setZero(); sizes.setZero();
PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm; PermutationMatrix<Dynamic,Dynamic,StorageIndex> perm_cpy;
if((Side==OnTheLeft) ^ Transposed) if((Side==OnTheLeft) ^ Transposed)
perm = m_permutation; perm_cpy = perm;
else else
perm = m_permutation.transpose(); perm_cpy = perm.transpose();
for(Index j=0; j<m_matrix.outerSize(); ++j) for(Index j=0; j<mat.outerSize(); ++j)
for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it) for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,j); it; ++it)
sizes[perm.indices().coeff(it.index())]++; sizes[perm_cpy.indices().coeff(it.index())]++;
tmp.reserve(sizes); tmp.reserve(sizes);
for(Index j=0; j<m_matrix.outerSize(); ++j) for(Index j=0; j<mat.outerSize(); ++j)
for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it) for(typename MatrixTypeNestedCleaned::InnerIterator it(mat,j); it; ++it)
tmp.insertByOuterInner(perm.indices().coeff(it.index()),j) = it.value(); tmp.insertByOuterInner(perm_cpy.indices().coeff(it.index()),j) = it.value();
dst = tmp; dst = tmp;
} }
} }
protected:
const PermutationType& m_permutation;
typename MatrixType::Nested m_matrix;
}; };
} }
@ -107,63 +84,17 @@ namespace internal {
template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; }; template <int ProductTag> struct product_promote_storage_type<Sparse, PermutationStorage, ProductTag> { typedef Sparse ret; };
template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse, ProductTag> { typedef Sparse ret; }; template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse, ProductTag> { typedef Sparse ret; };
// TODO, the following need cleaning, this is just a copy-past of the dense case
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, PermutationShape, SparseShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
permut_sparsematrix_product_retval<Lhs, Rhs, OnTheLeft, false> pmpr(lhs, rhs);
pmpr.evalTo(dst);
}
};
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, SparseShape, PermutationShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
{
permut_sparsematrix_product_retval<Rhs, Lhs, OnTheRight, false> pmpr(rhs, lhs);
pmpr.evalTo(dst);
}
};
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Transpose<Lhs>, Rhs, PermutationShape, SparseShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
{
permut_sparsematrix_product_retval<Lhs, Rhs, OnTheLeft, true> pmpr(lhs.nestedPermutation(), rhs);
pmpr.evalTo(dst);
}
};
template<typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Transpose<Rhs>, SparseShape, PermutationShape, ProductTag>
{
template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
{
permut_sparsematrix_product_retval<Rhs, Lhs, OnTheRight, true> pmpr(rhs.nestedPermutation(), lhs);
pmpr.evalTo(dst);
}
};
// TODO, the following two overloads are only needed to define the right temporary type through // TODO, the following two overloads are only needed to define the right temporary type through
// typename traits<permut_sparsematrix_product_retval<Rhs,Lhs,OnTheRight,false> >::ReturnType // typename traits<permutation_sparse_matrix_product<Rhs,Lhs,OnTheRight,false> >::ReturnType
// while it should be correctly handled by traits<Product<> >::PlainObject // whereas it should be correctly handled by traits<Product<> >::PlainObject
template<typename Lhs, typename Rhs, int ProductTag> template<typename Lhs, typename Rhs, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, PermutationShape, SparseShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar> struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, PermutationShape, SparseShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar>
: public evaluator<typename traits<permut_sparsematrix_product_retval<Lhs,Rhs,OnTheRight,false> >::ReturnType>::type : public evaluator<typename permutation_matrix_product<Rhs,OnTheRight,false,SparseShape>::ReturnType>::type
{ {
typedef Product<Lhs, Rhs, DefaultProduct> XprType; typedef Product<Lhs, Rhs, DefaultProduct> XprType;
typedef typename traits<permut_sparsematrix_product_retval<Lhs,Rhs,OnTheRight,false> >::ReturnType PlainObject; typedef typename permutation_matrix_product<Rhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
typedef typename evaluator<PlainObject>::type Base; typedef typename evaluator<PlainObject>::type Base;
explicit product_evaluator(const XprType& xpr) explicit product_evaluator(const XprType& xpr)
@ -179,10 +110,10 @@ protected:
template<typename Lhs, typename Rhs, int ProductTag> template<typename Lhs, typename Rhs, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, PermutationShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar> struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, PermutationShape, typename traits<Lhs>::Scalar, typename traits<Rhs>::Scalar>
: public evaluator<typename traits<permut_sparsematrix_product_retval<Rhs,Lhs,OnTheRight,false> >::ReturnType>::type : public evaluator<typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType>::type
{ {
typedef Product<Lhs, Rhs, DefaultProduct> XprType; typedef Product<Lhs, Rhs, DefaultProduct> XprType;
typedef typename traits<permut_sparsematrix_product_retval<Rhs,Lhs,OnTheRight,false> >::ReturnType PlainObject; typedef typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
typedef typename evaluator<PlainObject>::type Base; typedef typename evaluator<PlainObject>::type Base;
explicit product_evaluator(const XprType& xpr) explicit product_evaluator(const XprType& xpr)