1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_XPRHELPER_H
12 #define EIGEN_XPRHELPER_H
13 
14 // just a workaround because GCC seems to not really like empty structs
15 // FIXME: gcc 4.3 generates bad code when strict-aliasing is enabled
16 // so currently we simply disable this optimization for gcc 4.3
17 #if EIGEN_COMP_GNUC && !EIGEN_GNUC_AT(4,3)
18   #define EIGEN_EMPTY_STRUCT_CTOR(X) \
19     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE X() {} \
20     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE X(const X& ) {}
21 #else
22   #define EIGEN_EMPTY_STRUCT_CTOR(X)
23 #endif
24 
25 namespace Eigen {
26 
27 namespace internal {
28 
29 template<typename IndexDest, typename IndexSrc>
30 EIGEN_DEVICE_FUNC
convert_index(const IndexSrc & idx)31 inline IndexDest convert_index(const IndexSrc& idx) {
32   // for sizeof(IndexDest)>=sizeof(IndexSrc) compilers should be able to optimize this away:
33   eigen_internal_assert(idx <= NumTraits<IndexDest>::highest() && "Index value to big for target type");
34   return IndexDest(idx);
35 }
36 
37 
38 // promote_scalar_arg is an helper used in operation between an expression and a scalar, like:
39 //    expression * scalar
40 // Its role is to determine how the type T of the scalar operand should be promoted given the scalar type ExprScalar of the given expression.
41 // The IsSupported template parameter must be provided by the caller as: internal::has_ReturnType<ScalarBinaryOpTraits<ExprScalar,T,op> >::value using the proper order for ExprScalar and T.
42 // Then the logic is as follows:
43 //  - if the operation is natively supported as defined by IsSupported, then the scalar type is not promoted, and T is returned.
44 //  - otherwise, NumTraits<ExprScalar>::Literal is returned if T is implicitly convertible to NumTraits<ExprScalar>::Literal AND that this does not imply a float to integer conversion.
45 //  - otherwise, ExprScalar is returned if T is implicitly convertible to ExprScalar AND that this does not imply a float to integer conversion.
46 //  - In all other cases, the promoted type is not defined, and the respective operation is thus invalid and not available (SFINAE).
47 template<typename ExprScalar,typename T, bool IsSupported>
48 struct promote_scalar_arg;
49 
50 template<typename S,typename T>
51 struct promote_scalar_arg<S,T,true>
52 {
53   typedef T type;
54 };
55 
56 // Recursively check safe conversion to PromotedType, and then ExprScalar if they are different.
57 template<typename ExprScalar,typename T,typename PromotedType,
58   bool ConvertibleToLiteral = internal::is_convertible<T,PromotedType>::value,
59   bool IsSafe = NumTraits<T>::IsInteger || !NumTraits<PromotedType>::IsInteger>
60 struct promote_scalar_arg_unsupported;
61 
62 // Start recursion with NumTraits<ExprScalar>::Literal
63 template<typename S,typename T>
64 struct promote_scalar_arg<S,T,false> : promote_scalar_arg_unsupported<S,T,typename NumTraits<S>::Literal> {};
65 
66 // We found a match!
67 template<typename S,typename T, typename PromotedType>
68 struct promote_scalar_arg_unsupported<S,T,PromotedType,true,true>
69 {
70   typedef PromotedType type;
71 };
72 
73 // No match, but no real-to-integer issues, and ExprScalar and current PromotedType are different,
74 // so let's try to promote to ExprScalar
75 template<typename ExprScalar,typename T, typename PromotedType>
76 struct promote_scalar_arg_unsupported<ExprScalar,T,PromotedType,false,true>
77    : promote_scalar_arg_unsupported<ExprScalar,T,ExprScalar>
78 {};
79 
80 // Unsafe real-to-integer, let's stop.
81 template<typename S,typename T, typename PromotedType, bool ConvertibleToLiteral>
82 struct promote_scalar_arg_unsupported<S,T,PromotedType,ConvertibleToLiteral,false> {};
83 
84 // T is not even convertible to ExprScalar, let's stop.
85 template<typename S,typename T>
86 struct promote_scalar_arg_unsupported<S,T,S,false,true> {};
87 
88 //classes inheriting no_assignment_operator don't generate a default operator=.
89 class no_assignment_operator
90 {
91   private:
92     no_assignment_operator& operator=(const no_assignment_operator&);
93 };
94 
95 /** \internal return the index type with the largest number of bits */
96 template<typename I1, typename I2>
97 struct promote_index_type
98 {
99   typedef typename conditional<(sizeof(I1)<sizeof(I2)), I2, I1>::type type;
100 };
101 
102 /** \internal If the template parameter Value is Dynamic, this class is just a wrapper around a T variable that
103   * can be accessed using value() and setValue().
104   * Otherwise, this class is an empty structure and value() just returns the template parameter Value.
105   */
106 template<typename T, int Value> class variable_if_dynamic
107 {
108   public:
109     EIGEN_EMPTY_STRUCT_CTOR(variable_if_dynamic)
110     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit variable_if_dynamic(T v) { EIGEN_ONLY_USED_FOR_DEBUG(v); eigen_assert(v == T(Value)); }
111     EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T value() { return T(Value); }
112     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void setValue(T) {}
113 };
114 
115 template<typename T> class variable_if_dynamic<T, Dynamic>
116 {
117     T m_value;
118     EIGEN_DEVICE_FUNC variable_if_dynamic() { eigen_assert(false); }
119   public:
120     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit variable_if_dynamic(T value) : m_value(value) {}
121     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T value() const { return m_value; }
122     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void setValue(T value) { m_value = value; }
123 };
124 
125 /** \internal like variable_if_dynamic but for DynamicIndex
126   */
127 template<typename T, int Value> class variable_if_dynamicindex
128 {
129   public:
130     EIGEN_EMPTY_STRUCT_CTOR(variable_if_dynamicindex)
131     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit variable_if_dynamicindex(T v) { EIGEN_ONLY_USED_FOR_DEBUG(v); eigen_assert(v == T(Value)); }
132     EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE T value() { return T(Value); }
133     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void setValue(T) {}
134 };
135 
136 template<typename T> class variable_if_dynamicindex<T, DynamicIndex>
137 {
138     T m_value;
139     EIGEN_DEVICE_FUNC variable_if_dynamicindex() { eigen_assert(false); }
140   public:
141     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit variable_if_dynamicindex(T value) : m_value(value) {}
142     EIGEN_DEVICE_FUNC T EIGEN_STRONG_INLINE value() const { return m_value; }
143     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void setValue(T value) { m_value = value; }
144 };
145 
146 template<typename T> struct functor_traits
147 {
148   enum
149   {
150     Cost = 10,
151     PacketAccess = false,
152     IsRepeatable = false
153   };
154 };
155 
156 template<typename T> struct packet_traits;
157 
158 template<typename T> struct unpacket_traits
159 {
160   typedef T type;
161   typedef T half;
162   enum
163   {
164     size = 1,
165     alignment = 1
166   };
167 };
168 
169 template<int Size, typename PacketType,
170          bool Stop = Size==Dynamic || (Size%unpacket_traits<PacketType>::size)==0 || is_same<PacketType,typename unpacket_traits<PacketType>::half>::value>
171 struct find_best_packet_helper;
172 
173 template< int Size, typename PacketType>
174 struct find_best_packet_helper<Size,PacketType,true>
175 {
176   typedef PacketType type;
177 };
178 
179 template<int Size, typename PacketType>
180 struct find_best_packet_helper<Size,PacketType,false>
181 {
182   typedef typename find_best_packet_helper<Size,typename unpacket_traits<PacketType>::half>::type type;
183 };
184 
185 template<typename T, int Size>
186 struct find_best_packet
187 {
188   typedef typename find_best_packet_helper<Size,typename packet_traits<T>::type>::type type;
189 };
190 
191 #if EIGEN_MAX_STATIC_ALIGN_BYTES>0
192 template<int ArrayBytes, int AlignmentBytes,
193          bool Match     =  bool((ArrayBytes%AlignmentBytes)==0),
194          bool TryHalf   =  bool(EIGEN_MIN_ALIGN_BYTES<AlignmentBytes) >
195 struct compute_default_alignment_helper
196 {
197   enum { value = 0 };
198 };
199 
200 template<int ArrayBytes, int AlignmentBytes, bool TryHalf>
201 struct compute_default_alignment_helper<ArrayBytes, AlignmentBytes, true, TryHalf> // Match
202 {
203   enum { value = AlignmentBytes };
204 };
205 
206 template<int ArrayBytes, int AlignmentBytes>
207 struct compute_default_alignment_helper<ArrayBytes, AlignmentBytes, false, true> // Try-half
208 {
209   // current packet too large, try with an half-packet
210   enum { value = compute_default_alignment_helper<ArrayBytes, AlignmentBytes/2>::value };
211 };
212 #else
213 // If static alignment is disabled, no need to bother.
214 // This also avoids a division by zero in "bool Match =  bool((ArrayBytes%AlignmentBytes)==0)"
215 template<int ArrayBytes, int AlignmentBytes>
216 struct compute_default_alignment_helper
217 {
218   enum { value = 0 };
219 };
220 #endif
221 
222 template<typename T, int Size> struct compute_default_alignment {
223   enum { value = compute_default_alignment_helper<Size*sizeof(T),EIGEN_MAX_STATIC_ALIGN_BYTES>::value };
224 };
225 
226 template<typename T> struct compute_default_alignment<T,Dynamic> {
227   enum { value = EIGEN_MAX_ALIGN_BYTES };
228 };
229 
230 template<typename _Scalar, int _Rows, int _Cols,
231          int _Options = AutoAlign |
232                           ( (_Rows==1 && _Cols!=1) ? RowMajor
233                           : (_Cols==1 && _Rows!=1) ? ColMajor
234                           : EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION ),
235          int _MaxRows = _Rows,
236          int _MaxCols = _Cols
237 > class make_proper_matrix_type
238 {
239     enum {
240       IsColVector = _Cols==1 && _Rows!=1,
241       IsRowVector = _Rows==1 && _Cols!=1,
242       Options = IsColVector ? (_Options | ColMajor) & ~RowMajor
243               : IsRowVector ? (_Options | RowMajor) & ~ColMajor
244               : _Options
245     };
246   public:
247     typedef Matrix<_Scalar, _Rows, _Cols, Options, _MaxRows, _MaxCols> type;
248 };
249 
250 template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
251 class compute_matrix_flags
252 {
253     enum { row_major_bit = Options&RowMajor ? RowMajorBit : 0 };
254   public:
255     // FIXME currently we still have to handle DirectAccessBit at the expression level to handle DenseCoeffsBase<>
256     // and then propagate this information to the evaluator's flags.
257     // However, I (Gael) think that DirectAccessBit should only matter at the evaluation stage.
258     enum { ret = DirectAccessBit | LvalueBit | NestByRefBit | row_major_bit };
259 };
260 
261 template<int _Rows, int _Cols> struct size_at_compile_time
262 {
263   enum { ret = (_Rows==Dynamic || _Cols==Dynamic) ? Dynamic : _Rows * _Cols };
264 };
265 
266 template<typename XprType> struct size_of_xpr_at_compile_time
267 {
268   enum { ret = size_at_compile_time<traits<XprType>::RowsAtCompileTime,traits<XprType>::ColsAtCompileTime>::ret };
269 };
270 
271 /* plain_matrix_type : the difference from eval is that plain_matrix_type is always a plain matrix type,
272  * whereas eval is a const reference in the case of a matrix
273  */
274 
275 template<typename T, typename StorageKind = typename traits<T>::StorageKind> struct plain_matrix_type;
276 template<typename T, typename BaseClassType, int Flags> struct plain_matrix_type_dense;
277 template<typename T> struct plain_matrix_type<T,Dense>
278 {
279   typedef typename plain_matrix_type_dense<T,typename traits<T>::XprKind, traits<T>::Flags>::type type;
280 };
281 template<typename T> struct plain_matrix_type<T,DiagonalShape>
282 {
283   typedef typename T::PlainObject type;
284 };
285 
286 template<typename T, int Flags> struct plain_matrix_type_dense<T,MatrixXpr,Flags>
287 {
288   typedef Matrix<typename traits<T>::Scalar,
289                 traits<T>::RowsAtCompileTime,
290                 traits<T>::ColsAtCompileTime,
291                 AutoAlign | (Flags&RowMajorBit ? RowMajor : ColMajor),
292                 traits<T>::MaxRowsAtCompileTime,
293                 traits<T>::MaxColsAtCompileTime
294           > type;
295 };
296 
297 template<typename T, int Flags> struct plain_matrix_type_dense<T,ArrayXpr,Flags>
298 {
299   typedef Array<typename traits<T>::Scalar,
300                 traits<T>::RowsAtCompileTime,
301                 traits<T>::ColsAtCompileTime,
302                 AutoAlign | (Flags&RowMajorBit ? RowMajor : ColMajor),
303                 traits<T>::MaxRowsAtCompileTime,
304                 traits<T>::MaxColsAtCompileTime
305           > type;
306 };
307 
308 /* eval : the return type of eval(). For matrices, this is just a const reference
309  * in order to avoid a useless copy
310  */
311 
312 template<typename T, typename StorageKind = typename traits<T>::StorageKind> struct eval;
313 
314 template<typename T> struct eval<T,Dense>
315 {
316   typedef typename plain_matrix_type<T>::type type;
317 //   typedef typename T::PlainObject type;
318 //   typedef T::Matrix<typename traits<T>::Scalar,
319 //                 traits<T>::RowsAtCompileTime,
320 //                 traits<T>::ColsAtCompileTime,
321 //                 AutoAlign | (traits<T>::Flags&RowMajorBit ? RowMajor : ColMajor),
322 //                 traits<T>::MaxRowsAtCompileTime,
323 //                 traits<T>::MaxColsAtCompileTime
324 //           > type;
325 };
326 
327 template<typename T> struct eval<T,DiagonalShape>
328 {
329   typedef typename plain_matrix_type<T>::type type;
330 };
331 
332 // for matrices, no need to evaluate, just use a const reference to avoid a useless copy
333 template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
334 struct eval<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>, Dense>
335 {
336   typedef const Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& type;
337 };
338 
339 template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
340 struct eval<Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>, Dense>
341 {
342   typedef const Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& type;
343 };
344 
345 
346 /* similar to plain_matrix_type, but using the evaluator's Flags */
347 template<typename T, typename StorageKind = typename traits<T>::StorageKind> struct plain_object_eval;
348 
349 template<typename T>
350 struct plain_object_eval<T,Dense>
351 {
352   typedef typename plain_matrix_type_dense<T,typename traits<T>::XprKind, evaluator<T>::Flags>::type type;
353 };
354 
355 
356 /* plain_matrix_type_column_major : same as plain_matrix_type but guaranteed to be column-major
357  */
358 template<typename T> struct plain_matrix_type_column_major
359 {
360   enum { Rows = traits<T>::RowsAtCompileTime,
361          Cols = traits<T>::ColsAtCompileTime,
362          MaxRows = traits<T>::MaxRowsAtCompileTime,
363          MaxCols = traits<T>::MaxColsAtCompileTime
364   };
365   typedef Matrix<typename traits<T>::Scalar,
366                 Rows,
367                 Cols,
368                 (MaxRows==1&&MaxCols!=1) ? RowMajor : ColMajor,
369                 MaxRows,
370                 MaxCols
371           > type;
372 };
373 
374 /* plain_matrix_type_row_major : same as plain_matrix_type but guaranteed to be row-major
375  */
376 template<typename T> struct plain_matrix_type_row_major
377 {
378   enum { Rows = traits<T>::RowsAtCompileTime,
379          Cols = traits<T>::ColsAtCompileTime,
380          MaxRows = traits<T>::MaxRowsAtCompileTime,
381          MaxCols = traits<T>::MaxColsAtCompileTime
382   };
383   typedef Matrix<typename traits<T>::Scalar,
384                 Rows,
385                 Cols,
386                 (MaxCols==1&&MaxRows!=1) ? RowMajor : ColMajor,
387                 MaxRows,
388                 MaxCols
389           > type;
390 };
391 
392 /** \internal The reference selector for template expressions. The idea is that we don't
393   * need to use references for expressions since they are light weight proxy
394   * objects which should generate no copying overhead. */
395 template <typename T>
396 struct ref_selector
397 {
398   typedef typename conditional<
399     bool(traits<T>::Flags & NestByRefBit),
400     T const&,
401     const T
402   >::type type;
403 
404   typedef typename conditional<
405     bool(traits<T>::Flags & NestByRefBit),
406     T &,
407     T
408   >::type non_const_type;
409 };
410 
411 /** \internal Adds the const qualifier on the value-type of T2 if and only if T1 is a const type */
412 template<typename T1, typename T2>
413 struct transfer_constness
414 {
415   typedef typename conditional<
416     bool(internal::is_const<T1>::value),
417     typename internal::add_const_on_value_type<T2>::type,
418     T2
419   >::type type;
420 };
421 
422 
423 // However, we still need a mechanism to detect whether an expression which is evaluated multiple time
424 // has to be evaluated into a temporary.
425 // That's the purpose of this new nested_eval helper:
426 /** \internal Determines how a given expression should be nested when evaluated multiple times.
427   * For example, when you do a * (b+c), Eigen will determine how the expression b+c should be
428   * evaluated into the bigger product expression. The choice is between nesting the expression b+c as-is, or
429   * evaluating that expression b+c into a temporary variable d, and nest d so that the resulting expression is
430   * a*d. Evaluating can be beneficial for example if every coefficient access in the resulting expression causes
431   * many coefficient accesses in the nested expressions -- as is the case with matrix product for example.
432   *
433   * \tparam T the type of the expression being nested.
434   * \tparam n the number of coefficient accesses in the nested expression for each coefficient access in the bigger expression.
435   * \tparam PlainObject the type of the temporary if needed.
436   */
437 template<typename T, int n, typename PlainObject = typename plain_object_eval<T>::type> struct nested_eval
438 {
439   enum {
440     ScalarReadCost = NumTraits<typename traits<T>::Scalar>::ReadCost,
441     CoeffReadCost = evaluator<T>::CoeffReadCost,  // NOTE What if an evaluator evaluate itself into a tempory?
442                                                   //      Then CoeffReadCost will be small (e.g., 1) but we still have to evaluate, especially if n>1.
443                                                   //      This situation is already taken care by the EvalBeforeNestingBit flag, which is turned ON
444                                                   //      for all evaluator creating a temporary. This flag is then propagated by the parent evaluators.
445                                                   //      Another solution could be to count the number of temps?
446     NAsInteger = n == Dynamic ? HugeCost : n,
447     CostEval   = (NAsInteger+1) * ScalarReadCost + CoeffReadCost,
448     CostNoEval = NAsInteger * CoeffReadCost,
449     Evaluate = (int(evaluator<T>::Flags) & EvalBeforeNestingBit) || (int(CostEval) < int(CostNoEval))
450   };
451 
452   typedef typename conditional<Evaluate, PlainObject, typename ref_selector<T>::type>::type type;
453 };
454 
455 template<typename T>
456 EIGEN_DEVICE_FUNC
457 inline T* const_cast_ptr(const T* ptr)
458 {
459   return const_cast<T*>(ptr);
460 }
461 
462 template<typename Derived, typename XprKind = typename traits<Derived>::XprKind>
463 struct dense_xpr_base
464 {
465   /* dense_xpr_base should only ever be used on dense expressions, thus falling either into the MatrixXpr or into the ArrayXpr cases */
466 };
467 
468 template<typename Derived>
469 struct dense_xpr_base<Derived, MatrixXpr>
470 {
471   typedef MatrixBase<Derived> type;
472 };
473 
474 template<typename Derived>
475 struct dense_xpr_base<Derived, ArrayXpr>
476 {
477   typedef ArrayBase<Derived> type;
478 };
479 
480 template<typename Derived, typename XprKind = typename traits<Derived>::XprKind, typename StorageKind = typename traits<Derived>::StorageKind>
481 struct generic_xpr_base;
482 
483 template<typename Derived, typename XprKind>
484 struct generic_xpr_base<Derived, XprKind, Dense>
485 {
486   typedef typename dense_xpr_base<Derived,XprKind>::type type;
487 };
488 
489 template<typename XprType, typename CastType> struct cast_return_type
490 {
491   typedef typename XprType::Scalar CurrentScalarType;
492   typedef typename remove_all<CastType>::type _CastType;
493   typedef typename _CastType::Scalar NewScalarType;
494   typedef typename conditional<is_same<CurrentScalarType,NewScalarType>::value,
495                               const XprType&,CastType>::type type;
496 };
497 
498 template <typename A, typename B> struct promote_storage_type;
499 
500 template <typename A> struct promote_storage_type<A,A>
501 {
502   typedef A ret;
503 };
504 template <typename A> struct promote_storage_type<A, const A>
505 {
506   typedef A ret;
507 };
508 template <typename A> struct promote_storage_type<const A, A>
509 {
510   typedef A ret;
511 };
512 
513 /** \internal Specify the "storage kind" of applying a coefficient-wise
514   * binary operations between two expressions of kinds A and B respectively.
515   * The template parameter Functor permits to specialize the resulting storage kind wrt to
516   * the functor.
517   * The default rules are as follows:
518   * \code
519   * A      op A      -> A
520   * A      op dense  -> dense
521   * dense  op B      -> dense
522   * sparse op dense  -> sparse
523   * dense  op sparse -> sparse
524   * \endcode
525   */
526 template <typename A, typename B, typename Functor> struct cwise_promote_storage_type;
527 
528 template <typename A, typename Functor>                   struct cwise_promote_storage_type<A,A,Functor>                                      { typedef A      ret; };
529 template <typename Functor>                               struct cwise_promote_storage_type<Dense,Dense,Functor>                              { typedef Dense  ret; };
530 template <typename A, typename Functor>                   struct cwise_promote_storage_type<A,Dense,Functor>                                  { typedef Dense  ret; };
531 template <typename B, typename Functor>                   struct cwise_promote_storage_type<Dense,B,Functor>                                  { typedef Dense  ret; };
532 template <typename Functor>                               struct cwise_promote_storage_type<Sparse,Dense,Functor>                             { typedef Sparse ret; };
533 template <typename Functor>                               struct cwise_promote_storage_type<Dense,Sparse,Functor>                             { typedef Sparse ret; };
534 
535 template <typename LhsKind, typename RhsKind, int LhsOrder, int RhsOrder> struct cwise_promote_storage_order {
536   enum { value = LhsOrder };
537 };
538 
539 template <typename LhsKind, int LhsOrder, int RhsOrder>   struct cwise_promote_storage_order<LhsKind,Sparse,LhsOrder,RhsOrder>                { enum { value = RhsOrder }; };
540 template <typename RhsKind, int LhsOrder, int RhsOrder>   struct cwise_promote_storage_order<Sparse,RhsKind,LhsOrder,RhsOrder>                { enum { value = LhsOrder }; };
541 template <int Order>                                      struct cwise_promote_storage_order<Sparse,Sparse,Order,Order>                       { enum { value = Order }; };
542 
543 
544 /** \internal Specify the "storage kind" of multiplying an expression of kind A with kind B.
545   * The template parameter ProductTag permits to specialize the resulting storage kind wrt to
546   * some compile-time properties of the product: GemmProduct, GemvProduct, OuterProduct, InnerProduct.
547   * The default rules are as follows:
548   * \code
549   *  K * K            -> K
550   *  dense * K        -> dense
551   *  K * dense        -> dense
552   *  diag * K         -> K
553   *  K * diag         -> K
554   *  Perm * K         -> K
555   * K * Perm          -> K
556   * \endcode
557   */
558 template <typename A, typename B, int ProductTag> struct product_promote_storage_type;
559 
560 template <typename A, int ProductTag> struct product_promote_storage_type<A,                  A,                  ProductTag> { typedef A     ret;};
561 template <int ProductTag>             struct product_promote_storage_type<Dense,              Dense,              ProductTag> { typedef Dense ret;};
562 template <typename A, int ProductTag> struct product_promote_storage_type<A,                  Dense,              ProductTag> { typedef Dense ret; };
563 template <typename B, int ProductTag> struct product_promote_storage_type<Dense,              B,                  ProductTag> { typedef Dense ret; };
564 
565 template <typename A, int ProductTag> struct product_promote_storage_type<A,                  DiagonalShape,      ProductTag> { typedef A ret; };
566 template <typename B, int ProductTag> struct product_promote_storage_type<DiagonalShape,      B,                  ProductTag> { typedef B ret; };
567 template <int ProductTag>             struct product_promote_storage_type<Dense,              DiagonalShape,      ProductTag> { typedef Dense ret; };
568 template <int ProductTag>             struct product_promote_storage_type<DiagonalShape,      Dense,              ProductTag> { typedef Dense ret; };
569 
570 template <typename A, int ProductTag> struct product_promote_storage_type<A,                  PermutationStorage, ProductTag> { typedef A ret; };
571 template <typename B, int ProductTag> struct product_promote_storage_type<PermutationStorage, B,                  ProductTag> { typedef B ret; };
572 template <int ProductTag>             struct product_promote_storage_type<Dense,              PermutationStorage, ProductTag> { typedef Dense ret; };
573 template <int ProductTag>             struct product_promote_storage_type<PermutationStorage, Dense,              ProductTag> { typedef Dense ret; };
574 
575 /** \internal gives the plain matrix or array type to store a row/column/diagonal of a matrix type.
576   * \tparam Scalar optional parameter allowing to pass a different scalar type than the one of the MatrixType.
577   */
578 template<typename ExpressionType, typename Scalar = typename ExpressionType::Scalar>
579 struct plain_row_type
580 {
581   typedef Matrix<Scalar, 1, ExpressionType::ColsAtCompileTime,
582                  ExpressionType::PlainObject::Options | RowMajor, 1, ExpressionType::MaxColsAtCompileTime> MatrixRowType;
583   typedef Array<Scalar, 1, ExpressionType::ColsAtCompileTime,
584                  ExpressionType::PlainObject::Options | RowMajor, 1, ExpressionType::MaxColsAtCompileTime> ArrayRowType;
585 
586   typedef typename conditional<
587     is_same< typename traits<ExpressionType>::XprKind, MatrixXpr >::value,
588     MatrixRowType,
589     ArrayRowType
590   >::type type;
591 };
592 
593 template<typename ExpressionType, typename Scalar = typename ExpressionType::Scalar>
594 struct plain_col_type
595 {
596   typedef Matrix<Scalar, ExpressionType::RowsAtCompileTime, 1,
597                  ExpressionType::PlainObject::Options & ~RowMajor, ExpressionType::MaxRowsAtCompileTime, 1> MatrixColType;
598   typedef Array<Scalar, ExpressionType::RowsAtCompileTime, 1,
599                  ExpressionType::PlainObject::Options & ~RowMajor, ExpressionType::MaxRowsAtCompileTime, 1> ArrayColType;
600 
601   typedef typename conditional<
602     is_same< typename traits<ExpressionType>::XprKind, MatrixXpr >::value,
603     MatrixColType,
604     ArrayColType
605   >::type type;
606 };
607 
608 template<typename ExpressionType, typename Scalar = typename ExpressionType::Scalar>
609 struct plain_diag_type
610 {
611   enum { diag_size = EIGEN_SIZE_MIN_PREFER_DYNAMIC(ExpressionType::RowsAtCompileTime, ExpressionType::ColsAtCompileTime),
612          max_diag_size = EIGEN_SIZE_MIN_PREFER_FIXED(ExpressionType::MaxRowsAtCompileTime, ExpressionType::MaxColsAtCompileTime)
613   };
614   typedef Matrix<Scalar, diag_size, 1, ExpressionType::PlainObject::Options & ~RowMajor, max_diag_size, 1> MatrixDiagType;
615   typedef Array<Scalar, diag_size, 1, ExpressionType::PlainObject::Options & ~RowMajor, max_diag_size, 1> ArrayDiagType;
616 
617   typedef typename conditional<
618     is_same< typename traits<ExpressionType>::XprKind, MatrixXpr >::value,
619     MatrixDiagType,
620     ArrayDiagType
621   >::type type;
622 };
623 
624 template<typename Expr,typename Scalar = typename Expr::Scalar>
625 struct plain_constant_type
626 {
627   enum { Options = (traits<Expr>::Flags&RowMajorBit)?RowMajor:0 };
628 
629   typedef Array<Scalar,  traits<Expr>::RowsAtCompileTime,   traits<Expr>::ColsAtCompileTime,
630                 Options, traits<Expr>::MaxRowsAtCompileTime,traits<Expr>::MaxColsAtCompileTime> array_type;
631 
632   typedef Matrix<Scalar,  traits<Expr>::RowsAtCompileTime,   traits<Expr>::ColsAtCompileTime,
633                  Options, traits<Expr>::MaxRowsAtCompileTime,traits<Expr>::MaxColsAtCompileTime> matrix_type;
634 
635   typedef CwiseNullaryOp<scalar_constant_op<Scalar>, const typename conditional<is_same< typename traits<Expr>::XprKind, MatrixXpr >::value, matrix_type, array_type>::type > type;
636 };
637 
638 template<typename ExpressionType>
639 struct is_lvalue
640 {
641   enum { value = (!bool(is_const<ExpressionType>::value)) &&
642                  bool(traits<ExpressionType>::Flags & LvalueBit) };
643 };
644 
645 template<typename T> struct is_diagonal
646 { enum { ret = false }; };
647 
648 template<typename T> struct is_diagonal<DiagonalBase<T> >
649 { enum { ret = true }; };
650 
651 template<typename T> struct is_diagonal<DiagonalWrapper<T> >
652 { enum { ret = true }; };
653 
654 template<typename T, int S> struct is_diagonal<DiagonalMatrix<T,S> >
655 { enum { ret = true }; };
656 
657 template<typename S1, typename S2> struct glue_shapes;
658 template<> struct glue_shapes<DenseShape,TriangularShape> { typedef TriangularShape type;  };
659 
660 template<typename T1, typename T2>
661 bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if<has_direct_access<T1>::ret&&has_direct_access<T2>::ret, T1>::type * = 0)
662 {
663   return (mat1.data()==mat2.data()) && (mat1.innerStride()==mat2.innerStride()) && (mat1.outerStride()==mat2.outerStride());
664 }
665 
666 template<typename T1, typename T2>
667 bool is_same_dense(const T1 &, const T2 &, typename enable_if<!(has_direct_access<T1>::ret&&has_direct_access<T2>::ret), T1>::type * = 0)
668 {
669   return false;
670 }
671 
672 // Internal helper defining the cost of a scalar division for the type T.
673 // The default heuristic can be specialized for each scalar type and architecture.
674 template<typename T,bool Vectorized=false,typename EnaleIf = void>
675 struct scalar_div_cost {
676   enum { value = 8*NumTraits<T>::MulCost };
677 };
678 
679 template<typename T,bool Vectorized>
680 struct scalar_div_cost<std::complex<T>, Vectorized> {
681   enum { value = 2*scalar_div_cost<T>::value
682                + 6*NumTraits<T>::MulCost
683                + 3*NumTraits<T>::AddCost
684   };
685 };
686 
687 
688 template<bool Vectorized>
689 struct scalar_div_cost<signed long,Vectorized,typename conditional<sizeof(long)==8,void,false_type>::type> { enum { value = 24 }; };
690 template<bool Vectorized>
691 struct scalar_div_cost<unsigned long,Vectorized,typename conditional<sizeof(long)==8,void,false_type>::type> { enum { value = 21 }; };
692 
693 
694 #ifdef EIGEN_DEBUG_ASSIGN
695 std::string demangle_traversal(int t)
696 {
697   if(t==DefaultTraversal) return "DefaultTraversal";
698   if(t==LinearTraversal) return "LinearTraversal";
699   if(t==InnerVectorizedTraversal) return "InnerVectorizedTraversal";
700   if(t==LinearVectorizedTraversal) return "LinearVectorizedTraversal";
701   if(t==SliceVectorizedTraversal) return "SliceVectorizedTraversal";
702   return "?";
703 }
704 std::string demangle_unrolling(int t)
705 {
706   if(t==NoUnrolling) return "NoUnrolling";
707   if(t==InnerUnrolling) return "InnerUnrolling";
708   if(t==CompleteUnrolling) return "CompleteUnrolling";
709   return "?";
710 }
711 std::string demangle_flags(int f)
712 {
713   std::string res;
714   if(f&RowMajorBit)                 res += " | RowMajor";
715   if(f&PacketAccessBit)             res += " | Packet";
716   if(f&LinearAccessBit)             res += " | Linear";
717   if(f&LvalueBit)                   res += " | Lvalue";
718   if(f&DirectAccessBit)             res += " | Direct";
719   if(f&NestByRefBit)                res += " | NestByRef";
720   if(f&NoPreferredStorageOrderBit)  res += " | NoPreferredStorageOrderBit";
721 
722   return res;
723 }
724 #endif
725 
726 } // end namespace internal
727 
728 
729 /** \class ScalarBinaryOpTraits
730   * \ingroup Core_Module
731   *
732   * \brief Determines whether the given binary operation of two numeric types is allowed and what the scalar return type is.
733   *
734   * This class permits to control the scalar return type of any binary operation performed on two different scalar types through (partial) template specializations.
735   *
736   * For instance, let \c U1, \c U2 and \c U3 be three user defined scalar types for which most operations between instances of \c U1 and \c U2 returns an \c U3.
737   * You can let %Eigen knows that by defining:
738     \code
739     template<typename BinaryOp>
740     struct ScalarBinaryOpTraits<U1,U2,BinaryOp> { typedef U3 ReturnType;  };
741     template<typename BinaryOp>
742     struct ScalarBinaryOpTraits<U2,U1,BinaryOp> { typedef U3 ReturnType;  };
743     \endcode
744   * You can then explicitly disable some particular operations to get more explicit error messages:
745     \code
746     template<>
747     struct ScalarBinaryOpTraits<U1,U2,internal::scalar_max_op<U1,U2> > {};
748     \endcode
749   * Or customize the return type for individual operation:
750     \code
751     template<>
752     struct ScalarBinaryOpTraits<U1,U2,internal::scalar_sum_op<U1,U2> > { typedef U1 ReturnType; };
753     \endcode
754   *
755   * By default, the following generic combinations are supported:
756   <table class="manual">
757   <tr><th>ScalarA</th><th>ScalarB</th><th>BinaryOp</th><th>ReturnType</th><th>Note</th></tr>
758   <tr            ><td>\c T </td><td>\c T </td><td>\c * </td><td>\c T </td><td></td></tr>
759   <tr class="alt"><td>\c NumTraits<T>::Real </td><td>\c T </td><td>\c * </td><td>\c T </td><td>Only if \c NumTraits<T>::IsComplex </td></tr>
760   <tr            ><td>\c T </td><td>\c NumTraits<T>::Real </td><td>\c * </td><td>\c T </td><td>Only if \c NumTraits<T>::IsComplex </td></tr>
761   </table>
762   *
763   * \sa CwiseBinaryOp
764   */
765 template<typename ScalarA, typename ScalarB, typename BinaryOp=internal::scalar_product_op<ScalarA,ScalarB> >
766 struct ScalarBinaryOpTraits
767 #ifndef EIGEN_PARSED_BY_DOXYGEN
768   // for backward compatibility, use the hints given by the (deprecated) internal::scalar_product_traits class.
769   : internal::scalar_product_traits<ScalarA,ScalarB>
770 #endif // EIGEN_PARSED_BY_DOXYGEN
771 {};
772 
773 template<typename T, typename BinaryOp>
774 struct ScalarBinaryOpTraits<T,T,BinaryOp>
775 {
776   typedef T ReturnType;
777 };
778 
779 template <typename T, typename BinaryOp>
780 struct ScalarBinaryOpTraits<T, typename NumTraits<typename internal::enable_if<NumTraits<T>::IsComplex,T>::type>::Real, BinaryOp>
781 {
782   typedef T ReturnType;
783 };
784 template <typename T, typename BinaryOp>
785 struct ScalarBinaryOpTraits<typename NumTraits<typename internal::enable_if<NumTraits<T>::IsComplex,T>::type>::Real, T, BinaryOp>
786 {
787   typedef T ReturnType;
788 };
789 
790 // For Matrix * Permutation
791 template<typename T, typename BinaryOp>
792 struct ScalarBinaryOpTraits<T,void,BinaryOp>
793 {
794   typedef T ReturnType;
795 };
796 
797 // For Permutation * Matrix
798 template<typename T, typename BinaryOp>
799 struct ScalarBinaryOpTraits<void,T,BinaryOp>
800 {
801   typedef T ReturnType;
802 };
803 
804 // for Permutation*Permutation
805 template<typename BinaryOp>
806 struct ScalarBinaryOpTraits<void,void,BinaryOp>
807 {
808   typedef void ReturnType;
809 };
810 
811 // We require Lhs and Rhs to have "compatible" scalar types.
812 // It is tempting to always allow mixing different types but remember that this is often impossible in the vectorized paths.
813 // So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to
814 // add together a float matrix and a double matrix.
815 #define EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) \
816   EIGEN_STATIC_ASSERT((Eigen::internal::has_ReturnType<ScalarBinaryOpTraits<LHS, RHS,BINOP> >::value), \
817     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
818 
819 } // end namespace Eigen
820 
821 #endif // EIGEN_XPRHELPER_H
822