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