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

Make SparseSelfAdjointView, twists, and SparseQR more evaluator friendly

Browse Source
This commit is contained in:
Gael Guennebaud 2015-06-24 17:54:09 +02:00
parent 36643eec0c
commit 763c833637
4 changed files with 166 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Eigen/src/SparseCore/SparseAssign.h
View File

@ -16,8 +16,7 @@ template<typename Derived>
template<typename OtherDerived> template<typename OtherDerived>
Derived& SparseMatrixBase<Derived>::operator=(const EigenBase<OtherDerived> &other) Derived& SparseMatrixBase<Derived>::operator=(const EigenBase<OtherDerived> &other)
{ {
// TODO use the evaluator mechanism internal::call_assignment_no_alias(derived(), other.derived());
other.derived().evalTo(derived());
return derived(); return derived();
} }
@ -137,8 +136,8 @@ struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Sparse, Scalar>
}; };
// Sparse to Dense assignment // Sparse to Dense assignment
template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> template< typename DstXprType, typename SrcXprType, typename Functor>
struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Scalar> struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense>
{ {
static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
{ {
@ -153,8 +152,8 @@ struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Scalar>
} }
}; };
template< typename DstXprType, typename SrcXprType, typename Scalar> template< typename DstXprType, typename SrcXprType>
struct Assignment<DstXprType, SrcXprType, internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense, Scalar> struct Assignment<DstXprType, SrcXprType, internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense>
{ {
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &) static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &)
{ {
@ -173,7 +172,7 @@ struct Assignment<DstXprType, SrcXprType, internal::assign_op<typename DstXprTyp
// Specialization for "dst = dec.solve(rhs)" // Specialization for "dst = dec.solve(rhs)"
// NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error
template<typename DstXprType, typename DecType, typename RhsType, typename Scalar> template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar>, Sparse2Sparse, Scalar> struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar>, Sparse2Sparse>
{ {
typedef Solve<DecType,RhsType> SrcXprType; typedef Solve<DecType,RhsType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &) static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)

16
Eigen/src/SparseCore/SparseCwiseBinaryOp.h
View File

@ -390,6 +390,22 @@ SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& othe
return derived() = derived() + other.derived(); return derived() = derived() + other.derived();
} }
template<typename Derived>
template<typename OtherDerived>
Derived& SparseMatrixBase<Derived>::operator+=(const DiagonalBase<OtherDerived>& other)
{
call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar>());
return derived();
}
template<typename Derived>
template<typename OtherDerived>
Derived& SparseMatrixBase<Derived>::operator-=(const DiagonalBase<OtherDerived>& other)
{
call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar>());
return derived();
}
template<typename Derived> template<typename Derived>
template<typename OtherDerived> template<typename OtherDerived>
EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE

131
Eigen/src/SparseCore/SparseSelfAdjointView.h
View File

@ -45,8 +45,13 @@ template<typename MatrixType, unsigned int _Mode> class SparseSelfAdjointView
{ {
public: public:
enum { Mode = _Mode }; enum {
Mode = _Mode,
RowsAtCompileTime = internal::traits<SparseSelfAdjointView>::RowsAtCompileTime,
ColsAtCompileTime = internal::traits<SparseSelfAdjointView>::ColsAtCompileTime
};
typedef EigenBase<SparseSelfAdjointView> Base;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::StorageIndex StorageIndex;
typedef Matrix<StorageIndex,Dynamic,1> VectorI; typedef Matrix<StorageIndex,Dynamic,1> VectorI;
@ -116,20 +121,6 @@ template<typename MatrixType, unsigned int _Mode> class SparseSelfAdjointView
template<typename DerivedU> template<typename DerivedU>
SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
/** \internal triggered by sparse_matrix = SparseSelfadjointView; */
template<typename DestScalar,int StorageOrder> void evalTo(SparseMatrix<DestScalar,StorageOrder,StorageIndex>& _dest) const
{
internal::permute_symm_to_fullsymm<Mode>(m_matrix, _dest);
}
template<typename DestScalar> void evalTo(DynamicSparseMatrix<DestScalar,ColMajor,StorageIndex>& _dest) const
{
// TODO directly evaluate into _dest;
SparseMatrix<DestScalar,ColMajor,StorageIndex> tmp(_dest.rows(),_dest.cols());
internal::permute_symm_to_fullsymm<Mode>(m_matrix, tmp);
_dest = tmp;
}
/** \returns an expression of P H P^-1 */ /** \returns an expression of P H P^-1 */
// TODO implement twists in a more evaluator friendly fashion // TODO implement twists in a more evaluator friendly fashion
SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm) const SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm) const
@ -140,7 +131,7 @@ template<typename MatrixType, unsigned int _Mode> class SparseSelfAdjointView
template<typename SrcMatrixType,int SrcMode> template<typename SrcMatrixType,int SrcMode>
SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcMode>& permutedMatrix) SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct<SrcMatrixType,SrcMode>& permutedMatrix)
{ {
permutedMatrix.evalTo(*this); internal::call_assignment_no_alias_no_transpose(*this, permutedMatrix);
return *this; return *this;
} }
@ -157,11 +148,21 @@ template<typename MatrixType, unsigned int _Mode> class SparseSelfAdjointView
return *this = src.twistedBy(pnull); return *this = src.twistedBy(pnull);
} }
void resize(Index rows, Index cols)
{
EIGEN_ONLY_USED_FOR_DEBUG(rows);
EIGEN_ONLY_USED_FOR_DEBUG(cols);
eigen_assert(rows == this->rows() && cols == this->cols()
&& "SparseSelfadjointView::resize() does not actually allow to resize.");
}
protected: protected:
typename MatrixType::Nested m_matrix; typename MatrixType::Nested m_matrix;
//mutable VectorI m_countPerRow; //mutable VectorI m_countPerRow;
//mutable VectorI m_countPerCol; //mutable VectorI m_countPerCol;
private:
template<typename Dest> void evalTo(Dest &) const;
}; };
/*************************************************************************** /***************************************************************************
@ -200,6 +201,47 @@ SparseSelfAdjointView<MatrixType,Mode>::rankUpdate(const SparseMatrixBase<Derive
return *this; return *this;
} }
namespace internal {
// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.>
// in the future selfadjoint-ness should be defined by the expression traits
// such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
template<typename MatrixType, unsigned int Mode>
struct evaluator_traits<SparseSelfAdjointView<MatrixType,Mode> >
{
typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
typedef SparseSelfAdjointShape Shape;
static const int AssumeAliasing = 0;
};
struct SparseSelfAdjoint2Sparse {};
template<> struct AssignmentKind<SparseShape,SparseSelfAdjointShape> { typedef SparseSelfAdjoint2Sparse Kind; };
template<> struct AssignmentKind<SparseSelfAdjointShape,SparseShape> { typedef Sparse2Sparse Kind; };
template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
struct Assignment<DstXprType, SrcXprType, Functor, SparseSelfAdjoint2Sparse, Scalar>
{
typedef typename DstXprType::StorageIndex StorageIndex;
template<typename DestScalar,int StorageOrder>
static void run(SparseMatrix<DestScalar,StorageOrder,StorageIndex> &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &/*func*/)
{
internal::permute_symm_to_fullsymm<SrcXprType::Mode>(src.matrix(), dst);
}
template<typename DestScalar>
static void run(DynamicSparseMatrix<DestScalar,ColMajor,StorageIndex>& dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar> &/*func*/)
{
// TODO directly evaluate into dst;
SparseMatrix<DestScalar,ColMajor,StorageIndex> tmp(dst.rows(),dst.cols());
internal::permute_symm_to_fullsymm<SrcXprType::Mode>(src.matrix(), tmp);
dst = tmp;
}
};
} // end namespace internal
/*************************************************************************** /***************************************************************************
* Implementation of sparse self-adjoint time dense matrix * Implementation of sparse self-adjoint time dense matrix
***************************************************************************/ ***************************************************************************/
@ -253,17 +295,6 @@ inline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, cons
} }
} }
// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.>
// in the future selfadjoint-ness should be defined by the expression traits
// such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
template<typename MatrixType, unsigned int Mode>
struct evaluator_traits<SparseSelfAdjointView<MatrixType,Mode> >
{
typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
typedef SparseSelfAdjointShape Shape;
static const int AssumeAliasing = 0;
};
template<typename LhsView, typename Rhs, int ProductType> template<typename LhsView, typename Rhs, int ProductType>
struct generic_product_impl<LhsView, Rhs, SparseSelfAdjointShape, DenseShape, ProductType> struct generic_product_impl<LhsView, Rhs, SparseSelfAdjointShape, DenseShape, ProductType>
@ -519,12 +550,16 @@ class SparseSymmetricPermutationProduct
public: public:
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::StorageIndex StorageIndex;
enum {
RowsAtCompileTime = internal::traits<SparseSymmetricPermutationProduct>::RowsAtCompileTime,
ColsAtCompileTime = internal::traits<SparseSymmetricPermutationProduct>::ColsAtCompileTime
};
protected: protected:
typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> Perm; typedef PermutationMatrix<Dynamic,Dynamic,StorageIndex> Perm;
public: public:
typedef Matrix<StorageIndex,Dynamic,1> VectorI; typedef Matrix<StorageIndex,Dynamic,1> VectorI;
typedef typename MatrixType::Nested MatrixTypeNested; typedef typename MatrixType::Nested MatrixTypeNested;
typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested; typedef typename internal::remove_all<MatrixTypeNested>::type NestedExpression;
SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm)
: m_matrix(mat), m_perm(perm) : m_matrix(mat), m_perm(perm)
@ -533,19 +568,8 @@ class SparseSymmetricPermutationProduct
inline Index rows() const { return m_matrix.rows(); } inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); } inline Index cols() const { return m_matrix.cols(); }
template<typename DestScalar, int Options, typename DstIndex> const NestedExpression& matrix() const { return m_matrix; }
void evalTo(SparseMatrix<DestScalar,Options,DstIndex>& _dest) const const Perm& perm() const { return m_perm; }
{
// internal::permute_symm_to_fullsymm<Mode>(m_matrix,_dest,m_perm.indices().data());
SparseMatrix<DestScalar,(Options&RowMajor)==RowMajor ? ColMajor : RowMajor, DstIndex> tmp;
internal::permute_symm_to_fullsymm<Mode>(m_matrix,tmp,m_perm.indices().data());
_dest = tmp;
}
template<typename DestType,unsigned int DestMode> void evalTo(SparseSelfAdjointView<DestType,DestMode>& dest) const
{
internal::permute_symm_to_symm<Mode,DestMode>(m_matrix,dest.matrix(),m_perm.indices().data());
}
protected: protected:
MatrixTypeNested m_matrix; MatrixTypeNested m_matrix;
@ -553,6 +577,31 @@ class SparseSymmetricPermutationProduct
}; };
namespace internal {
template<typename DstXprType, typename MatrixType, int Mode, typename Scalar>
struct Assignment<DstXprType, SparseSymmetricPermutationProduct<MatrixType,Mode>, internal::assign_op<Scalar>, Sparse2Sparse>
{
typedef SparseSymmetricPermutationProduct<MatrixType,Mode> SrcXprType;
typedef typename DstXprType::StorageIndex DstIndex;
template<int Options>
static void run(SparseMatrix<Scalar,Options,DstIndex> &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
// internal::permute_symm_to_fullsymm<Mode>(m_matrix,_dest,m_perm.indices().data());
SparseMatrix<Scalar,(Options&RowMajor)==RowMajor ? ColMajor : RowMajor, DstIndex> tmp;
internal::permute_symm_to_fullsymm<Mode>(src.matrix(),tmp,src.perm().indices().data());
dst = tmp;
}
template<typename DestType,unsigned int DestMode>
static void run(SparseSelfAdjointView<DestType,DestMode>& dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
internal::permute_symm_to_symm<Mode,DestMode>(src.matrix(),dst.matrix(),src.perm().indices().data());
}
};
} // end namespace internal
} // end namespace Eigen } // end namespace Eigen
#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H

68
Eigen/src/SparseQR/SparseQR.h
View File

@ -23,6 +23,10 @@ namespace internal {
typedef typename SparseQRType::MatrixType ReturnType; typedef typename SparseQRType::MatrixType ReturnType;
typedef typename ReturnType::StorageIndex StorageIndex; typedef typename ReturnType::StorageIndex StorageIndex;
typedef typename ReturnType::StorageKind StorageKind; typedef typename ReturnType::StorageKind StorageKind;
enum {
RowsAtCompileTime = Dynamic,
ColsAtCompileTime = Dynamic
};
}; };
template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> > template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
{ {
@ -235,8 +239,9 @@ class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
return m_info; return m_info;
} }
protected:
inline void sort_matrix_Q() /** \internal */
inline void _sort_matrix_Q()
{ {
if(this->m_isQSorted) return; if(this->m_isQSorted) return;
// The matrix Q is sorted during the transposition // The matrix Q is sorted during the transposition
@ -267,7 +272,6 @@ class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> >
bool m_isEtreeOk; // whether the elimination tree match the initial input matrix bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
template <typename, typename > friend struct SparseQR_QProduct; template <typename, typename > friend struct SparseQR_QProduct;
template <typename > friend struct SparseQRMatrixQReturnType;
}; };
@ -635,6 +639,10 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
{ {
typedef typename SparseQRType::Scalar Scalar; typedef typename SparseQRType::Scalar Scalar;
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
enum {
RowsAtCompileTime = Dynamic,
ColsAtCompileTime = Dynamic
};
explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
template<typename Derived> template<typename Derived>
SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
@ -652,19 +660,6 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
{ {
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
} }
template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const
{
dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows());
}
template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const
{
Dest idMat(m_qr.rows(), m_qr.rows());
idMat.setIdentity();
// Sort the sparse householder reflectors if needed
const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q();
dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false);
}
const SparseQRType& m_qr; const SparseQRType& m_qr;
}; };
@ -680,6 +675,47 @@ struct SparseQRMatrixQTransposeReturnType
const SparseQRType& m_qr; const SparseQRType& m_qr;
}; };
namespace internal {
template<typename SparseQRType>
struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> >
{
typedef typename SparseQRType::MatrixType MatrixType;
typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
typedef SparseShape Shape;
static const int AssumeAliasing = 0;
};
template< typename DstXprType, typename SparseQRType>
struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar>, Sparse2Sparse>
{
typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
typedef typename DstXprType::Scalar Scalar;
typedef typename DstXprType::StorageIndex StorageIndex;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &/*func*/)
{
typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows());
idMat.setIdentity();
// Sort the sparse householder reflectors if needed
const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q();
dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false);
}
};
template< typename DstXprType, typename SparseQRType>
struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar>, Sparse2Dense>
{
typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType;
typedef typename DstXprType::Scalar Scalar;
typedef typename DstXprType::StorageIndex StorageIndex;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &/*func*/)
{
dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows());
}
};
} // end namespace internal
} // end namespace Eigen } // end namespace Eigen
#endif #endif