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

add mixed sparse-dense outer product

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-25 11:36:38 +02:00
parent 1927b4dff5
commit e313826890
4 changed files with 129 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

116
Eigen/src/Sparse/SparseDenseProduct.h
View File

@ -25,6 +25,118 @@
#ifndef EIGEN_SPARSEDENSEPRODUCT_H
#define EIGEN_SPARSEDENSEPRODUCT_H
template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductReturnType
{
typedef SparseTimeDenseProduct<Lhs,Rhs> Type;
};
template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
{
typedef SparseDenseOuterProduct<Lhs,Rhs,false> Type;
};
template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
{
typedef DenseTimeSparseProduct<Lhs,Rhs> Type;
};
template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
{
typedef SparseDenseOuterProduct<Rhs,Lhs,true> Type;
};
template<typename Lhs, typename Rhs, bool Tr>
struct ei_traits<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
{
typedef Sparse StorageKind;
typedef typename ei_scalar_product_traits<typename ei_traits<Lhs>::Scalar,
typename ei_traits<Rhs>::Scalar>::ReturnType Scalar;
typedef typename Lhs::Index Index;
typedef typename Lhs::Nested LhsNested;
typedef typename Rhs::Nested RhsNested;
typedef typename ei_cleantype<LhsNested>::type _LhsNested;
typedef typename ei_cleantype<RhsNested>::type _RhsNested;
enum {
LhsCoeffReadCost = ei_traits<_LhsNested>::CoeffReadCost,
RhsCoeffReadCost = ei_traits<_RhsNested>::CoeffReadCost,
RowsAtCompileTime = Tr ? int(ei_traits<Rhs>::RowsAtCompileTime) : int(ei_traits<Lhs>::RowsAtCompileTime),
ColsAtCompileTime = Tr ? int(ei_traits<Lhs>::ColsAtCompileTime) : int(ei_traits<Rhs>::ColsAtCompileTime),
MaxRowsAtCompileTime = Tr ? int(ei_traits<Rhs>::MaxRowsAtCompileTime) : int(ei_traits<Lhs>::MaxRowsAtCompileTime),
MaxColsAtCompileTime = Tr ? int(ei_traits<Lhs>::MaxColsAtCompileTime) : int(ei_traits<Rhs>::MaxColsAtCompileTime),
Flags = Tr ? RowMajorBit : 0,
CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + NumTraits<Scalar>::MulCost
};
};
template<typename Lhs, typename Rhs, bool Tr>
class SparseDenseOuterProduct
: public SparseMatrixBase<SparseDenseOuterProduct<Lhs,Rhs,Tr> >
{
public:
typedef SparseMatrixBase<SparseDenseOuterProduct> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(SparseDenseOuterProduct)
typedef ei_traits<SparseDenseOuterProduct> Traits;
private:
typedef typename Traits::LhsNested LhsNested;
typedef typename Traits::RhsNested RhsNested;
typedef typename Traits::_LhsNested _LhsNested;
typedef typename Traits::_RhsNested _RhsNested;
public:
class InnerIterator;
EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Lhs& lhs, const Rhs& rhs)
: m_lhs(lhs), m_rhs(rhs)
{
EIGEN_STATIC_ASSERT(!Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
}
EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Rhs& rhs, const Lhs& lhs)
: m_lhs(lhs), m_rhs(rhs)
{
EIGEN_STATIC_ASSERT(Tr,YOU_MADE_A_PROGRAMMING_MISTAKE);
}
EIGEN_STRONG_INLINE Index rows() const { return Tr ? m_rhs.rows() : m_lhs.rows(); }
EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : m_rhs.cols(); }
EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
protected:
LhsNested m_lhs;
RhsNested m_rhs;
};
template<typename Lhs, typename Rhs, bool Transpose>
class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNested::InnerIterator
{
typedef typename _LhsNested::InnerIterator Base;
public:
EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer)
: Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer))
{
}
inline Index outer() const { return m_outer; }
inline Index row() const { return Transpose ? Base::row() : m_outer; }
inline Index col() const { return Transpose ? m_outer : Base::row(); }
inline Scalar value() const { return Base::value() * m_factor; }
protected:
int m_outer;
Scalar m_factor;
};
template<typename Lhs, typename Rhs>
struct ei_traits<SparseTimeDenseProduct<Lhs,Rhs> >
: ei_traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
@ -102,10 +214,10 @@ class DenseTimeSparseProduct
// sparse * dense
template<typename Derived>
template<typename OtherDerived>
inline const SparseTimeDenseProduct<Derived,OtherDerived>
inline const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
{
return SparseTimeDenseProduct<Derived,OtherDerived>(derived(), other.derived());
return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
}
#endif // EIGEN_SPARSEDENSEPRODUCT_H

10
Eigen/src/Sparse/SparseMatrixBase.h
View File

@ -362,15 +362,15 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
{ return SparseDiagonalProduct<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
// dense * sparse (return a dense object)
/** dense * sparse (return a dense object unless it is an outer product) */
template<typename OtherDerived> friend
const DenseTimeSparseProduct<OtherDerived,Derived>
const typename DenseSparseProductReturnType<OtherDerived,Derived>::Type
operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs)
{ return DenseTimeSparseProduct<OtherDerived,Derived>(lhs.derived(),rhs); }
{ return typename DenseSparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); }
// sparse * dense (returns a dense object)
/** sparse * dense (returns a dense object unless it is an outer product) */
template<typename OtherDerived>
const SparseTimeDenseProduct<Derived,OtherDerived>
const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
operator*(const MatrixBase<OtherDerived> &other) const;
template<typename OtherDerived>

3
Eigen/src/Sparse/SparseUtil.h
View File

@ -95,8 +95,11 @@ template<typename MatrixType> class SparseView;
template<typename Lhs, typename Rhs> class SparseSparseProduct;
template<typename Lhs, typename Rhs> class SparseTimeDenseProduct;
template<typename Lhs, typename Rhs> class DenseTimeSparseProduct;
template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct;
template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType;
template<typename Lhs, typename Rhs, int InnerSize = ei_traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType;
template<typename Lhs, typename Rhs, int InnerSize = ei_traits<Lhs>::ColsAtCompileTime> struct SparseDenseProductReturnType;
template<typename T> struct ei_eval<T,Sparse>
{

7
test/sparse_product.cpp
View File

@ -46,6 +46,7 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
DenseMatrix refMat5 = DenseMatrix::Random(rows, rows);
DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
DenseVector dv1 = DenseVector::Random(rows);
SparseMatrixType m2(rows, rows);
SparseMatrixType m3(rows, rows);
SparseMatrixType m4(rows, rows);
@ -53,6 +54,8 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
initSparse<Scalar>(density, refMat3, m3);
initSparse<Scalar>(density, refMat4, m4);
int c = ei_random<int>(0,rows-1);
VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
@ -77,6 +80,10 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
// sparse * dense and dense * sparse outer product
VERIFY_IS_APPROX(m4=m2.col(c)*dv1.transpose(), refMat4=refMat2.col(c)*dv1.transpose());
VERIFY_IS_APPROX(m4=dv1*m2.col(c).transpose(), refMat4=dv1*refMat2.col(c).transpose());
VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3);
}