Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
SparseDenseProduct.h
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>
24 struct sparse_time_dense_product_impl;
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 {
29  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
30  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
31  typedef typename internal::remove_all<DenseResType>::type Res;
32  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
33  typedef evaluator<Lhs> LhsEval;
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
40  Eigen::initParallel();
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 {
88  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
89  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
90  typedef typename internal::remove_all<DenseResType>::type Res;
91  typedef evaluator<Lhs> LhsEval;
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 {
112  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
113  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
114  typedef typename internal::remove_all<DenseResType>::type Res;
115  typedef evaluator<Lhs> LhsEval;
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
123  Eigen::initParallel();
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 {
152  typedef typename internal::remove_all<SparseLhsType>::type Lhs;
153  typedef typename internal::remove_all<DenseRhsType>::type Rhs;
154  typedef typename internal::remove_all<DenseResType>::type Res;
155  typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
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 {
171  sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType>::run(lhs, rhs, res, alpha);
172 }
173 
174 } // end namespace internal
175 
176 namespace internal {
177 
178 template<typename Lhs, typename Rhs, int ProductType>
179 struct generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
180  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SparseShape,DenseShape,ProductType> >
181 {
182  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
183 
184  template<typename Dest>
185  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
186  {
187  typedef typename nested_eval<Lhs,((Rhs::Flags&RowMajorBit)==0) ? 1 : Rhs::ColsAtCompileTime>::type LhsNested;
188  typedef typename nested_eval<Rhs,((Lhs::Flags&RowMajorBit)==0) ? 1 : Dynamic>::type RhsNested;
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>
196 struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, DenseShape, ProductType>
197  : generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
198 {};
199 
200 template<typename Lhs, typename Rhs, int ProductType>
201 struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
202  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SparseShape,ProductType> >
203 {
204  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
205 
206  template<typename Dst>
207  static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
208  {
209  typedef typename nested_eval<Lhs,((Rhs::Flags&RowMajorBit)==0) ? Dynamic : 1>::type LhsNested;
210  typedef typename nested_eval<Rhs,((Lhs::Flags&RowMajorBit)==RowMajorBit) ? 1 : Lhs::RowsAtCompileTime>::type RhsNested;
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>
221 struct generic_product_impl<Lhs, Rhs, DenseShape, SparseTriangularShape, ProductType>
222  : generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
223 {};
224 
225 template<typename LhsT, typename RhsT, bool NeedToTranspose>
226 struct sparse_dense_outer_product_evaluator
227 {
228 protected:
229  typedef typename conditional<NeedToTranspose,RhsT,LhsT>::type Lhs1;
230  typedef typename conditional<NeedToTranspose,LhsT,RhsT>::type ActualRhs;
231  typedef Product<LhsT,RhsT,DefaultProduct> ProdXprType;
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.
235  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
236  Lhs1, SparseView<Lhs1> >::type ActualLhs;
237  typedef typename conditional<is_same<typename internal::traits<Lhs1>::StorageKind,Sparse>::value,
238  Lhs1 const&, SparseView<Lhs1> >::type LhsArg;
239 
240  typedef evaluator<ActualLhs> LhsEval;
241  typedef evaluator<ActualRhs> RhsEval;
242  typedef typename evaluator<ActualLhs>::InnerIterator LhsIterator;
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:
254  InnerIterator(const sparse_dense_outer_product_evaluator &xprEval, Index outer)
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 
283  Index m_outer;
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;
303  evaluator<ActualLhs> m_lhsXprImpl;
304  evaluator<ActualRhs> m_rhsXprImpl;
305 };
306 
307 // sparse * dense outer product
308 template<typename Lhs, typename Rhs>
309 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, OuterProduct, SparseShape, DenseShape>
310  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor>
311 {
312  typedef sparse_dense_outer_product_evaluator<Lhs,Rhs, Lhs::IsRowMajor> Base;
313 
314  typedef Product<Lhs, Rhs> XprType;
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>
324 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, OuterProduct, DenseShape, SparseShape>
325  : sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor>
326 {
327  typedef sparse_dense_outer_product_evaluator<Lhs,Rhs, Rhs::IsRowMajor> Base;
328 
329  typedef Product<Lhs, Rhs> XprType;
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
Definition: Constants.h:319
const int HugeCost
Definition: Constants.h:44
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const unsigned int RowMajorBit
Definition: Constants.h:66
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Eigen_Colamd.h:50
Definition: Constants.h:321