mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #632: add specializations for res ?= dense +/- sparse and res ?= sparse +/- dense.

Browse Source 
They are rewritten as two compound assignment to by-pass hybrid dense-sparse iterator.
This commit is contained in:
Gael Guennebaud 2018-10-10 22:50:15 +02:00
parent 76ceae49c1
commit eec0dfd688
2 changed files with 95 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -134,8 +134,8 @@ struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Sparse>
};
// Generic Sparse to Dense assignment
template< typename DstXprType, typename SrcXprType, typename Functor>
struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense>
template< typename DstXprType, typename SrcXprType, typename Functor, typename Weak>
struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense, Weak>
{
static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
{
@ -153,6 +153,73 @@ struct Assignment<DstXprType, SrcXprType, Functor, Sparse2Dense>
}
};
// Specialization for dense ?= dense +/- sparse and dense ?= sparse +/- dense
template<typename DstXprType, typename Func1, typename Func2>
struct assignment_from_dense_op_sparse
{
template<typename SrcXprType, typename InitialFunc>
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
{
#ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN
EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN
#endif
call_assignment_no_alias(dst, src.lhs(), Func1());
call_assignment_no_alias(dst, src.rhs(), Func2());
}
// Specialization for dense1 = sparse + dense2; -> dense1 = dense2; dense1 += sparse;
template<typename Lhs, typename Rhs, typename Scalar>
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type
run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_sum_op<Scalar,Scalar>, const Lhs, const Rhs> &src,
const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/)
{
#ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN
EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN
#endif
// Apply the dense matrix first, then the sparse one.
call_assignment_no_alias(dst, src.rhs(), Func1());
call_assignment_no_alias(dst, src.lhs(), Func2());
}
// Specialization for dense1 = sparse - dense2; -> dense1 = -dense2; dense1 += sparse;
template<typename Lhs, typename Rhs, typename Scalar>
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
typename internal::enable_if<internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type
run(DstXprType &dst, const CwiseBinaryOp<internal::scalar_difference_op<Scalar,Scalar>, const Lhs, const Rhs> &src,
const internal::assign_op<typename DstXprType::Scalar,Scalar>& /*func*/)
{
#ifdef EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN
EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN
#endif
// Apply the dense matrix first, then the sparse one.
call_assignment_no_alias(dst, -src.rhs(), Func1());
call_assignment_no_alias(dst, src.lhs(), add_assign_op<typename DstXprType::Scalar,typename Lhs::Scalar>());
}
};
#define EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(ASSIGN_OP,BINOP,ASSIGN_OP2) \
template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> \
struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<Scalar,Scalar>, const Lhs, const Rhs>, internal::ASSIGN_OP<typename DstXprType::Scalar,Scalar>, \
Sparse2Dense, \
typename internal::enable_if< internal::is_same<typename internal::evaluator_traits<Lhs>::Shape,DenseShape>::value \
|| internal::is_same<typename internal::evaluator_traits<Rhs>::Shape,DenseShape>::value>::type> \
: assignment_from_dense_op_sparse<DstXprType, internal::ASSIGN_OP<typename DstXprType::Scalar,typename Lhs::Scalar>, internal::ASSIGN_OP2<typename DstXprType::Scalar,typename Rhs::Scalar> > \
{}
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op, scalar_sum_op,add_assign_op);
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_sum_op,add_assign_op);
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_sum_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(assign_op, scalar_difference_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(add_assign_op,scalar_difference_op,sub_assign_op);
EIGEN_CATCH_ASSIGN_DENSE_OP_SPARSE(sub_assign_op,scalar_difference_op,add_assign_op);
// Specialization for "dst = dec.solve(rhs)"
// NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error
template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>

26
test/sparse_basic.cpp
View File

@ -12,6 +12,11 @@
static long g_realloc_count = 0;
#define EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN g_realloc_count++;
static long g_dense_op_sparse_count = 0;
#define EIGEN_SPARSE_ASSIGNMENT_FROM_DENSE_OP_SPARSE_PLUGIN g_dense_op_sparse_count++;
#define EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_ADD_DENSE_PLUGIN g_dense_op_sparse_count+=10;
#define EIGEN_SPARSE_ASSIGNMENT_FROM_SPARSE_SUB_DENSE_PLUGIN g_dense_op_sparse_count+=20;
#include "sparse.h"
template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
@ -194,6 +199,7 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
VERIFY_IS_APPROX(refM4.cwiseProduct(m3), refM4.cwiseProduct(refM3));
// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
// mixed sparse-dense
VERIFY_IS_APPROX(refM4 + m3, refM4 + refM3);
VERIFY_IS_APPROX(m3 + refM4, refM3 + refM4);
VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3);
@ -222,6 +228,26 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
refM3 = refM1;
VERIFY_IS_APPROX(refM1+=m2, refM3+=refM2);
VERIFY_IS_APPROX(refM1-=m2, refM3-=refM2);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1 =m2+refM4, refM3 =refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,10);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1+=m2+refM4, refM3+=refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1-=m2+refM4, refM3-=refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1 =refM4+m2, refM3 =refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1+=refM4+m2, refM3+=refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1-=refM4+m2, refM3-=refM2+refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1 =m2-refM4, refM3 =refM2-refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,20);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1+=m2-refM4, refM3+=refM2-refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1-=m2-refM4, refM3-=refM2-refM4); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1 =refM4-m2, refM3 =refM4-refM2); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1+=refM4-m2, refM3+=refM4-refM2); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
g_dense_op_sparse_count=0; VERIFY_IS_APPROX(refM1-=refM4-m2, refM3-=refM4-refM2); VERIFY_IS_EQUAL(g_dense_op_sparse_count,1);
refM3 = m3;
if (rows>=2 && cols>=2)
{
VERIFY_RAISES_ASSERT( m1 += m1.innerVector(0) );