Eigen  3.2.93
ProductEvaluators.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
6 // Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 
13 #ifndef EIGEN_PRODUCTEVALUATORS_H
14 #define EIGEN_PRODUCTEVALUATORS_H
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
28 template<typename Lhs, typename Rhs, int Options>
29 struct evaluator<Product<Lhs, Rhs, Options> >
30  : public product_evaluator<Product<Lhs, Rhs, Options> >
31 {
32  typedef Product<Lhs, Rhs, Options> XprType;
33  typedef product_evaluator<XprType> Base;
34 
35  EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr) : Base(xpr) {}
36 };
37 
38 // Catch "scalar * ( A * B )" and transform it to "(A*scalar) * B"
39 // TODO we should apply that rule only if that's really helpful
40 template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
41 struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
42  const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
43  const Product<Lhs, Rhs, DefaultProduct> > >
44 {
45  static const bool value = true;
46 };
47 template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
48 struct evaluator<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
49  const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
50  const Product<Lhs, Rhs, DefaultProduct> > >
51  : public evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> >
52 {
53  typedef CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
54  const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
55  const Product<Lhs, Rhs, DefaultProduct> > XprType;
56  typedef evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> > Base;
57 
58  EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr)
59  : Base(xpr.lhs().functor().m_other * xpr.rhs().lhs() * xpr.rhs().rhs())
60  {}
61 };
62 
63 
64 template<typename Lhs, typename Rhs, int DiagIndex>
65 struct evaluator<Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> >
66  : public evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> >
67 {
68  typedef Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> XprType;
69  typedef evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> > Base;
70 
71  EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr)
72  : Base(Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex>(
73  Product<Lhs, Rhs, LazyProduct>(xpr.nestedExpression().lhs(), xpr.nestedExpression().rhs()),
74  xpr.index() ))
75  {}
76 };
77 
78 
79 // Helper class to perform a matrix product with the destination at hand.
80 // Depending on the sizes of the factors, there are different evaluation strategies
81 // as controlled by internal::product_type.
82 template< typename Lhs, typename Rhs,
83  typename LhsShape = typename evaluator_traits<Lhs>::Shape,
84  typename RhsShape = typename evaluator_traits<Rhs>::Shape,
85  int ProductType = internal::product_type<Lhs,Rhs>::value>
86 struct generic_product_impl;
87 
88 template<typename Lhs, typename Rhs>
89 struct evaluator_assume_aliasing<Product<Lhs, Rhs, DefaultProduct> > {
90  static const bool value = true;
91 };
92 
93 // This is the default evaluator implementation for products:
94 // It creates a temporary and call generic_product_impl
95 template<typename Lhs, typename Rhs, int Options, int ProductTag, typename LhsShape, typename RhsShape>
96 struct product_evaluator<Product<Lhs, Rhs, Options>, ProductTag, LhsShape, RhsShape>
97  : public evaluator<typename Product<Lhs, Rhs, Options>::PlainObject>
98 {
99  typedef Product<Lhs, Rhs, Options> XprType;
100  typedef typename XprType::PlainObject PlainObject;
101  typedef evaluator<PlainObject> Base;
102  enum {
103  Flags = Base::Flags | EvalBeforeNestingBit
104  };
105 
106  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
107  explicit product_evaluator(const XprType& xpr)
108  : m_result(xpr.rows(), xpr.cols())
109  {
110  ::new (static_cast<Base*>(this)) Base(m_result);
111 
112 // FIXME shall we handle nested_eval here?,
113 // if so, then we must take care at removing the call to nested_eval in the specializations (e.g., in permutation_matrix_product, transposition_matrix_product, etc.)
114 // typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
115 // typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
116 // typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
117 // typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
118 //
119 // const LhsNested lhs(xpr.lhs());
120 // const RhsNested rhs(xpr.rhs());
121 //
122 // generic_product_impl<LhsNestedCleaned, RhsNestedCleaned>::evalTo(m_result, lhs, rhs);
123 
124  generic_product_impl<Lhs, Rhs, LhsShape, RhsShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
125  }
126 
127 protected:
128  PlainObject m_result;
129 };
130 
131 // The following three shortcuts are enabled only if the scalar types match excatly.
132 // TODO: we could enable them for different scalar types when the product is not vectorized.
133 
134 // Dense = Product
135 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
136 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::assign_op<Scalar,Scalar>, Dense2Dense,
137  typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
138 {
139  typedef Product<Lhs,Rhs,Options> SrcXprType;
140  static EIGEN_STRONG_INLINE
141  void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
142  {
143  // FIXME shall we handle nested_eval here?
144  generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
145  }
146 };
147 
148 // Dense += Product
149 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
150 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::add_assign_op<Scalar,Scalar>, Dense2Dense,
151  typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
152 {
153  typedef Product<Lhs,Rhs,Options> SrcXprType;
154  static EIGEN_STRONG_INLINE
155  void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,Scalar> &)
156  {
157  // FIXME shall we handle nested_eval here?
158  generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs());
159  }
160 };
161 
162 // Dense -= Product
163 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
164 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::sub_assign_op<Scalar,Scalar>, Dense2Dense,
165  typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
166 {
167  typedef Product<Lhs,Rhs,Options> SrcXprType;
168  static EIGEN_STRONG_INLINE
169  void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,Scalar> &)
170  {
171  // FIXME shall we handle nested_eval here?
172  generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs());
173  }
174 };
175 
176 
177 // Dense ?= scalar * Product
178 // TODO we should apply that rule if that's really helpful
179 // for instance, this is not good for inner products
180 template< typename DstXprType, typename Lhs, typename Rhs, typename AssignFunc, typename Scalar, typename ScalarBis, typename Plain>
181 struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>, const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
182  const Product<Lhs,Rhs,DefaultProduct> >, AssignFunc, Dense2Dense>
183 {
184  typedef CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>,
185  const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
186  const Product<Lhs,Rhs,DefaultProduct> > SrcXprType;
187  static EIGEN_STRONG_INLINE
188  void run(DstXprType &dst, const SrcXprType &src, const AssignFunc& func)
189  {
190  call_assignment_no_alias(dst, (src.lhs().functor().m_other * src.rhs().lhs())*src.rhs().rhs(), func);
191  }
192 };
193 
194 //----------------------------------------
195 // Catch "Dense ?= xpr + Product<>" expression to save one temporary
196 // FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct
197 // TODO enable it for "Dense ?= xpr - Product<>" as well.
198 
199 template<typename OtherXpr, typename Lhs, typename Rhs>
200 struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
201  const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
202  static const bool value = true;
203 };
204 
205 template<typename DstXprType, typename OtherXpr, typename ProductType, typename Func1, typename Func2>
206 struct assignment_from_xpr_plus_product
207 {
208  typedef CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar,typename ProductType::Scalar>, const OtherXpr, const ProductType> SrcXprType;
209  template<typename InitialFunc>
210  static EIGEN_STRONG_INLINE
211  void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
212  {
213  call_assignment_no_alias(dst, src.lhs(), Func1());
214  call_assignment_no_alias(dst, src.rhs(), Func2());
215  }
216 };
217 
218 template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar>
219 struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_sum_op<OtherScalar,ProdScalar>, const OtherXpr,
220  const Product<Lhs,Rhs,DefaultProduct> >, internal::assign_op<DstScalar,SrcScalar>, Dense2Dense>
221  : assignment_from_xpr_plus_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<DstScalar,OtherScalar>, internal::add_assign_op<DstScalar,ProdScalar> >
222 {};
223 template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar>
224 struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_sum_op<OtherScalar,ProdScalar>, const OtherXpr,
225  const Product<Lhs,Rhs,DefaultProduct> >, internal::add_assign_op<DstScalar,SrcScalar>, Dense2Dense>
226  : assignment_from_xpr_plus_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<DstScalar,OtherScalar>, internal::add_assign_op<DstScalar,ProdScalar> >
227 {};
228 template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar>
229 struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_sum_op<OtherScalar,ProdScalar>, const OtherXpr,
230  const Product<Lhs,Rhs,DefaultProduct> >, internal::sub_assign_op<DstScalar,SrcScalar>, Dense2Dense>
231  : assignment_from_xpr_plus_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<DstScalar,OtherScalar>, internal::sub_assign_op<DstScalar,ProdScalar> >
232 {};
233 //----------------------------------------
234 
235 template<typename Lhs, typename Rhs>
236 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,InnerProduct>
237 {
238  template<typename Dst>
239  static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
240  {
241  dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
242  }
243 
244  template<typename Dst>
245  static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
246  {
247  dst.coeffRef(0,0) += (lhs.transpose().cwiseProduct(rhs)).sum();
248  }
249 
250  template<typename Dst>
251  static void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
252  { dst.coeffRef(0,0) -= (lhs.transpose().cwiseProduct(rhs)).sum(); }
253 };
254 
255 
256 /***********************************************************************
257 * Implementation of outer dense * dense vector product
258 ***********************************************************************/
259 
260 // Column major result
261 template<typename Dst, typename Lhs, typename Rhs, typename Func>
262 void outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const false_type&)
263 {
264  evaluator<Rhs> rhsEval(rhs);
265  typename nested_eval<Lhs,Rhs::SizeAtCompileTime>::type actual_lhs(lhs);
266  // FIXME if cols is large enough, then it might be useful to make sure that lhs is sequentially stored
267  // FIXME not very good if rhs is real and lhs complex while alpha is real too
268  const Index cols = dst.cols();
269  for (Index j=0; j<cols; ++j)
270  func(dst.col(j), rhsEval.coeff(0,j) * actual_lhs);
271 }
272 
273 // Row major result
274 template<typename Dst, typename Lhs, typename Rhs, typename Func>
275 void outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const true_type&)
276 {
277  evaluator<Lhs> lhsEval(lhs);
278  typename nested_eval<Rhs,Lhs::SizeAtCompileTime>::type actual_rhs(rhs);
279  // FIXME if rows is large enough, then it might be useful to make sure that rhs is sequentially stored
280  // FIXME not very good if lhs is real and rhs complex while alpha is real too
281  const Index rows = dst.rows();
282  for (Index i=0; i<rows; ++i)
283  func(dst.row(i), lhsEval.coeff(i,0) * actual_rhs);
284 }
285 
286 template<typename Lhs, typename Rhs>
287 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,OuterProduct>
288 {
289  template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
290  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
291 
292  // TODO it would be nice to be able to exploit our *_assign_op functors for that purpose
293  struct set { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } };
294  struct add { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
295  struct sub { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
296  struct adds {
297  Scalar m_scale;
298  explicit adds(const Scalar& s) : m_scale(s) {}
299  template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const {
300  dst.const_cast_derived() += m_scale * src;
301  }
302  };
303 
304  template<typename Dst>
305  static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
306  {
307  internal::outer_product_selector_run(dst, lhs, rhs, set(), is_row_major<Dst>());
308  }
309 
310  template<typename Dst>
311  static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
312  {
313  internal::outer_product_selector_run(dst, lhs, rhs, add(), is_row_major<Dst>());
314  }
315 
316  template<typename Dst>
317  static inline void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
318  {
319  internal::outer_product_selector_run(dst, lhs, rhs, sub(), is_row_major<Dst>());
320  }
321 
322  template<typename Dst>
323  static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
324  {
325  internal::outer_product_selector_run(dst, lhs, rhs, adds(alpha), is_row_major<Dst>());
326  }
327 
328 };
329 
330 
331 // This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo
332 template<typename Lhs, typename Rhs, typename Derived>
333 struct generic_product_impl_base
334 {
335  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
336 
337  template<typename Dst>
338  static EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
339  { dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
340 
341  template<typename Dst>
342  static EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
343  { scaleAndAddTo(dst,lhs, rhs, Scalar(1)); }
344 
345  template<typename Dst>
346  static EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
347  { scaleAndAddTo(dst, lhs, rhs, Scalar(-1)); }
348 
349  template<typename Dst>
350  static EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
351  { Derived::scaleAndAddTo(dst,lhs,rhs,alpha); }
352 
353 };
354 
355 template<typename Lhs, typename Rhs>
356 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct>
357  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct> >
358 {
359  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
360  enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
361  typedef typename internal::conditional<int(Side)==OnTheRight,Lhs,Rhs>::type MatrixType;
362 
363  template<typename Dest>
364  static EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
365  {
366  internal::gemv_dense_selector<Side,
367  (int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
368  bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)
369  >::run(lhs, rhs, dst, alpha);
370  }
371 };
372 
373 template<typename Lhs, typename Rhs>
374 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode>
375 {
376  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
377 
378  template<typename Dst>
379  static EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
380  {
381  // Same as: dst.noalias() = lhs.lazyProduct(rhs);
382  // but easier on the compiler side
383  call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op<typename Dst::Scalar,Scalar>());
384  }
385 
386  template<typename Dst>
387  static EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
388  {
389  // dst.noalias() += lhs.lazyProduct(rhs);
390  call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op<typename Dst::Scalar,Scalar>());
391  }
392 
393  template<typename Dst>
394  static EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
395  {
396  // dst.noalias() -= lhs.lazyProduct(rhs);
397  call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op<typename Dst::Scalar,Scalar>());
398  }
399 
400 // template<typename Dst>
401 // static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
402 // { dst.noalias() += alpha * lhs.lazyProduct(rhs); }
403 };
404 
405 // This specialization enforces the use of a coefficient-based evaluation strategy
406 template<typename Lhs, typename Rhs>
407 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,LazyCoeffBasedProductMode>
408  : generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode> {};
409 
410 // Case 2: Evaluate coeff by coeff
411 //
412 // This is mostly taken from CoeffBasedProduct.h
413 // The main difference is that we add an extra argument to the etor_product_*_impl::run() function
414 // for the inner dimension of the product, because evaluator object do not know their size.
415 
416 template<int Traversal, int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar>
417 struct etor_product_coeff_impl;
418 
419 template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
420 struct etor_product_packet_impl;
421 
422 template<typename Lhs, typename Rhs, int ProductTag>
423 struct product_evaluator<Product<Lhs, Rhs, LazyProduct>, ProductTag, DenseShape, DenseShape>
424  : evaluator_base<Product<Lhs, Rhs, LazyProduct> >
425 {
426  typedef Product<Lhs, Rhs, LazyProduct> XprType;
427  typedef typename XprType::Scalar Scalar;
428  typedef typename XprType::CoeffReturnType CoeffReturnType;
429 
430  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
431  explicit product_evaluator(const XprType& xpr)
432  : m_lhs(xpr.lhs()),
433  m_rhs(xpr.rhs()),
434  m_lhsImpl(m_lhs), // FIXME the creation of the evaluator objects should result in a no-op, but check that!
435  m_rhsImpl(m_rhs), // Moreover, they are only useful for the packet path, so we could completely disable them when not needed,
436  // or perhaps declare them on the fly on the packet method... We have experiment to check what's best.
437  m_innerDim(xpr.lhs().cols())
438  {
439  EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
440  EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::AddCost);
441  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
442  }
443 
444  // Everything below here is taken from CoeffBasedProduct.h
445 
446  typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
447  typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
448 
449  typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
450  typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
451 
452  typedef evaluator<LhsNestedCleaned> LhsEtorType;
453  typedef evaluator<RhsNestedCleaned> RhsEtorType;
454 
455  enum {
456  RowsAtCompileTime = LhsNestedCleaned::RowsAtCompileTime,
457  ColsAtCompileTime = RhsNestedCleaned::ColsAtCompileTime,
458  InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsNestedCleaned::ColsAtCompileTime, RhsNestedCleaned::RowsAtCompileTime),
459  MaxRowsAtCompileTime = LhsNestedCleaned::MaxRowsAtCompileTime,
460  MaxColsAtCompileTime = RhsNestedCleaned::MaxColsAtCompileTime
461  };
462 
463  typedef typename find_best_packet<Scalar,RowsAtCompileTime>::type LhsVecPacketType;
464  typedef typename find_best_packet<Scalar,ColsAtCompileTime>::type RhsVecPacketType;
465 
466  enum {
467 
468  LhsCoeffReadCost = LhsEtorType::CoeffReadCost,
469  RhsCoeffReadCost = RhsEtorType::CoeffReadCost,
470  CoeffReadCost = InnerSize==0 ? NumTraits<Scalar>::ReadCost
471  : InnerSize == Dynamic ? HugeCost
472  : InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
473  + (InnerSize - 1) * NumTraits<Scalar>::AddCost,
474 
475  Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
476 
477  LhsFlags = LhsEtorType::Flags,
478  RhsFlags = RhsEtorType::Flags,
479 
480  LhsRowMajor = LhsFlags & RowMajorBit,
481  RhsRowMajor = RhsFlags & RowMajorBit,
482 
483  LhsVecPacketSize = unpacket_traits<LhsVecPacketType>::size,
484  RhsVecPacketSize = unpacket_traits<RhsVecPacketType>::size,
485 
486  // Here, we don't care about alignment larger than the usable packet size.
487  LhsAlignment = EIGEN_PLAIN_ENUM_MIN(LhsEtorType::Alignment,LhsVecPacketSize*int(sizeof(typename LhsNestedCleaned::Scalar))),
488  RhsAlignment = EIGEN_PLAIN_ENUM_MIN(RhsEtorType::Alignment,RhsVecPacketSize*int(sizeof(typename RhsNestedCleaned::Scalar))),
489 
490  SameType = is_same<typename LhsNestedCleaned::Scalar,typename RhsNestedCleaned::Scalar>::value,
491 
492  CanVectorizeRhs = bool(RhsRowMajor) && (RhsFlags & PacketAccessBit),
493  CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit),
494 
495  EvalToRowMajor = (MaxRowsAtCompileTime==1&&MaxColsAtCompileTime!=1) ? 1
496  : (MaxColsAtCompileTime==1&&MaxRowsAtCompileTime!=1) ? 0
497  : (bool(RhsRowMajor) && !CanVectorizeLhs),
498 
499  Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & ~RowMajorBit)
500  | (EvalToRowMajor ? RowMajorBit : 0)
501  // TODO enable vectorization for mixed types
502  | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0)
503  | (XprType::IsVectorAtCompileTime ? LinearAccessBit : 0),
504 
505  LhsOuterStrideBytes = int(LhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename LhsNestedCleaned::Scalar)),
506  RhsOuterStrideBytes = int(RhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename RhsNestedCleaned::Scalar)),
507 
508  Alignment = bool(CanVectorizeLhs) ? (LhsOuterStrideBytes<0 || (int(LhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,LhsAlignment))!=0 ? 0 : LhsAlignment)
509  : bool(CanVectorizeRhs) ? (RhsOuterStrideBytes<0 || (int(RhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,RhsAlignment))!=0 ? 0 : RhsAlignment)
510  : 0,
511 
512  /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside
513  * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner
514  * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
515  * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
516  */
517  CanVectorizeInner = SameType
518  && LhsRowMajor
519  && (!RhsRowMajor)
520  && (LhsFlags & RhsFlags & ActualPacketAccessBit)
521  && (InnerSize % packet_traits<Scalar>::size == 0)
522  };
523 
524  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index row, Index col) const
525  {
526  return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
527  }
528 
529  /* Allow index-based non-packet access. It is impossible though to allow index-based packed access,
530  * which is why we don't set the LinearAccessBit.
531  * TODO: this seems possible when the result is a vector
532  */
533  EIGEN_DEVICE_FUNC const CoeffReturnType coeff(Index index) const
534  {
535  const Index row = RowsAtCompileTime == 1 ? 0 : index;
536  const Index col = RowsAtCompileTime == 1 ? index : 0;
537  return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
538  }
539 
540  template<int LoadMode, typename PacketType>
541  const PacketType packet(Index row, Index col) const
542  {
543  PacketType res;
544  typedef etor_product_packet_impl<bool(int(Flags)&RowMajorBit) ? RowMajor : ColMajor,
545  Unroll ? int(InnerSize) : Dynamic,
546  LhsEtorType, RhsEtorType, PacketType, LoadMode> PacketImpl;
547  PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res);
548  return res;
549  }
550 
551  template<int LoadMode, typename PacketType>
552  const PacketType packet(Index index) const
553  {
554  const Index row = RowsAtCompileTime == 1 ? 0 : index;
555  const Index col = RowsAtCompileTime == 1 ? index : 0;
556  return packet<LoadMode,PacketType>(row,col);
557  }
558 
559 protected:
560  const LhsNested m_lhs;
561  const RhsNested m_rhs;
562 
563  LhsEtorType m_lhsImpl;
564  RhsEtorType m_rhsImpl;
565 
566  // TODO: Get rid of m_innerDim if known at compile time
567  Index m_innerDim;
568 };
569 
570 template<typename Lhs, typename Rhs>
571 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, LazyCoeffBasedProductMode, DenseShape, DenseShape>
572  : product_evaluator<Product<Lhs, Rhs, LazyProduct>, CoeffBasedProductMode, DenseShape, DenseShape>
573 {
574  typedef Product<Lhs, Rhs, DefaultProduct> XprType;
575  typedef Product<Lhs, Rhs, LazyProduct> BaseProduct;
576  typedef product_evaluator<BaseProduct, CoeffBasedProductMode, DenseShape, DenseShape> Base;
577  enum {
578  Flags = Base::Flags | EvalBeforeNestingBit
579  };
580  EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
581  : Base(BaseProduct(xpr.lhs(),xpr.rhs()))
582  {}
583 };
584 
585 /****************************************
586 *** Coeff based product, Packet path ***
587 ****************************************/
588 
589 template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
590 struct etor_product_packet_impl<RowMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
591 {
592  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
593  {
594  etor_product_packet_impl<RowMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
595  res = pmadd(pset1<Packet>(lhs.coeff(row, UnrollingIndex-1)), rhs.template packet<LoadMode,Packet>(UnrollingIndex-1, col), res);
596  }
597 };
598 
599 template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
600 struct etor_product_packet_impl<ColMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
601 {
602  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
603  {
604  etor_product_packet_impl<ColMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
605  res = pmadd(lhs.template packet<LoadMode,Packet>(row, UnrollingIndex-1), pset1<Packet>(rhs.coeff(UnrollingIndex-1, col)), res);
606  }
607 };
608 
609 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
610 struct etor_product_packet_impl<RowMajor, 1, Lhs, Rhs, Packet, LoadMode>
611 {
612  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
613  {
614  res = pmul(pset1<Packet>(lhs.coeff(row, 0)),rhs.template packet<LoadMode,Packet>(0, col));
615  }
616 };
617 
618 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
619 struct etor_product_packet_impl<ColMajor, 1, Lhs, Rhs, Packet, LoadMode>
620 {
621  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
622  {
623  res = pmul(lhs.template packet<LoadMode,Packet>(row, 0), pset1<Packet>(rhs.coeff(0, col)));
624  }
625 };
626 
627 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
628 struct etor_product_packet_impl<RowMajor, 0, Lhs, Rhs, Packet, LoadMode>
629 {
630  static EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
631  {
632  res = pset1<Packet>(0);
633  }
634 };
635 
636 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
637 struct etor_product_packet_impl<ColMajor, 0, Lhs, Rhs, Packet, LoadMode>
638 {
639  static EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
640  {
641  res = pset1<Packet>(0);
642  }
643 };
644 
645 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
646 struct etor_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
647 {
648  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
649  {
650  res = pset1<Packet>(0);
651  for(Index i = 0; i < innerDim; ++i)
652  res = pmadd(pset1<Packet>(lhs.coeff(row, i)), rhs.template packet<LoadMode,Packet>(i, col), res);
653  }
654 };
655 
656 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
657 struct etor_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
658 {
659  static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
660  {
661  res = pset1<Packet>(0);
662  for(Index i = 0; i < innerDim; ++i)
663  res = pmadd(lhs.template packet<LoadMode,Packet>(row, i), pset1<Packet>(rhs.coeff(i, col)), res);
664  }
665 };
666 
667 
668 /***************************************************************************
669 * Triangular products
670 ***************************************************************************/
671 template<int Mode, bool LhsIsTriangular,
672  typename Lhs, bool LhsIsVector,
673  typename Rhs, bool RhsIsVector>
674 struct triangular_product_impl;
675 
676 template<typename Lhs, typename Rhs, int ProductTag>
677 struct generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag>
678  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag> >
679 {
680  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
681 
682  template<typename Dest>
683  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
684  {
685  triangular_product_impl<Lhs::Mode,true,typename Lhs::MatrixType,false,Rhs, Rhs::ColsAtCompileTime==1>
686  ::run(dst, lhs.nestedExpression(), rhs, alpha);
687  }
688 };
689 
690 template<typename Lhs, typename Rhs, int ProductTag>
691 struct generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag>
692 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag> >
693 {
694  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
695 
696  template<typename Dest>
697  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
698  {
699  triangular_product_impl<Rhs::Mode,false,Lhs,Lhs::RowsAtCompileTime==1, typename Rhs::MatrixType, false>::run(dst, lhs, rhs.nestedExpression(), alpha);
700  }
701 };
702 
703 
704 /***************************************************************************
705 * SelfAdjoint products
706 ***************************************************************************/
707 template <typename Lhs, int LhsMode, bool LhsIsVector,
708  typename Rhs, int RhsMode, bool RhsIsVector>
709 struct selfadjoint_product_impl;
710 
711 template<typename Lhs, typename Rhs, int ProductTag>
712 struct generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag>
713  : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag> >
714 {
715  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
716 
717  template<typename Dest>
718  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
719  {
720  selfadjoint_product_impl<typename Lhs::MatrixType,Lhs::Mode,false,Rhs,0,Rhs::IsVectorAtCompileTime>::run(dst, lhs.nestedExpression(), rhs, alpha);
721  }
722 };
723 
724 template<typename Lhs, typename Rhs, int ProductTag>
725 struct generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag>
726 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag> >
727 {
728  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
729 
730  template<typename Dest>
731  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
732  {
733  selfadjoint_product_impl<Lhs,0,Lhs::IsVectorAtCompileTime,typename Rhs::MatrixType,Rhs::Mode,false>::run(dst, lhs, rhs.nestedExpression(), alpha);
734  }
735 };
736 
737 
738 /***************************************************************************
739 * Diagonal products
740 ***************************************************************************/
741 
742 template<typename MatrixType, typename DiagonalType, typename Derived, int ProductOrder>
743 struct diagonal_product_evaluator_base
744  : evaluator_base<Derived>
745 {
746  typedef typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar;
747 public:
748  enum {
749  CoeffReadCost = NumTraits<Scalar>::MulCost + evaluator<MatrixType>::CoeffReadCost + evaluator<DiagonalType>::CoeffReadCost,
750 
751  MatrixFlags = evaluator<MatrixType>::Flags,
752  DiagFlags = evaluator<DiagonalType>::Flags,
753  _StorageOrder = MatrixFlags & RowMajorBit ? RowMajor : ColMajor,
754  _ScalarAccessOnDiag = !((int(_StorageOrder) == ColMajor && int(ProductOrder) == OnTheLeft)
755  ||(int(_StorageOrder) == RowMajor && int(ProductOrder) == OnTheRight)),
756  _SameTypes = is_same<typename MatrixType::Scalar, typename DiagonalType::Scalar>::value,
757  // FIXME currently we need same types, but in the future the next rule should be the one
758  //_Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && ((!_PacketOnDiag) || (_SameTypes && bool(int(DiagFlags)&PacketAccessBit))),
759  _Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && _SameTypes && (_ScalarAccessOnDiag || (bool(int(DiagFlags)&PacketAccessBit))),
760  _LinearAccessMask = (MatrixType::RowsAtCompileTime==1 || MatrixType::ColsAtCompileTime==1) ? LinearAccessBit : 0,
761  Flags = ((HereditaryBits|_LinearAccessMask) & (unsigned int)(MatrixFlags)) | (_Vectorizable ? PacketAccessBit : 0),
762  Alignment = evaluator<MatrixType>::Alignment
763  };
764 
765  diagonal_product_evaluator_base(const MatrixType &mat, const DiagonalType &diag)
766  : m_diagImpl(diag), m_matImpl(mat)
767  {
768  EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
769  EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
770  }
771 
772  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index idx) const
773  {
774  return m_diagImpl.coeff(idx) * m_matImpl.coeff(idx);
775  }
776 
777 protected:
778  template<int LoadMode,typename PacketType>
779  EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::true_type) const
780  {
781  return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
782  internal::pset1<PacketType>(m_diagImpl.coeff(id)));
783  }
784 
785  template<int LoadMode,typename PacketType>
786  EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::false_type) const
787  {
788  enum {
789  InnerSize = (MatrixType::Flags & RowMajorBit) ? MatrixType::ColsAtCompileTime : MatrixType::RowsAtCompileTime,
790  DiagonalPacketLoadMode = EIGEN_PLAIN_ENUM_MIN(LoadMode,((InnerSize%16) == 0) ? int(Aligned16) : int(evaluator<DiagonalType>::Alignment)) // FIXME hardcoded 16!!
791  };
792  return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
793  m_diagImpl.template packet<DiagonalPacketLoadMode,PacketType>(id));
794  }
795 
796  evaluator<DiagonalType> m_diagImpl;
797  evaluator<MatrixType> m_matImpl;
798 };
799 
800 // diagonal * dense
801 template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
802 struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, DenseShape>
803  : diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft>
804 {
805  typedef diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft> Base;
806  using Base::m_diagImpl;
807  using Base::m_matImpl;
808  using Base::coeff;
809  typedef typename Base::Scalar Scalar;
810 
811  typedef Product<Lhs, Rhs, ProductKind> XprType;
812  typedef typename XprType::PlainObject PlainObject;
813 
814  enum {
815  StorageOrder = int(Rhs::Flags) & RowMajorBit ? RowMajor : ColMajor
816  };
817 
818  EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
819  : Base(xpr.rhs(), xpr.lhs().diagonal())
820  {
821  }
822 
823  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
824  {
825  return m_diagImpl.coeff(row) * m_matImpl.coeff(row, col);
826  }
827 
828 #ifndef __CUDACC__
829  template<int LoadMode,typename PacketType>
830  EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
831  {
832  // FIXME: NVCC used to complain about the template keyword, but we have to check whether this is still the case.
833  // See also similar calls below.
834  return this->template packet_impl<LoadMode,PacketType>(row,col, row,
835  typename internal::conditional<int(StorageOrder)==RowMajor, internal::true_type, internal::false_type>::type());
836  }
837 
838  template<int LoadMode,typename PacketType>
839  EIGEN_STRONG_INLINE PacketType packet(Index idx) const
840  {
841  return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
842  }
843 #endif
844 };
845 
846 // dense * diagonal
847 template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
848 struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DenseShape, DiagonalShape>
849  : diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight>
850 {
851  typedef diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight> Base;
852  using Base::m_diagImpl;
853  using Base::m_matImpl;
854  using Base::coeff;
855  typedef typename Base::Scalar Scalar;
856 
857  typedef Product<Lhs, Rhs, ProductKind> XprType;
858  typedef typename XprType::PlainObject PlainObject;
859 
860  enum { StorageOrder = int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor };
861 
862  EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
863  : Base(xpr.lhs(), xpr.rhs().diagonal())
864  {
865  }
866 
867  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
868  {
869  return m_matImpl.coeff(row, col) * m_diagImpl.coeff(col);
870  }
871 
872 #ifndef __CUDACC__
873  template<int LoadMode,typename PacketType>
874  EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
875  {
876  return this->template packet_impl<LoadMode,PacketType>(row,col, col,
877  typename internal::conditional<int(StorageOrder)==ColMajor, internal::true_type, internal::false_type>::type());
878  }
879 
880  template<int LoadMode,typename PacketType>
881  EIGEN_STRONG_INLINE PacketType packet(Index idx) const
882  {
883  return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
884  }
885 #endif
886 };
887 
888 /***************************************************************************
889 * Products with permutation matrices
890 ***************************************************************************/
891 
897 template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
898 struct permutation_matrix_product;
899 
900 template<typename ExpressionType, int Side, bool Transposed>
901 struct permutation_matrix_product<ExpressionType, Side, Transposed, DenseShape>
902 {
903  typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
904  typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
905 
906  template<typename Dest, typename PermutationType>
907  static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
908  {
909  MatrixType mat(xpr);
910  const Index n = Side==OnTheLeft ? mat.rows() : mat.cols();
911  // FIXME we need an is_same for expression that is not sensitive to constness. For instance
912  // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
913  //if(is_same<MatrixTypeCleaned,Dest>::value && extract_data(dst) == extract_data(mat))
914  if(is_same_dense(dst, mat))
915  {
916  // apply the permutation inplace
917  Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(perm.size());
918  mask.fill(false);
919  Index r = 0;
920  while(r < perm.size())
921  {
922  // search for the next seed
923  while(r<perm.size() && mask[r]) r++;
924  if(r>=perm.size())
925  break;
926  // we got one, let's follow it until we are back to the seed
927  Index k0 = r++;
928  Index kPrev = k0;
929  mask.coeffRef(k0) = true;
930  for(Index k=perm.indices().coeff(k0); k!=k0; k=perm.indices().coeff(k))
931  {
932  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
933  .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
934  (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
935 
936  mask.coeffRef(k) = true;
937  kPrev = k;
938  }
939  }
940  }
941  else
942  {
943  for(Index i = 0; i < n; ++i)
944  {
945  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
946  (dst, ((Side==OnTheLeft) ^ Transposed) ? perm.indices().coeff(i) : i)
947 
948  =
949 
950  Block<const MatrixTypeCleaned,Side==OnTheLeft ? 1 : MatrixTypeCleaned::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixTypeCleaned::ColsAtCompileTime>
951  (mat, ((Side==OnTheRight) ^ Transposed) ? perm.indices().coeff(i) : i);
952  }
953  }
954  }
955 };
956 
957 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
958 struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>
959 {
960  template<typename Dest>
961  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
962  {
963  permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
964  }
965 };
966 
967 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
968 struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>
969 {
970  template<typename Dest>
971  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
972  {
973  permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
974  }
975 };
976 
977 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
978 struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>
979 {
980  template<typename Dest>
981  static void evalTo(Dest& dst, const Inverse<Lhs>& lhs, const Rhs& rhs)
982  {
983  permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
984  }
985 };
986 
987 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
988 struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag>
989 {
990  template<typename Dest>
991  static void evalTo(Dest& dst, const Lhs& lhs, const Inverse<Rhs>& rhs)
992  {
993  permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
994  }
995 };
996 
997 
998 /***************************************************************************
999 * Products with transpositions matrices
1000 ***************************************************************************/
1001 
1002 // FIXME could we unify Transpositions and Permutation into a single "shape"??
1003 
1008 template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
1009 struct transposition_matrix_product
1010 {
1011  typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
1012  typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
1013 
1014  template<typename Dest, typename TranspositionType>
1015  static inline void run(Dest& dst, const TranspositionType& tr, const ExpressionType& xpr)
1016  {
1017  MatrixType mat(xpr);
1018  typedef typename TranspositionType::StorageIndex StorageIndex;
1019  const Index size = tr.size();
1020  StorageIndex j = 0;
1021 
1022  if(!is_same_dense(dst,mat))
1023  dst = mat;
1024 
1025  for(Index k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
1026  if(Index(j=tr.coeff(k))!=k)
1027  {
1028  if(Side==OnTheLeft) dst.row(k).swap(dst.row(j));
1029  else if(Side==OnTheRight) dst.col(k).swap(dst.col(j));
1030  }
1031  }
1032 };
1033 
1034 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1035 struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1036 {
1037  template<typename Dest>
1038  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1039  {
1040  transposition_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
1041  }
1042 };
1043 
1044 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1045 struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag>
1046 {
1047  template<typename Dest>
1048  static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1049  {
1050  transposition_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
1051  }
1052 };
1053 
1054 
1055 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1056 struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1057 {
1058  template<typename Dest>
1059  static void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
1060  {
1061  transposition_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
1062  }
1063 };
1064 
1065 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1066 struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag>
1067 {
1068  template<typename Dest>
1069  static void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
1070  {
1071  transposition_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
1072  }
1073 };
1074 
1075 } // end namespace internal
1076 
1077 } // end namespace Eigen
1078 
1079 #endif // EIGEN_PRODUCT_EVALUATORS_H
Definition: Constants.h:320
const int HugeCost
Definition: Constants.h:39
Definition: Constants.h:335
Definition: Constants.h:230
Namespace containing all symbols from the Eigen library.
Definition: Core:271
const unsigned int RowMajorBit
Definition: Constants.h:61
const unsigned int PacketAccessBit
Definition: Constants.h:89
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
Definition: Constants.h:333
Definition: Eigen_Colamd.h:50
Definition: Constants.h:322
const int Dynamic
Definition: Constants.h:21
const unsigned int EvalBeforeNestingBit
Definition: Constants.h:65
const unsigned int ActualPacketAccessBit
Definition: Constants.h:100
const unsigned int LinearAccessBit
Definition: Constants.h:125