SparseDenseProduct.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SPARSEDENSEPRODUCT_H
11 #define EIGEN_SPARSEDENSEPRODUCT_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template <> struct product_promote_storage_type<Sparse,Dense, OuterProduct> { typedef Sparse ret; };
18 template <> struct product_promote_storage_type<Dense,Sparse, OuterProduct> { typedef Sparse ret; };
19 
20 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,
21  typename AlphaType,
22  int LhsStorageOrder = ((SparseLhsType::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor,
23  bool ColPerCol = ((DenseRhsType::Flags&RowMajorBit)==0) || DenseRhsType::ColsAtCompileTime==1>
25 
26 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
27 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, true>
28 {
34  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
35  {
36  LhsEval lhsEval(lhs);
37 
38  Index n = lhs.outerSize();
39 #ifdef EIGEN_HAS_OPENMP
41  Index threads = Eigen::nbThreads();
42 #endif
43 
44  for(Index c=0; c<rhs.cols(); ++c)
45  {
46 #ifdef EIGEN_HAS_OPENMP
47  // This 20000 threshold has been found experimentally on 2D and 3D Poisson problems.
48  // It basically represents the minimal amount of work to be done to be worth it.
49  if(threads>1 && lhsEval.nonZerosEstimate() > 20000)
50  {
51  #pragma omp parallel for schedule(dynamic,(n+threads*4-1)/(threads*4)) num_threads(threads)
52  for(Index i=0; i<n; ++i)
53  processRow(lhsEval,rhs,res,alpha,i,c);
54  }
55  else
56 #endif
57  {
58  for(Index i=0; i<n; ++i)
59  processRow(lhsEval,rhs,res,alpha,i,c);
60  }
61  }
62  }
63 
64  static void processRow(const LhsEval& lhsEval, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha, Index i, Index col)
65  {
66  typename Res::Scalar tmp(0);
67  for(LhsInnerIterator it(lhsEval,i); it ;++it)
68  tmp += it.value() * rhs.coeff(it.index(),col);
69  res.coeffRef(i,col) += alpha * tmp;
70  }
71 
72 };
73 
74 // FIXME: what is the purpose of the following specialization? Is it for the BlockedSparse format?
75 // -> let's disable it for now as it is conflicting with generic scalar*matrix and matrix*scalar operators
76 // template<typename T1, typename T2/*, int _Options, typename _StrideType*/>
77 // struct ScalarBinaryOpTraits<T1, Ref<T2/*, _Options, _StrideType*/> >
78 // {
79 // enum {
80 // Defined = 1
81 // };
82 // typedef typename CwiseUnaryOp<scalar_multiple2_op<T1, typename T2::Scalar>, T2>::PlainObject ReturnType;
83 // };
84 
85 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
86 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType, ColMajor, true>
87 {
92  typedef typename LhsEval::InnerIterator LhsInnerIterator;
93  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
94  {
95  LhsEval lhsEval(lhs);
96  for(Index c=0; c<rhs.cols(); ++c)
97  {
98  for(Index j=0; j<lhs.outerSize(); ++j)
99  {
100 // typename Res::Scalar rhs_j = alpha * rhs.coeff(j,c);
101  typename ScalarBinaryOpTraits<AlphaType, typename Rhs::Scalar>::ReturnType rhs_j(alpha * rhs.coeff(j,c));
102  for(LhsInnerIterator it(lhsEval,j); it ;++it)
103  res.coeffRef(it.index(),c) += it.value() * rhs_j;
104  }
105  }
106  }
107 };
108 
109 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
110 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, false>
111 {
116  typedef typename LhsEval::InnerIterator LhsInnerIterator;
117  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
118  {
119  Index n = lhs.rows();
120  LhsEval lhsEval(lhs);
121 
122 #ifdef EIGEN_HAS_OPENMP
124  Index threads = Eigen::nbThreads();
125  // This 20000 threshold has been found experimentally on 2D and 3D Poisson problems.
126  // It basically represents the minimal amount of work to be done to be worth it.
127  if(threads>1 && lhsEval.nonZerosEstimate()*rhs.cols() > 20000)
128  {
129  #pragma omp parallel for schedule(dynamic,(n+threads*4-1)/(threads*4)) num_threads(threads)
130  for(Index i=0; i<n; ++i)
131  processRow(lhsEval,rhs,res,alpha,i);
132  }
133  else
134 #endif
135  {
136  for(Index i=0; i<n; ++i)
137  processRow(lhsEval, rhs, res, alpha, i);
138  }
139  }
140 
141  static void processRow(const LhsEval& lhsEval, const DenseRhsType& rhs, Res& res, const typename Res::Scalar& alpha, Index i)
142  {
143  typename Res::RowXpr res_i(res.row(i));
144  for(LhsInnerIterator it(lhsEval,i); it ;++it)
145  res_i += (alpha*it.value()) * rhs.row(it.index());
146  }
147 };
148 
149 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
150 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, ColMajor, false>
151 {
156  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
157  {
158  evaluator<Lhs> lhsEval(lhs);
159  for(Index j=0; j<lhs.outerSize(); ++j)
160  {
161  typename Rhs::ConstRowXpr rhs_j(rhs.row(j));
162  for(LhsInnerIterator it(lhsEval,j); it ;++it)
163  res.row(it.index()) += (alpha*it.value()) * rhs_j;
164  }
165  }
166 };
167 
168 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
169 inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
170 {
172 }
173 
174 } // end namespace internal
175 
176 namespace internal {
177 
178 template<typename Lhs, typename Rhs, int ProductType>
180  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SparseShape,DenseShape,ProductType> >
181 {
183 
184  template<typename Dest>
185  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
186  {
189  LhsNested lhsNested(lhs);
190  RhsNested rhsNested(rhs);
191  internal::sparse_time_dense_product(lhsNested, rhsNested, dst, alpha);
192  }
193 };
194 
195 template<typename Lhs, typename Rhs, int ProductType>
197  : generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
198 {};
199 
200 template<typename Lhs, typename Rhs, int ProductType>
202  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SparseShape,ProductType> >
203 {
205 
206  template<typename Dst>
207  static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
208  {
211  LhsNested lhsNested(lhs);
212  RhsNested rhsNested(rhs);
213 
214  // transpose everything
215  Transpose<Dst> dstT(dst);
216  internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha);
217  }
218 };
219 
220 template<typename Lhs, typename Rhs, int ProductType>
222  : generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
223 {};
224 
225 template<typename LhsT, typename RhsT, bool NeedToTranspose>
227 {
228 protected:
232 
233  // if the actual left-hand side is a dense vector,
234  // then build a sparse-view so that we can seamlessly iterate over it.
237  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
239 
243  typedef typename ProdXprType::Scalar Scalar;
244 
245 public:
246  enum {
247  Flags = NeedToTranspose ? RowMajorBit : 0,
248  CoeffReadCost = HugeCost
249  };
250 
251  class InnerIterator : public LhsIterator
252  {
253  public:
255  : LhsIterator(xprEval.m_lhsXprImpl, 0),
256  m_outer(outer),
257  m_empty(false),
258  m_factor(get(xprEval.m_rhsXprImpl, outer, typename internal::traits<ActualRhs>::StorageKind() ))
259  {}
260 
261  EIGEN_STRONG_INLINE Index outer() const { return m_outer; }
262  EIGEN_STRONG_INLINE Index row() const { return NeedToTranspose ? m_outer : LhsIterator::index(); }
263  EIGEN_STRONG_INLINE Index col() const { return NeedToTranspose ? LhsIterator::index() : m_outer; }
264 
265  EIGEN_STRONG_INLINE Scalar value() const { return LhsIterator::value() * m_factor; }
266  EIGEN_STRONG_INLINE operator bool() const { return LhsIterator::operator bool() && (!m_empty); }
267 
268  protected:
269  Scalar get(const RhsEval &rhs, Index outer, Dense = Dense()) const
270  {
271  return rhs.coeff(outer);
272  }
273 
274  Scalar get(const RhsEval &rhs, Index outer, Sparse = Sparse())
275  {
276  typename RhsEval::InnerIterator it(rhs, outer);
277  if (it && it.index()==0 && it.value()!=Scalar(0))
278  return it.value();
279  m_empty = true;
280  return Scalar(0);
281  }
282 
284  bool m_empty;
285  Scalar m_factor;
286  };
287 
288  sparse_dense_outer_product_evaluator(const Lhs1 &lhs, const ActualRhs &rhs)
289  : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
290  {
291  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
292  }
293 
294  // transpose case
295  sparse_dense_outer_product_evaluator(const ActualRhs &rhs, const Lhs1 &lhs)
296  : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
297  {
298  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
299  }
300 
301 protected:
302  const LhsArg m_lhs;
305 };
306 
307 // sparse * dense outer product
308 template<typename Lhs, typename Rhs>
310  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor>
311 {
313 
315  typedef typename XprType::PlainObject PlainObject;
316 
317  explicit product_evaluator(const XprType& xpr)
318  : Base(xpr.lhs(), xpr.rhs())
319  {}
320 
321 };
322 
323 template<typename Lhs, typename Rhs>
325  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor>
326 {
328 
330  typedef typename XprType::PlainObject PlainObject;
331 
332  explicit product_evaluator(const XprType& xpr)
333  : Base(xpr.lhs(), xpr.rhs())
334  {}
335 
336 };
337 
338 } // end namespace internal
339 
340 } // end namespace Eigen
341 
342 #endif // EIGEN_SPARSEDENSEPRODUCT_H
Block< Derived, 1, internal::traits< Derived >::ColsAtCompileTime, IsRowMajor > RowXpr
Definition: BlockMethods.h:17
conditional< NeedToTranspose, RhsT, LhsT >::type Lhs1
static void processRow(const LhsEval &lhsEval, const DenseRhsType &rhs, DenseResType &res, const typename Res::Scalar &alpha, Index i, Index col)
void sparse_time_dense_product(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const AlphaType &alpha)
static void run(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const typename Res::Scalar &alpha)
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define EIGEN_STRONG_INLINE
Definition: Macros.h:917
void initParallel()
Definition: Parallelizer.h:53
const int HugeCost
Definition: Constants.h:44
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
static void run(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const typename Res::Scalar &alpha)
conditional< NeedToTranspose, LhsT, RhsT >::type ActualRhs
sparse_dense_outer_product_evaluator(const ActualRhs &rhs, const Lhs1 &lhs)
Expression of the transpose of a matrix.
Definition: Transpose.h:52
const StorageIndex & row() const
Definition: SparseUtil.h:172
int n
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Product< LhsT, RhsT, DefaultProduct > ProdXprType
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
static void run(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const typename Res::Scalar &alpha)
const unsigned int RowMajorBit
Definition: Constants.h:66
Eigen::SparseMatrix< double > Sparse
int nbThreads()
Definition: Parallelizer.h:63
static void scaleAndAddTo(Dst &dst, const Lhs &lhs, const Rhs &rhs, const Scalar &alpha)
static void processRow(const LhsEval &lhsEval, const DenseRhsType &rhs, Res &res, const typename Res::Scalar &alpha, Index i)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
#define EIGEN_INTERNAL_CHECK_COST_VALUE(C)
Definition: StaticAssert.h:218
static void scaleAndAddTo(Dest &dst, const Lhs &lhs, const Rhs &rhs, const Scalar &alpha)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
RealScalar alpha
evaluator< ActualLhs >::InnerIterator LhsIterator
Expression of a dense or sparse matrix with zero or too small values removed.
const Block< const Derived, 1, internal::traits< Derived >::ColsAtCompileTime, IsRowMajor > ConstRowXpr
Definition: BlockMethods.h:18
conditional< is_same< typename internal::traits< Lhs1 >::StorageKind, Sparse >::value, Lhs1 const &, SparseView< Lhs1 > >::type LhsArg
static void run(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const AlphaType &alpha)
sparse_dense_outer_product_evaluator(const Lhs1 &lhs, const ActualRhs &rhs)
m col(1)
Determines whether the given binary operation of two numeric types is allowed and what the scalar ret...
Definition: XprHelper.h:801
InnerIterator(const sparse_dense_outer_product_evaluator &xprEval, Index outer)
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
conditional< is_same< typename internal::traits< Lhs1 >::StorageKind, Sparse >::value, Lhs1, SparseView< Lhs1 > >::type ActualLhs
std::ptrdiff_t j
Definition: pytypes.h:1370


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:57