Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
16namespace Eigen {
17
18namespace internal {
19
28template<typename Lhs, typename Rhs, int Options>
29struct 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 EIGEN_STRONG_INLINE 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
40template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
41struct 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};
47template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
48struct 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 EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr)
59 : Base(xpr.lhs().functor().m_other * xpr.rhs().lhs() * xpr.rhs().rhs())
60 {}
61};
62
63
64template<typename Lhs, typename Rhs, int DiagIndex>
65struct 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 EIGEN_STRONG_INLINE 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.
82template< 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>
86struct generic_product_impl;
87
88template<typename Lhs, typename Rhs>
89struct 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
95template<typename Lhs, typename Rhs, int Options, int ProductTag, typename LhsShape, typename RhsShape>
96struct 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
127protected:
128 PlainObject m_result;
129};
130
131// The following three shortcuts are enabled only if the scalar types match exactly.
132// TODO: we could enable them for different scalar types when the product is not vectorized.
133
134// Dense = Product
135template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
136struct 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_DEVICE_FUNC EIGEN_STRONG_INLINE
141 void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
142 {
143 Index dstRows = src.rows();
144 Index dstCols = src.cols();
145 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
146 dst.resize(dstRows, dstCols);
147 // FIXME shall we handle nested_eval here?
148 generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
149 }
150};
151
152// Dense += Product
153template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
154struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::add_assign_op<Scalar,Scalar>, Dense2Dense,
155 typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
156{
157 typedef Product<Lhs,Rhs,Options> SrcXprType;
158 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
159 void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,Scalar> &)
160 {
161 eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
162 // FIXME shall we handle nested_eval here?
163 generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs());
164 }
165};
166
167// Dense -= Product
168template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
169struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::sub_assign_op<Scalar,Scalar>, Dense2Dense,
170 typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
171{
172 typedef Product<Lhs,Rhs,Options> SrcXprType;
173 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
174 void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,Scalar> &)
175 {
176 eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
177 // FIXME shall we handle nested_eval here?
178 generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs());
179 }
180};
181
182
183// Dense ?= scalar * Product
184// TODO we should apply that rule if that's really helpful
185// for instance, this is not good for inner products
186template< typename DstXprType, typename Lhs, typename Rhs, typename AssignFunc, typename Scalar, typename ScalarBis, typename Plain>
187struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>, const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
188 const Product<Lhs,Rhs,DefaultProduct> >, AssignFunc, Dense2Dense>
189{
190 typedef CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>,
191 const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
192 const Product<Lhs,Rhs,DefaultProduct> > SrcXprType;
193 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
194 void run(DstXprType &dst, const SrcXprType &src, const AssignFunc& func)
195 {
196 call_assignment_no_alias(dst, (src.lhs().functor().m_other * src.rhs().lhs())*src.rhs().rhs(), func);
197 }
198};
199
200//----------------------------------------
201// Catch "Dense ?= xpr + Product<>" expression to save one temporary
202// FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct
203
204template<typename OtherXpr, typename Lhs, typename Rhs>
205struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
206 const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
207 static const bool value = true;
208};
209
210template<typename OtherXpr, typename Lhs, typename Rhs>
211struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_difference_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
212 const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
213 static const bool value = true;
214};
215
216template<typename DstXprType, typename OtherXpr, typename ProductType, typename Func1, typename Func2>
217struct assignment_from_xpr_op_product
218{
219 template<typename SrcXprType, typename InitialFunc>
220 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
221 void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
222 {
223 call_assignment_no_alias(dst, src.lhs(), Func1());
224 call_assignment_no_alias(dst, src.rhs(), Func2());
225 }
226};
227
228#define EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(ASSIGN_OP,BINOP,ASSIGN_OP2) \
229 template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> \
230 struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<OtherScalar,ProdScalar>, const OtherXpr, \
231 const Product<Lhs,Rhs,DefaultProduct> >, internal::ASSIGN_OP<DstScalar,SrcScalar>, Dense2Dense> \
232 : assignment_from_xpr_op_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::ASSIGN_OP<DstScalar,OtherScalar>, internal::ASSIGN_OP2<DstScalar,ProdScalar> > \
233 {}
234
235EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op, scalar_sum_op,add_assign_op);
236EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_sum_op,add_assign_op);
237EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_sum_op,sub_assign_op);
238
239EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op, scalar_difference_op,sub_assign_op);
240EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_difference_op,sub_assign_op);
241EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_difference_op,add_assign_op);
242
243//----------------------------------------
244
245template<typename Lhs, typename Rhs>
246struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,InnerProduct>
247{
248 template<typename Dst>
249 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
250 {
251 dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
252 }
253
254 template<typename Dst>
255 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
256 {
257 dst.coeffRef(0,0) += (lhs.transpose().cwiseProduct(rhs)).sum();
258 }
259
260 template<typename Dst>
261 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
262 { dst.coeffRef(0,0) -= (lhs.transpose().cwiseProduct(rhs)).sum(); }
263};
264
265
266/***********************************************************************
267* Implementation of outer dense * dense vector product
268***********************************************************************/
269
270// Column major result
271template<typename Dst, typename Lhs, typename Rhs, typename Func>
272void EIGEN_DEVICE_FUNC outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const false_type&)
273{
274 evaluator<Rhs> rhsEval(rhs);
275 ei_declare_local_nested_eval(Lhs,lhs,Rhs::SizeAtCompileTime,actual_lhs);
276 // FIXME if cols is large enough, then it might be useful to make sure that lhs is sequentially stored
277 // FIXME not very good if rhs is real and lhs complex while alpha is real too
278 const Index cols = dst.cols();
279 for (Index j=0; j<cols; ++j)
280 func(dst.col(j), rhsEval.coeff(Index(0),j) * actual_lhs);
281}
282
283// Row major result
284template<typename Dst, typename Lhs, typename Rhs, typename Func>
285void EIGEN_DEVICE_FUNC outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const true_type&)
286{
287 evaluator<Lhs> lhsEval(lhs);
288 ei_declare_local_nested_eval(Rhs,rhs,Lhs::SizeAtCompileTime,actual_rhs);
289 // FIXME if rows is large enough, then it might be useful to make sure that rhs is sequentially stored
290 // FIXME not very good if lhs is real and rhs complex while alpha is real too
291 const Index rows = dst.rows();
292 for (Index i=0; i<rows; ++i)
293 func(dst.row(i), lhsEval.coeff(i,Index(0)) * actual_rhs);
294}
295
296template<typename Lhs, typename Rhs>
297struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,OuterProduct>
298{
299 template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
300 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
301
302 // TODO it would be nice to be able to exploit our *_assign_op functors for that purpose
303 struct set { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } };
304 struct add { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
305 struct sub { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
306 struct adds {
307 Scalar m_scale;
308 explicit adds(const Scalar& s) : m_scale(s) {}
309 template<typename Dst, typename Src> void EIGEN_DEVICE_FUNC operator()(const Dst& dst, const Src& src) const {
310 dst.const_cast_derived() += m_scale * src;
311 }
312 };
313
314 template<typename Dst>
315 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
316 {
317 internal::outer_product_selector_run(dst, lhs, rhs, set(), is_row_major<Dst>());
318 }
319
320 template<typename Dst>
321 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
322 {
323 internal::outer_product_selector_run(dst, lhs, rhs, add(), is_row_major<Dst>());
324 }
325
326 template<typename Dst>
327 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
328 {
329 internal::outer_product_selector_run(dst, lhs, rhs, sub(), is_row_major<Dst>());
330 }
331
332 template<typename Dst>
333 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
334 {
335 internal::outer_product_selector_run(dst, lhs, rhs, adds(alpha), is_row_major<Dst>());
336 }
337
338};
339
340
341// This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo
342template<typename Lhs, typename Rhs, typename Derived>
343struct generic_product_impl_base
344{
345 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
346
347 template<typename Dst>
348 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
349 { dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
350
351 template<typename Dst>
352 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
353 { scaleAndAddTo(dst,lhs, rhs, Scalar(1)); }
354
355 template<typename Dst>
356 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
357 { scaleAndAddTo(dst, lhs, rhs, Scalar(-1)); }
358
359 template<typename Dst>
360 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
361 { Derived::scaleAndAddTo(dst,lhs,rhs,alpha); }
362
363};
364
365template<typename Lhs, typename Rhs>
366struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct>
367 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct> >
368{
369 typedef typename nested_eval<Lhs,1>::type LhsNested;
370 typedef typename nested_eval<Rhs,1>::type RhsNested;
371 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
372 enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
373 typedef typename internal::remove_all<typename internal::conditional<int(Side)==OnTheRight,LhsNested,RhsNested>::type>::type MatrixType;
374
375 template<typename Dest>
376 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
377 {
378 // Fallback to inner product if both the lhs and rhs is a runtime vector.
379 if (lhs.rows() == 1 && rhs.cols() == 1) {
380 dst.coeffRef(0,0) += alpha * lhs.row(0).conjugate().dot(rhs.col(0));
381 return;
382 }
383 LhsNested actual_lhs(lhs);
384 RhsNested actual_rhs(rhs);
385 internal::gemv_dense_selector<Side,
386 (int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
387 bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)
388 >::run(actual_lhs, actual_rhs, dst, alpha);
389 }
390};
391
392template<typename Lhs, typename Rhs>
393struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode>
394{
395 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
396
397 template<typename Dst>
398 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
399 {
400 // Same as: dst.noalias() = lhs.lazyProduct(rhs);
401 // but easier on the compiler side
402 call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op<typename Dst::Scalar,Scalar>());
403 }
404
405 template<typename Dst>
406 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
407 {
408 // dst.noalias() += lhs.lazyProduct(rhs);
409 call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op<typename Dst::Scalar,Scalar>());
410 }
411
412 template<typename Dst>
413 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
414 {
415 // dst.noalias() -= lhs.lazyProduct(rhs);
416 call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op<typename Dst::Scalar,Scalar>());
417 }
418
419 // This is a special evaluation path called from generic_product_impl<...,GemmProduct> in file GeneralMatrixMatrix.h
420 // This variant tries to extract scalar multiples from both the LHS and RHS and factor them out. For instance:
421 // dst {,+,-}= (s1*A)*(B*s2)
422 // will be rewritten as:
423 // dst {,+,-}= (s1*s2) * (A.lazyProduct(B))
424 // There are at least four benefits of doing so:
425 // 1 - huge performance gain for heap-allocated matrix types as it save costly allocations.
426 // 2 - it is faster than simply by-passing the heap allocation through stack allocation.
427 // 3 - it makes this fallback consistent with the heavy GEMM routine.
428 // 4 - it fully by-passes huge stack allocation attempts when multiplying huge fixed-size matrices.
429 // (see https://stackoverflow.com/questions/54738495)
430 // For small fixed sizes matrices, howver, the gains are less obvious, it is sometimes x2 faster, but sometimes x3 slower,
431 // and the behavior depends also a lot on the compiler... This is why this re-writting strategy is currently
432 // enabled only when falling back from the main GEMM.
433 template<typename Dst, typename Func>
434 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
435 void eval_dynamic(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Func &func)
436 {
437 enum {
438 HasScalarFactor = blas_traits<Lhs>::HasScalarFactor || blas_traits<Rhs>::HasScalarFactor,
439 ConjLhs = blas_traits<Lhs>::NeedToConjugate,
440 ConjRhs = blas_traits<Rhs>::NeedToConjugate
441 };
442 // FIXME: in c++11 this should be auto, and extractScalarFactor should also return auto
443 // this is important for real*complex_mat
444 Scalar actualAlpha = combine_scalar_factors<Scalar>(lhs, rhs);
445
446 eval_dynamic_impl(dst,
447 blas_traits<Lhs>::extract(lhs).template conjugateIf<ConjLhs>(),
448 blas_traits<Rhs>::extract(rhs).template conjugateIf<ConjRhs>(),
449 func,
450 actualAlpha,
451 typename conditional<HasScalarFactor,true_type,false_type>::type());
452 }
453
454protected:
455
456 template<typename Dst, typename LhsT, typename RhsT, typename Func, typename Scalar>
457 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
458 void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar& s /* == 1 */, false_type)
459 {
460 EIGEN_UNUSED_VARIABLE(s);
461 eigen_internal_assert(s==Scalar(1));
462 call_restricted_packet_assignment_no_alias(dst, lhs.lazyProduct(rhs), func);
463 }
464
465 template<typename Dst, typename LhsT, typename RhsT, typename Func, typename Scalar>
466 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
467 void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar& s, true_type)
468 {
469 call_restricted_packet_assignment_no_alias(dst, s * lhs.lazyProduct(rhs), func);
470 }
471};
472
473// This specialization enforces the use of a coefficient-based evaluation strategy
474template<typename Lhs, typename Rhs>
475struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,LazyCoeffBasedProductMode>
476 : generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode> {};
477
478// Case 2: Evaluate coeff by coeff
479//
480// This is mostly taken from CoeffBasedProduct.h
481// The main difference is that we add an extra argument to the etor_product_*_impl::run() function
482// for the inner dimension of the product, because evaluator object do not know their size.
483
484template<int Traversal, int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar>
485struct etor_product_coeff_impl;
486
487template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
488struct etor_product_packet_impl;
489
490template<typename Lhs, typename Rhs, int ProductTag>
491struct product_evaluator<Product<Lhs, Rhs, LazyProduct>, ProductTag, DenseShape, DenseShape>
492 : evaluator_base<Product<Lhs, Rhs, LazyProduct> >
493{
494 typedef Product<Lhs, Rhs, LazyProduct> XprType;
495 typedef typename XprType::Scalar Scalar;
496 typedef typename XprType::CoeffReturnType CoeffReturnType;
497
498 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
499 explicit product_evaluator(const XprType& xpr)
500 : m_lhs(xpr.lhs()),
501 m_rhs(xpr.rhs()),
502 m_lhsImpl(m_lhs), // FIXME the creation of the evaluator objects should result in a no-op, but check that!
503 m_rhsImpl(m_rhs), // Moreover, they are only useful for the packet path, so we could completely disable them when not needed,
504 // or perhaps declare them on the fly on the packet method... We have experiment to check what's best.
505 m_innerDim(xpr.lhs().cols())
506 {
507 EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
508 EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::AddCost);
509 EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
510#if 0
511 std::cerr << "LhsOuterStrideBytes= " << LhsOuterStrideBytes << "\n";
512 std::cerr << "RhsOuterStrideBytes= " << RhsOuterStrideBytes << "\n";
513 std::cerr << "LhsAlignment= " << LhsAlignment << "\n";
514 std::cerr << "RhsAlignment= " << RhsAlignment << "\n";
515 std::cerr << "CanVectorizeLhs= " << CanVectorizeLhs << "\n";
516 std::cerr << "CanVectorizeRhs= " << CanVectorizeRhs << "\n";
517 std::cerr << "CanVectorizeInner= " << CanVectorizeInner << "\n";
518 std::cerr << "EvalToRowMajor= " << EvalToRowMajor << "\n";
519 std::cerr << "Alignment= " << Alignment << "\n";
520 std::cerr << "Flags= " << Flags << "\n";
521#endif
522 }
523
524 // Everything below here is taken from CoeffBasedProduct.h
525
526 typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
527 typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
528
529 typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
530 typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
531
532 typedef evaluator<LhsNestedCleaned> LhsEtorType;
533 typedef evaluator<RhsNestedCleaned> RhsEtorType;
534
535 enum {
536 RowsAtCompileTime = LhsNestedCleaned::RowsAtCompileTime,
537 ColsAtCompileTime = RhsNestedCleaned::ColsAtCompileTime,
538 InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsNestedCleaned::ColsAtCompileTime, RhsNestedCleaned::RowsAtCompileTime),
539 MaxRowsAtCompileTime = LhsNestedCleaned::MaxRowsAtCompileTime,
540 MaxColsAtCompileTime = RhsNestedCleaned::MaxColsAtCompileTime
541 };
542
543 typedef typename find_best_packet<Scalar,RowsAtCompileTime>::type LhsVecPacketType;
544 typedef typename find_best_packet<Scalar,ColsAtCompileTime>::type RhsVecPacketType;
545
546 enum {
547
548 LhsCoeffReadCost = LhsEtorType::CoeffReadCost,
549 RhsCoeffReadCost = RhsEtorType::CoeffReadCost,
550 CoeffReadCost = InnerSize==0 ? NumTraits<Scalar>::ReadCost
551 : InnerSize == Dynamic ? HugeCost
552 : InnerSize * (NumTraits<Scalar>::MulCost + int(LhsCoeffReadCost) + int(RhsCoeffReadCost))
553 + (InnerSize - 1) * NumTraits<Scalar>::AddCost,
554
555 Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
556
557 LhsFlags = LhsEtorType::Flags,
558 RhsFlags = RhsEtorType::Flags,
559
560 LhsRowMajor = LhsFlags & RowMajorBit,
561 RhsRowMajor = RhsFlags & RowMajorBit,
562
563 LhsVecPacketSize = unpacket_traits<LhsVecPacketType>::size,
564 RhsVecPacketSize = unpacket_traits<RhsVecPacketType>::size,
565
566 // Here, we don't care about alignment larger than the usable packet size.
567 LhsAlignment = EIGEN_PLAIN_ENUM_MIN(LhsEtorType::Alignment,LhsVecPacketSize*int(sizeof(typename LhsNestedCleaned::Scalar))),
568 RhsAlignment = EIGEN_PLAIN_ENUM_MIN(RhsEtorType::Alignment,RhsVecPacketSize*int(sizeof(typename RhsNestedCleaned::Scalar))),
569
570 SameType = is_same<typename LhsNestedCleaned::Scalar,typename RhsNestedCleaned::Scalar>::value,
571
572 CanVectorizeRhs = bool(RhsRowMajor) && (RhsFlags & PacketAccessBit) && (ColsAtCompileTime!=1),
573 CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) && (RowsAtCompileTime!=1),
574
575 EvalToRowMajor = (MaxRowsAtCompileTime==1&&MaxColsAtCompileTime!=1) ? 1
576 : (MaxColsAtCompileTime==1&&MaxRowsAtCompileTime!=1) ? 0
577 : (bool(RhsRowMajor) && !CanVectorizeLhs),
578
579 Flags = ((int(LhsFlags) | int(RhsFlags)) & HereditaryBits & ~RowMajorBit)
580 | (EvalToRowMajor ? RowMajorBit : 0)
581 // TODO enable vectorization for mixed types
582 | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0)
583 | (XprType::IsVectorAtCompileTime ? LinearAccessBit : 0),
584
585 LhsOuterStrideBytes = int(LhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename LhsNestedCleaned::Scalar)),
586 RhsOuterStrideBytes = int(RhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename RhsNestedCleaned::Scalar)),
587
588 Alignment = bool(CanVectorizeLhs) ? (LhsOuterStrideBytes<=0 || (int(LhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,LhsAlignment))!=0 ? 0 : LhsAlignment)
589 : bool(CanVectorizeRhs) ? (RhsOuterStrideBytes<=0 || (int(RhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,RhsAlignment))!=0 ? 0 : RhsAlignment)
590 : 0,
591
592 /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside
593 * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner
594 * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
595 * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
596 */
597 CanVectorizeInner = SameType
598 && LhsRowMajor
599 && (!RhsRowMajor)
600 && (int(LhsFlags) & int(RhsFlags) & ActualPacketAccessBit)
601 && (int(InnerSize) % packet_traits<Scalar>::size == 0)
602 };
603
604 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index row, Index col) const
605 {
606 return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
607 }
608
609 /* Allow index-based non-packet access. It is impossible though to allow index-based packed access,
610 * which is why we don't set the LinearAccessBit.
611 * TODO: this seems possible when the result is a vector
612 */
613 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
614 const CoeffReturnType coeff(Index index) const
615 {
616 const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
617 const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
618 return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
619 }
620
621 template<int LoadMode, typename PacketType>
622 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
623 const PacketType packet(Index row, Index col) const
624 {
625 PacketType res;
626 typedef etor_product_packet_impl<bool(int(Flags)&RowMajorBit) ? RowMajor : ColMajor,
627 Unroll ? int(InnerSize) : Dynamic,
628 LhsEtorType, RhsEtorType, PacketType, LoadMode> PacketImpl;
629 PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res);
630 return res;
631 }
632
633 template<int LoadMode, typename PacketType>
634 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
635 const PacketType packet(Index index) const
636 {
637 const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
638 const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
639 return packet<LoadMode,PacketType>(row,col);
640 }
641
642protected:
643 typename internal::add_const_on_value_type<LhsNested>::type m_lhs;
644 typename internal::add_const_on_value_type<RhsNested>::type m_rhs;
645
646 LhsEtorType m_lhsImpl;
647 RhsEtorType m_rhsImpl;
648
649 // TODO: Get rid of m_innerDim if known at compile time
650 Index m_innerDim;
651};
652
653template<typename Lhs, typename Rhs>
654struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, LazyCoeffBasedProductMode, DenseShape, DenseShape>
655 : product_evaluator<Product<Lhs, Rhs, LazyProduct>, CoeffBasedProductMode, DenseShape, DenseShape>
656{
657 typedef Product<Lhs, Rhs, DefaultProduct> XprType;
658 typedef Product<Lhs, Rhs, LazyProduct> BaseProduct;
659 typedef product_evaluator<BaseProduct, CoeffBasedProductMode, DenseShape, DenseShape> Base;
660 enum {
661 Flags = Base::Flags | EvalBeforeNestingBit
662 };
663 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
664 explicit product_evaluator(const XprType& xpr)
665 : Base(BaseProduct(xpr.lhs(),xpr.rhs()))
666 {}
667};
668
669/****************************************
670*** Coeff based product, Packet path ***
671****************************************/
672
673template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
674struct etor_product_packet_impl<RowMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
675{
676 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
677 {
678 etor_product_packet_impl<RowMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
679 res = pmadd(pset1<Packet>(lhs.coeff(row, Index(UnrollingIndex-1))), rhs.template packet<LoadMode,Packet>(Index(UnrollingIndex-1), col), res);
680 }
681};
682
683template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
684struct etor_product_packet_impl<ColMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
685{
686 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
687 {
688 etor_product_packet_impl<ColMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
689 res = pmadd(lhs.template packet<LoadMode,Packet>(row, Index(UnrollingIndex-1)), pset1<Packet>(rhs.coeff(Index(UnrollingIndex-1), col)), res);
690 }
691};
692
693template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
694struct etor_product_packet_impl<RowMajor, 1, Lhs, Rhs, Packet, LoadMode>
695{
696 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
697 {
698 res = pmul(pset1<Packet>(lhs.coeff(row, Index(0))),rhs.template packet<LoadMode,Packet>(Index(0), col));
699 }
700};
701
702template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
703struct etor_product_packet_impl<ColMajor, 1, Lhs, Rhs, Packet, LoadMode>
704{
705 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
706 {
707 res = pmul(lhs.template packet<LoadMode,Packet>(row, Index(0)), pset1<Packet>(rhs.coeff(Index(0), col)));
708 }
709};
710
711template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
712struct etor_product_packet_impl<RowMajor, 0, Lhs, Rhs, Packet, LoadMode>
713{
714 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
715 {
716 res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
717 }
718};
719
720template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
721struct etor_product_packet_impl<ColMajor, 0, Lhs, Rhs, Packet, LoadMode>
722{
723 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
724 {
725 res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
726 }
727};
728
729template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
730struct etor_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
731{
732 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
733 {
734 res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
735 for(Index i = 0; i < innerDim; ++i)
736 res = pmadd(pset1<Packet>(lhs.coeff(row, i)), rhs.template packet<LoadMode,Packet>(i, col), res);
737 }
738};
739
740template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
741struct etor_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
742{
743 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
744 {
745 res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
746 for(Index i = 0; i < innerDim; ++i)
747 res = pmadd(lhs.template packet<LoadMode,Packet>(row, i), pset1<Packet>(rhs.coeff(i, col)), res);
748 }
749};
750
751
752/***************************************************************************
753* Triangular products
754***************************************************************************/
755template<int Mode, bool LhsIsTriangular,
756 typename Lhs, bool LhsIsVector,
757 typename Rhs, bool RhsIsVector>
758struct triangular_product_impl;
759
760template<typename Lhs, typename Rhs, int ProductTag>
761struct generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag>
762 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag> >
763{
764 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
765
766 template<typename Dest>
767 static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
768 {
769 triangular_product_impl<Lhs::Mode,true,typename Lhs::MatrixType,false,Rhs, Rhs::ColsAtCompileTime==1>
770 ::run(dst, lhs.nestedExpression(), rhs, alpha);
771 }
772};
773
774template<typename Lhs, typename Rhs, int ProductTag>
775struct generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag>
776: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag> >
777{
778 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
779
780 template<typename Dest>
781 static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
782 {
783 triangular_product_impl<Rhs::Mode,false,Lhs,Lhs::RowsAtCompileTime==1, typename Rhs::MatrixType, false>::run(dst, lhs, rhs.nestedExpression(), alpha);
784 }
785};
786
787
788/***************************************************************************
789* SelfAdjoint products
790***************************************************************************/
791template <typename Lhs, int LhsMode, bool LhsIsVector,
792 typename Rhs, int RhsMode, bool RhsIsVector>
793struct selfadjoint_product_impl;
794
795template<typename Lhs, typename Rhs, int ProductTag>
796struct generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag>
797 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag> >
798{
799 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
800
801 template<typename Dest>
802 static EIGEN_DEVICE_FUNC
803 void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
804 {
805 selfadjoint_product_impl<typename Lhs::MatrixType,Lhs::Mode,false,Rhs,0,Rhs::IsVectorAtCompileTime>::run(dst, lhs.nestedExpression(), rhs, alpha);
806 }
807};
808
809template<typename Lhs, typename Rhs, int ProductTag>
810struct generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag>
811: generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag> >
812{
813 typedef typename Product<Lhs,Rhs>::Scalar Scalar;
814
815 template<typename Dest>
816 static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
817 {
818 selfadjoint_product_impl<Lhs,0,Lhs::IsVectorAtCompileTime,typename Rhs::MatrixType,Rhs::Mode,false>::run(dst, lhs, rhs.nestedExpression(), alpha);
819 }
820};
821
822
823/***************************************************************************
824* Diagonal products
825***************************************************************************/
826
827template<typename MatrixType, typename DiagonalType, typename Derived, int ProductOrder>
828struct diagonal_product_evaluator_base
829 : evaluator_base<Derived>
830{
831 typedef typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar;
832public:
833 enum {
834 CoeffReadCost = int(NumTraits<Scalar>::MulCost) + int(evaluator<MatrixType>::CoeffReadCost) + int(evaluator<DiagonalType>::CoeffReadCost),
835
836 MatrixFlags = evaluator<MatrixType>::Flags,
837 DiagFlags = evaluator<DiagonalType>::Flags,
838
839 _StorageOrder = (Derived::MaxRowsAtCompileTime==1 && Derived::MaxColsAtCompileTime!=1) ? RowMajor
840 : (Derived::MaxColsAtCompileTime==1 && Derived::MaxRowsAtCompileTime!=1) ? ColMajor
841 : MatrixFlags & RowMajorBit ? RowMajor : ColMajor,
842 _SameStorageOrder = _StorageOrder == (MatrixFlags & RowMajorBit ? RowMajor : ColMajor),
843
844 _ScalarAccessOnDiag = !((int(_StorageOrder) == ColMajor && int(ProductOrder) == OnTheLeft)
845 ||(int(_StorageOrder) == RowMajor && int(ProductOrder) == OnTheRight)),
846 _SameTypes = is_same<typename MatrixType::Scalar, typename DiagonalType::Scalar>::value,
847 // FIXME currently we need same types, but in the future the next rule should be the one
848 //_Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && ((!_PacketOnDiag) || (_SameTypes && bool(int(DiagFlags)&PacketAccessBit))),
849 _Vectorizable = bool(int(MatrixFlags)&PacketAccessBit)
850 && _SameTypes
851 && (_SameStorageOrder || (MatrixFlags&LinearAccessBit)==LinearAccessBit)
852 && (_ScalarAccessOnDiag || (bool(int(DiagFlags)&PacketAccessBit))),
853 _LinearAccessMask = (MatrixType::RowsAtCompileTime==1 || MatrixType::ColsAtCompileTime==1) ? LinearAccessBit : 0,
854 Flags = ((HereditaryBits|_LinearAccessMask) & (unsigned int)(MatrixFlags)) | (_Vectorizable ? PacketAccessBit : 0),
855 Alignment = evaluator<MatrixType>::Alignment,
856
857 AsScalarProduct = (DiagonalType::SizeAtCompileTime==1)
858 || (DiagonalType::SizeAtCompileTime==Dynamic && MatrixType::RowsAtCompileTime==1 && ProductOrder==OnTheLeft)
859 || (DiagonalType::SizeAtCompileTime==Dynamic && MatrixType::ColsAtCompileTime==1 && ProductOrder==OnTheRight)
860 };
861
862 EIGEN_DEVICE_FUNC diagonal_product_evaluator_base(const MatrixType &mat, const DiagonalType &diag)
863 : m_diagImpl(diag), m_matImpl(mat)
864 {
865 EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
866 EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
867 }
868
869 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index idx) const
870 {
871 if(AsScalarProduct)
872 return m_diagImpl.coeff(0) * m_matImpl.coeff(idx);
873 else
874 return m_diagImpl.coeff(idx) * m_matImpl.coeff(idx);
875 }
876
877protected:
878 template<int LoadMode,typename PacketType>
879 EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::true_type) const
880 {
881 return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
882 internal::pset1<PacketType>(m_diagImpl.coeff(id)));
883 }
884
885 template<int LoadMode,typename PacketType>
886 EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::false_type) const
887 {
888 enum {
889 InnerSize = (MatrixType::Flags & RowMajorBit) ? MatrixType::ColsAtCompileTime : MatrixType::RowsAtCompileTime,
890 DiagonalPacketLoadMode = EIGEN_PLAIN_ENUM_MIN(LoadMode,((InnerSize%16) == 0) ? int(Aligned16) : int(evaluator<DiagonalType>::Alignment)) // FIXME hardcoded 16!!
891 };
892 return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
893 m_diagImpl.template packet<DiagonalPacketLoadMode,PacketType>(id));
894 }
895
896 evaluator<DiagonalType> m_diagImpl;
897 evaluator<MatrixType> m_matImpl;
898};
899
900// diagonal * dense
901template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
902struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, DenseShape>
903 : diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft>
904{
905 typedef diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft> Base;
906 using Base::m_diagImpl;
907 using Base::m_matImpl;
908 using Base::coeff;
909 typedef typename Base::Scalar Scalar;
910
911 typedef Product<Lhs, Rhs, ProductKind> XprType;
912 typedef typename XprType::PlainObject PlainObject;
913 typedef typename Lhs::DiagonalVectorType DiagonalType;
914
915
916 enum { StorageOrder = Base::_StorageOrder };
917
918 EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
919 : Base(xpr.rhs(), xpr.lhs().diagonal())
920 {
921 }
922
923 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
924 {
925 return m_diagImpl.coeff(row) * m_matImpl.coeff(row, col);
926 }
927
928#ifndef EIGEN_GPUCC
929 template<int LoadMode,typename PacketType>
930 EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
931 {
932 // FIXME: NVCC used to complain about the template keyword, but we have to check whether this is still the case.
933 // See also similar calls below.
934 return this->template packet_impl<LoadMode,PacketType>(row,col, row,
935 typename internal::conditional<int(StorageOrder)==RowMajor, internal::true_type, internal::false_type>::type());
936 }
937
938 template<int LoadMode,typename PacketType>
939 EIGEN_STRONG_INLINE PacketType packet(Index idx) const
940 {
941 return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
942 }
943#endif
944};
945
946// dense * diagonal
947template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
948struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DenseShape, DiagonalShape>
949 : diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight>
950{
951 typedef diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight> Base;
952 using Base::m_diagImpl;
953 using Base::m_matImpl;
954 using Base::coeff;
955 typedef typename Base::Scalar Scalar;
956
957 typedef Product<Lhs, Rhs, ProductKind> XprType;
958 typedef typename XprType::PlainObject PlainObject;
959
960 enum { StorageOrder = Base::_StorageOrder };
961
962 EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
963 : Base(xpr.lhs(), xpr.rhs().diagonal())
964 {
965 }
966
967 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
968 {
969 return m_matImpl.coeff(row, col) * m_diagImpl.coeff(col);
970 }
971
972#ifndef EIGEN_GPUCC
973 template<int LoadMode,typename PacketType>
974 EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
975 {
976 return this->template packet_impl<LoadMode,PacketType>(row,col, col,
977 typename internal::conditional<int(StorageOrder)==ColMajor, internal::true_type, internal::false_type>::type());
978 }
979
980 template<int LoadMode,typename PacketType>
981 EIGEN_STRONG_INLINE PacketType packet(Index idx) const
982 {
983 return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
984 }
985#endif
986};
987
988/***************************************************************************
989* Products with permutation matrices
990***************************************************************************/
991
997template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
998struct permutation_matrix_product;
999
1000template<typename ExpressionType, int Side, bool Transposed>
1001struct permutation_matrix_product<ExpressionType, Side, Transposed, DenseShape>
1002{
1003 typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
1004 typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
1005
1006 template<typename Dest, typename PermutationType>
1007 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
1008 {
1009 MatrixType mat(xpr);
1010 const Index n = Side==OnTheLeft ? mat.rows() : mat.cols();
1011 // FIXME we need an is_same for expression that is not sensitive to constness. For instance
1012 // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
1013 //if(is_same<MatrixTypeCleaned,Dest>::value && extract_data(dst) == extract_data(mat))
1014 if(is_same_dense(dst, mat))
1015 {
1016 // apply the permutation inplace
1017 Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(perm.size());
1018 mask.fill(false);
1019 Index r = 0;
1020 while(r < perm.size())
1021 {
1022 // search for the next seed
1023 while(r<perm.size() && mask[r]) r++;
1024 if(r>=perm.size())
1025 break;
1026 // we got one, let's follow it until we are back to the seed
1027 Index k0 = r++;
1028 Index kPrev = k0;
1029 mask.coeffRef(k0) = true;
1030 for(Index k=perm.indices().coeff(k0); k!=k0; k=perm.indices().coeff(k))
1031 {
1032 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
1033 .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
1034 (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
1035
1036 mask.coeffRef(k) = true;
1037 kPrev = k;
1038 }
1039 }
1040 }
1041 else
1042 {
1043 for(Index i = 0; i < n; ++i)
1044 {
1045 Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
1046 (dst, ((Side==OnTheLeft) ^ Transposed) ? perm.indices().coeff(i) : i)
1047
1048 =
1049
1050 Block<const MatrixTypeCleaned,Side==OnTheLeft ? 1 : MatrixTypeCleaned::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixTypeCleaned::ColsAtCompileTime>
1051 (mat, ((Side==OnTheRight) ^ Transposed) ? perm.indices().coeff(i) : i);
1052 }
1053 }
1054 }
1055};
1056
1057template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1058struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>
1059{
1060 template<typename Dest>
1061 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1062 {
1063 permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
1064 }
1065};
1066
1067template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1068struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>
1069{
1070 template<typename Dest>
1071 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1072 {
1073 permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
1074 }
1075};
1076
1077template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1078struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>
1079{
1080 template<typename Dest>
1081 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Inverse<Lhs>& lhs, const Rhs& rhs)
1082 {
1083 permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
1084 }
1085};
1086
1087template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1088struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag>
1089{
1090 template<typename Dest>
1091 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Inverse<Rhs>& rhs)
1092 {
1093 permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
1094 }
1095};
1096
1097
1098/***************************************************************************
1099* Products with transpositions matrices
1100***************************************************************************/
1101
1102// FIXME could we unify Transpositions and Permutation into a single "shape"??
1103
1108template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
1109struct transposition_matrix_product
1110{
1111 typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
1112 typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
1113
1114 template<typename Dest, typename TranspositionType>
1115 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const TranspositionType& tr, const ExpressionType& xpr)
1116 {
1117 MatrixType mat(xpr);
1118 typedef typename TranspositionType::StorageIndex StorageIndex;
1119 const Index size = tr.size();
1120 StorageIndex j = 0;
1121
1122 if(!is_same_dense(dst,mat))
1123 dst = mat;
1124
1125 for(Index k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
1126 if(Index(j=tr.coeff(k))!=k)
1127 {
1128 if(Side==OnTheLeft) dst.row(k).swap(dst.row(j));
1129 else if(Side==OnTheRight) dst.col(k).swap(dst.col(j));
1130 }
1131 }
1132};
1133
1134template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1135struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1136{
1137 template<typename Dest>
1138 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1139 {
1140 transposition_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
1141 }
1142};
1143
1144template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1145struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag>
1146{
1147 template<typename Dest>
1148 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1149 {
1150 transposition_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
1151 }
1152};
1153
1154
1155template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1156struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1157{
1158 template<typename Dest>
1159 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
1160 {
1161 transposition_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
1162 }
1163};
1164
1165template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1166struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag>
1167{
1168 template<typename Dest>
1169 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
1170 {
1171 transposition_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
1172 }
1173};
1174
1175} // end namespace internal
1176
1177} // end namespace Eigen
1178
1179#endif // EIGEN_PRODUCT_EVALUATORS_H
@ ColMajor
Definition: Constants.h:319
@ RowMajor
Definition: Constants.h:321
@ OnTheLeft
Definition: Constants.h:332
@ OnTheRight
Definition: Constants.h:334
const unsigned int ActualPacketAccessBit
Definition: Constants.h:105
const unsigned int PacketAccessBit
Definition: Constants.h:94
const unsigned int LinearAccessBit
Definition: Constants.h:130
const unsigned int EvalBeforeNestingBit
Definition: Constants.h:70
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const int HugeCost
Definition: Constants.h:44
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22