bug #650: fix sparse * dense wrt noalias and compound assignment

This commit is contained in:
Gael Guennebaud 2015-06-09 18:33:24 +02:00
parent 785b9c0127
commit f899aeb301
2 changed files with 20 additions and 10 deletions

View File

@ -41,7 +41,7 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t
typename Res::Scalar tmp(0); typename Res::Scalar tmp(0);
for(LhsInnerIterator it(lhsEval,j); it ;++it) for(LhsInnerIterator it(lhsEval,j); it ;++it)
tmp += it.value() * rhs.coeff(it.index(),c); tmp += it.value() * rhs.coeff(it.index(),c);
res.coeffRef(j,c) = alpha * tmp; res.coeffRef(j,c) += alpha * tmp;
} }
} }
} }
@ -128,17 +128,18 @@ namespace internal {
template<typename Lhs, typename Rhs, int ProductType> template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType> struct generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SparseShape,DenseShape,ProductType> >
{ {
typedef typename Product<Lhs,Rhs>::Scalar Scalar;
template<typename Dest> template<typename Dest>
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{ {
typedef typename nested_eval<Lhs,Dynamic>::type LhsNested; typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
typedef typename nested_eval<Rhs,Dynamic>::type RhsNested; typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
LhsNested lhsNested(lhs); LhsNested lhsNested(lhs);
RhsNested rhsNested(rhs); RhsNested rhsNested(rhs);
internal::sparse_time_dense_product(lhsNested, rhsNested, dst, alpha);
dst.setZero();
internal::sparse_time_dense_product(lhsNested, rhsNested, dst, typename Dest::Scalar(1));
} }
}; };
@ -149,19 +150,21 @@ struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, DenseShape, Product
template<typename Lhs, typename Rhs, int ProductType> template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType> struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SparseShape,ProductType> >
{ {
template<typename Dest> typedef typename Product<Lhs,Rhs>::Scalar Scalar;
static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
template<typename Dst>
static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
{ {
typedef typename nested_eval<Lhs,Dynamic>::type LhsNested; typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
typedef typename nested_eval<Rhs,Dynamic>::type RhsNested; typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
LhsNested lhsNested(lhs); LhsNested lhsNested(lhs);
RhsNested rhsNested(rhs); RhsNested rhsNested(rhs);
dst.setZero();
// transpose everything // transpose everything
Transpose<Dest> dstT(dst); Transpose<Dst> dstT(dst);
internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, typename Dest::Scalar(1)); internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha);
} }
}; };

View File

@ -88,6 +88,10 @@ template<typename SparseMatrixType> void sparse_product()
VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
VERIFY_IS_APPROX(dm4=dm4+m2*refMat3, refMat4=refMat4+refMat2*refMat3); VERIFY_IS_APPROX(dm4=dm4+m2*refMat3, refMat4=refMat4+refMat2*refMat3);
VERIFY_IS_APPROX(dm4+=m2*refMat3, refMat4+=refMat2*refMat3);
VERIFY_IS_APPROX(dm4-=m2*refMat3, refMat4-=refMat2*refMat3);
VERIFY_IS_APPROX(dm4.noalias()+=m2*refMat3, refMat4+=refMat2*refMat3);
VERIFY_IS_APPROX(dm4.noalias()-=m2*refMat3, refMat4-=refMat2*refMat3);
VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3)); VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3));
VERIFY_IS_APPROX(dm4=m2t.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2t.transpose()*(refMat3+refMat5)*0.5); VERIFY_IS_APPROX(dm4=m2t.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2t.transpose()*(refMat3+refMat5)*0.5);
@ -101,6 +105,9 @@ template<typename SparseMatrixType> void sparse_product()
VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
VERIFY_IS_APPROX(dm4=dm4+refMat2*m3, refMat4=refMat4+refMat2*refMat3); VERIFY_IS_APPROX(dm4=dm4+refMat2*m3, refMat4=refMat4+refMat2*refMat3);
VERIFY_IS_APPROX(dm4+=refMat2*m3, refMat4+=refMat2*refMat3); VERIFY_IS_APPROX(dm4+=refMat2*m3, refMat4+=refMat2*refMat3);
VERIFY_IS_APPROX(dm4-=refMat2*m3, refMat4-=refMat2*refMat3);
VERIFY_IS_APPROX(dm4.noalias()+=refMat2*m3, refMat4+=refMat2*refMat3);
VERIFY_IS_APPROX(dm4.noalias()-=refMat2*m3, refMat4-=refMat2*refMat3);
VERIFY_IS_APPROX(dm4=refMat2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose()); VERIFY_IS_APPROX(dm4=refMat2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose());
VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3); VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose()); VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());