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