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  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
93  {
94  evaluator<Lhs> lhsEval(lhs);
95  for(Index c=0; c<rhs.cols(); ++c)
96  {
97  for(Index j=0; j<lhs.outerSize(); ++j)
98  {
99 // typename Res::Scalar rhs_j = alpha * rhs.coeff(j,c);
100  typename ScalarBinaryOpTraits<AlphaType, typename Rhs::Scalar>::ReturnType rhs_j(alpha * rhs.coeff(j,c));
101  for(LhsInnerIterator it(lhsEval,j); it ;++it)
102  res.coeffRef(it.index(),c) += it.value() * rhs_j;
103  }
104  }
105  }
106 };
107 
108 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
109 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, false>
110 {
115  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
116  {
117  evaluator<Lhs> lhsEval(lhs);
118  for(Index j=0; j<lhs.outerSize(); ++j)
119  {
120  typename Res::RowXpr res_j(res.row(j));
121  for(LhsInnerIterator it(lhsEval,j); it ;++it)
122  res_j += (alpha*it.value()) * rhs.row(it.index());
123  }
124  }
125 };
126 
127 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
128 struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, ColMajor, false>
129 {
134  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
135  {
136  evaluator<Lhs> lhsEval(lhs);
137  for(Index j=0; j<lhs.outerSize(); ++j)
138  {
139  typename Rhs::ConstRowXpr rhs_j(rhs.row(j));
140  for(LhsInnerIterator it(lhsEval,j); it ;++it)
141  res.row(it.index()) += (alpha*it.value()) * rhs_j;
142  }
143  }
144 };
145 
146 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
147 inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
148 {
150 }
151 
152 } // end namespace internal
153 
154 namespace internal {
155 
156 template<typename Lhs, typename Rhs, int ProductType>
158  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SparseShape,DenseShape,ProductType> >
159 {
161 
162  template<typename Dest>
163  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
164  {
167  LhsNested lhsNested(lhs);
168  RhsNested rhsNested(rhs);
169  internal::sparse_time_dense_product(lhsNested, rhsNested, dst, alpha);
170  }
171 };
172 
173 template<typename Lhs, typename Rhs, int ProductType>
175  : generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
176 {};
177 
178 template<typename Lhs, typename Rhs, int ProductType>
180  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SparseShape,ProductType> >
181 {
183 
184  template<typename Dst>
185  static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
186  {
189  LhsNested lhsNested(lhs);
190  RhsNested rhsNested(rhs);
191 
192  // transpose everything
193  Transpose<Dst> dstT(dst);
194  internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha);
195  }
196 };
197 
198 template<typename Lhs, typename Rhs, int ProductType>
200  : generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
201 {};
202 
203 template<typename LhsT, typename RhsT, bool NeedToTranspose>
205 {
206 protected:
210 
211  // if the actual left-hand side is a dense vector,
212  // then build a sparse-view so that we can seamlessly iterate over it.
215  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
217 
221  typedef typename ProdXprType::Scalar Scalar;
222 
223 public:
224  enum {
225  Flags = NeedToTranspose ? RowMajorBit : 0,
226  CoeffReadCost = HugeCost
227  };
228 
229  class InnerIterator : public LhsIterator
230  {
231  public:
233  : LhsIterator(xprEval.m_lhsXprImpl, 0),
234  m_outer(outer),
235  m_empty(false),
236  m_factor(get(xprEval.m_rhsXprImpl, outer, typename internal::traits<ActualRhs>::StorageKind() ))
237  {}
238 
239  EIGEN_STRONG_INLINE Index outer() const { return m_outer; }
240  EIGEN_STRONG_INLINE Index row() const { return NeedToTranspose ? m_outer : LhsIterator::index(); }
241  EIGEN_STRONG_INLINE Index col() const { return NeedToTranspose ? LhsIterator::index() : m_outer; }
242 
243  EIGEN_STRONG_INLINE Scalar value() const { return LhsIterator::value() * m_factor; }
244  EIGEN_STRONG_INLINE operator bool() const { return LhsIterator::operator bool() && (!m_empty); }
245 
246  protected:
247  Scalar get(const RhsEval &rhs, Index outer, Dense = Dense()) const
248  {
249  return rhs.coeff(outer);
250  }
251 
252  Scalar get(const RhsEval &rhs, Index outer, Sparse = Sparse())
253  {
254  typename RhsEval::InnerIterator it(rhs, outer);
255  if (it && it.index()==0 && it.value()!=Scalar(0))
256  return it.value();
257  m_empty = true;
258  return Scalar(0);
259  }
260 
262  bool m_empty;
263  Scalar m_factor;
264  };
265 
266  sparse_dense_outer_product_evaluator(const Lhs1 &lhs, const ActualRhs &rhs)
267  : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
268  {
269  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
270  }
271 
272  // transpose case
273  sparse_dense_outer_product_evaluator(const ActualRhs &rhs, const Lhs1 &lhs)
274  : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs)
275  {
276  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
277  }
278 
279 protected:
280  const LhsArg m_lhs;
283 };
284 
285 // sparse * dense outer product
286 template<typename Lhs, typename Rhs>
288  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor>
289 {
291 
293  typedef typename XprType::PlainObject PlainObject;
294 
295  explicit product_evaluator(const XprType& xpr)
296  : Base(xpr.lhs(), xpr.rhs())
297  {}
298 
299 };
300 
301 template<typename Lhs, typename Rhs>
303  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor>
304 {
306 
308  typedef typename XprType::PlainObject PlainObject;
309 
310  explicit product_evaluator(const XprType& xpr)
311  : Base(xpr.lhs(), xpr.rhs())
312  {}
313 
314 };
315 
316 } // end namespace internal
317 
318 } // end namespace Eigen
319 
320 #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:33
#define EIGEN_STRONG_INLINE
Definition: Macros.h:494
void initParallel()
Definition: Parallelizer.h:48
const int HugeCost
Definition: Constants.h:39
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
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:61
Eigen::SparseMatrix< double > Sparse
int nbThreads()
Definition: Parallelizer.h:58
static void scaleAndAddTo(Dst &dst, const Lhs &lhs, const Rhs &rhs, const Scalar &alpha)
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:215
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:33
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:766
InnerIterator(const sparse_dense_outer_product_evaluator &xprEval, Index outer)
conditional< is_same< typename internal::traits< Lhs1 >::StorageKind, Sparse >::value, Lhs1, SparseView< Lhs1 > >::type ActualLhs
void run(Expr &expr, Dev &dev)
Definition: TensorSyclRun.h:33
std::ptrdiff_t j
Definition: pytypes.h:897


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:21