Eigen  3.2.92
GeneralProduct.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-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
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_GENERAL_PRODUCT_H
12 #define EIGEN_GENERAL_PRODUCT_H
13 
14 namespace Eigen {
15 
16 enum {
17  Large = 2,
18  Small = 3
19 };
20 
21 namespace internal {
22 
23 template<int Rows, int Cols, int Depth> struct product_type_selector;
24 
25 template<int Size, int MaxSize> struct product_size_category
26 {
27  enum { is_large = MaxSize == Dynamic ||
28  Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD,
29  value = is_large ? Large
30  : Size == 1 ? 1
31  : Small
32  };
33 };
34 
35 template<typename Lhs, typename Rhs> struct product_type
36 {
37  typedef typename remove_all<Lhs>::type _Lhs;
38  typedef typename remove_all<Rhs>::type _Rhs;
39  enum {
40  MaxRows = traits<_Lhs>::MaxRowsAtCompileTime,
41  Rows = traits<_Lhs>::RowsAtCompileTime,
42  MaxCols = traits<_Rhs>::MaxColsAtCompileTime,
43  Cols = traits<_Rhs>::ColsAtCompileTime,
44  MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(traits<_Lhs>::MaxColsAtCompileTime,
45  traits<_Rhs>::MaxRowsAtCompileTime),
46  Depth = EIGEN_SIZE_MIN_PREFER_FIXED(traits<_Lhs>::ColsAtCompileTime,
47  traits<_Rhs>::RowsAtCompileTime)
48  };
49 
50  // the splitting into different lines of code here, introducing the _select enums and the typedef below,
51  // is to work around an internal compiler error with gcc 4.1 and 4.2.
52 private:
53  enum {
54  rows_select = product_size_category<Rows,MaxRows>::value,
55  cols_select = product_size_category<Cols,MaxCols>::value,
56  depth_select = product_size_category<Depth,MaxDepth>::value
57  };
58  typedef product_type_selector<rows_select, cols_select, depth_select> selector;
59 
60 public:
61  enum {
62  value = selector::ret,
63  ret = selector::ret
64  };
65 #ifdef EIGEN_DEBUG_PRODUCT
66  static void debug()
67  {
68  EIGEN_DEBUG_VAR(Rows);
69  EIGEN_DEBUG_VAR(Cols);
70  EIGEN_DEBUG_VAR(Depth);
71  EIGEN_DEBUG_VAR(rows_select);
72  EIGEN_DEBUG_VAR(cols_select);
73  EIGEN_DEBUG_VAR(depth_select);
74  EIGEN_DEBUG_VAR(value);
75  }
76 #endif
77 };
78 
79 // template<typename Lhs, typename Rhs> struct product_tag
80 // {
81 // private:
82 //
83 // typedef typename remove_all<Lhs>::type _Lhs;
84 // typedef typename remove_all<Rhs>::type _Rhs;
85 // enum {
86 // Rows = _Lhs::RowsAtCompileTime,
87 // Cols = _Rhs::ColsAtCompileTime,
88 // Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime, _Rhs::RowsAtCompileTime)
89 // };
90 //
91 // enum {
92 // rows_select = Rows==1 ? int(Rows) : int(Large),
93 // cols_select = Cols==1 ? int(Cols) : int(Large),
94 // depth_select = Depth==1 ? int(Depth) : int(Large)
95 // };
96 // typedef product_type_selector<rows_select, cols_select, depth_select> selector;
97 //
98 // public:
99 // enum {
100 // ret = selector::ret
101 // };
102 //
103 // };
104 
105 /* The following allows to select the kind of product at compile time
106  * based on the three dimensions of the product.
107  * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
108 // FIXME I'm not sure the current mapping is the ideal one.
109 template<int M, int N> struct product_type_selector<M,N,1> { enum { ret = OuterProduct }; };
110 template<int Depth> struct product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; };
111 template<> struct product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; };
112 template<> struct product_type_selector<Small,1, Small> { enum { ret = CoeffBasedProductMode }; };
113 template<> struct product_type_selector<1, Small,Small> { enum { ret = CoeffBasedProductMode }; };
114 template<> struct product_type_selector<Small,Small,Small> { enum { ret = CoeffBasedProductMode }; };
115 template<> struct product_type_selector<Small, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; };
116 template<> struct product_type_selector<Small, Large, 1> { enum { ret = LazyCoeffBasedProductMode }; };
117 template<> struct product_type_selector<Large, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; };
118 template<> struct product_type_selector<1, Large,Small> { enum { ret = CoeffBasedProductMode }; };
119 template<> struct product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; };
120 template<> struct product_type_selector<1, Small,Large> { enum { ret = CoeffBasedProductMode }; };
121 template<> struct product_type_selector<Large,1, Small> { enum { ret = CoeffBasedProductMode }; };
122 template<> struct product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; };
123 template<> struct product_type_selector<Small,1, Large> { enum { ret = CoeffBasedProductMode }; };
124 template<> struct product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; };
125 template<> struct product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; };
126 template<> struct product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; };
127 template<> struct product_type_selector<Large,Large,Large> { enum { ret = GemmProduct }; };
128 template<> struct product_type_selector<Large,Small,Small> { enum { ret = GemmProduct }; };
129 template<> struct product_type_selector<Small,Large,Small> { enum { ret = GemmProduct }; };
130 template<> struct product_type_selector<Large,Large,Small> { enum { ret = GemmProduct }; };
131 
132 } // end namespace internal
133 
134 /***********************************************************************
135 * Implementation of Inner Vector Vector Product
136 ***********************************************************************/
137 
138 // FIXME : maybe the "inner product" could return a Scalar
139 // instead of a 1x1 matrix ??
140 // Pro: more natural for the user
141 // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
142 // product ends up to a row-vector times col-vector product... To tackle this use
143 // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
144 
145 /***********************************************************************
146 * Implementation of Outer Vector Vector Product
147 ***********************************************************************/
148 
149 /***********************************************************************
150 * Implementation of General Matrix Vector Product
151 ***********************************************************************/
152 
153 /* According to the shape/flags of the matrix we have to distinghish 3 different cases:
154  * 1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
155  * 2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
156  * 3 - all other cases are handled using a simple loop along the outer-storage direction.
157  * Therefore we need a lower level meta selector.
158  * Furthermore, if the matrix is the rhs, then the product has to be transposed.
159  */
160 namespace internal {
161 
162 template<int Side, int StorageOrder, bool BlasCompatible>
163 struct gemv_dense_selector;
164 
165 } // end namespace internal
166 
167 namespace internal {
168 
169 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if;
170 
171 template<typename Scalar,int Size,int MaxSize>
172 struct gemv_static_vector_if<Scalar,Size,MaxSize,false>
173 {
174  EIGEN_STRONG_INLINE Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; }
175 };
176 
177 template<typename Scalar,int Size>
178 struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
179 {
180  EIGEN_STRONG_INLINE Scalar* data() { return 0; }
181 };
182 
183 template<typename Scalar,int Size,int MaxSize>
184 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
185 {
186  #if EIGEN_MAX_STATIC_ALIGN_BYTES!=0
187  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
188  EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
189  #else
190  // Some architectures cannot align on the stack,
191  // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
192  enum {
193  ForceAlignment = internal::packet_traits<Scalar>::Vectorizable,
194  PacketSize = internal::packet_traits<Scalar>::size
195  };
196  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
197  EIGEN_STRONG_INLINE Scalar* data() {
198  return ForceAlignment
199  ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(EIGEN_MAX_ALIGN_BYTES-1))) + EIGEN_MAX_ALIGN_BYTES)
200  : m_data.array;
201  }
202  #endif
203 };
204 
205 // The vector is on the left => transposition
206 template<int StorageOrder, bool BlasCompatible>
207 struct gemv_dense_selector<OnTheLeft,StorageOrder,BlasCompatible>
208 {
209  template<typename Lhs, typename Rhs, typename Dest>
210  static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha)
211  {
212  Transpose<Dest> destT(dest);
213  enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
214  gemv_dense_selector<OnTheRight,OtherStorageOrder,BlasCompatible>
215  ::run(rhs.transpose(), lhs.transpose(), destT, alpha);
216  }
217 };
218 
219 template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
220 {
221  template<typename Lhs, typename Rhs, typename Dest>
222  static inline void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha)
223  {
224  typedef typename Lhs::Scalar LhsScalar;
225  typedef typename Rhs::Scalar RhsScalar;
226  typedef typename Dest::Scalar ResScalar;
227  typedef typename Dest::RealScalar RealScalar;
228 
229  typedef internal::blas_traits<Lhs> LhsBlasTraits;
230  typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
231  typedef internal::blas_traits<Rhs> RhsBlasTraits;
232  typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
233 
234  typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
235 
236  ActualLhsType actualLhs = LhsBlasTraits::extract(lhs);
237  ActualRhsType actualRhs = RhsBlasTraits::extract(rhs);
238 
239  ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs)
240  * RhsBlasTraits::extractScalarFactor(rhs);
241 
242  enum {
243  // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
244  // on, the other hand it is good for the cache to pack the vector anyways...
245  EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
246  ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
247  MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
248  };
249 
250  gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
251 
252  const bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
253  const bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
254 
255  RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
256 
257  ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
258  evalToDest ? dest.data() : static_dest.data());
259 
260  if(!evalToDest)
261  {
262  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
263  Index size = dest.size();
264  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
265  #endif
266  if(!alphaIsCompatible)
267  {
268  MappedDest(actualDestPtr, dest.size()).setZero();
269  compatibleAlpha = RhsScalar(1);
270  }
271  else
272  MappedDest(actualDestPtr, dest.size()) = dest;
273  }
274 
275  typedef const_blas_data_mapper<LhsScalar,Index,ColMajor> LhsMapper;
276  typedef const_blas_data_mapper<RhsScalar,Index,RowMajor> RhsMapper;
277  general_matrix_vector_product
278  <Index,LhsScalar,LhsMapper,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsMapper,RhsBlasTraits::NeedToConjugate>::run(
279  actualLhs.rows(), actualLhs.cols(),
280  LhsMapper(actualLhs.data(), actualLhs.outerStride()),
281  RhsMapper(actualRhs.data(), actualRhs.innerStride()),
282  actualDestPtr, 1,
283  compatibleAlpha);
284 
285  if (!evalToDest)
286  {
287  if(!alphaIsCompatible)
288  dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
289  else
290  dest = MappedDest(actualDestPtr, dest.size());
291  }
292  }
293 };
294 
295 template<> struct gemv_dense_selector<OnTheRight,RowMajor,true>
296 {
297  template<typename Lhs, typename Rhs, typename Dest>
298  static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha)
299  {
300  typedef typename Lhs::Scalar LhsScalar;
301  typedef typename Rhs::Scalar RhsScalar;
302  typedef typename Dest::Scalar ResScalar;
303 
304  typedef internal::blas_traits<Lhs> LhsBlasTraits;
305  typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
306  typedef internal::blas_traits<Rhs> RhsBlasTraits;
307  typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
308  typedef typename internal::remove_all<ActualRhsType>::type ActualRhsTypeCleaned;
309 
310  typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(lhs);
311  typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(rhs);
312 
313  ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs)
314  * RhsBlasTraits::extractScalarFactor(rhs);
315 
316  enum {
317  // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
318  // on, the other hand it is good for the cache to pack the vector anyways...
319  DirectlyUseRhs = ActualRhsTypeCleaned::InnerStrideAtCompileTime==1
320  };
321 
322  gemv_static_vector_if<RhsScalar,ActualRhsTypeCleaned::SizeAtCompileTime,ActualRhsTypeCleaned::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
323 
324  ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
325  DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
326 
327  if(!DirectlyUseRhs)
328  {
329  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
330  Index size = actualRhs.size();
331  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
332  #endif
333  Map<typename ActualRhsTypeCleaned::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
334  }
335 
336  typedef const_blas_data_mapper<LhsScalar,Index,RowMajor> LhsMapper;
337  typedef const_blas_data_mapper<RhsScalar,Index,ColMajor> RhsMapper;
338  general_matrix_vector_product
339  <Index,LhsScalar,LhsMapper,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsMapper,RhsBlasTraits::NeedToConjugate>::run(
340  actualLhs.rows(), actualLhs.cols(),
341  LhsMapper(actualLhs.data(), actualLhs.outerStride()),
342  RhsMapper(actualRhsPtr, 1),
343  dest.data(), dest.innerStride(),
344  actualAlpha);
345  }
346 };
347 
348 template<> struct gemv_dense_selector<OnTheRight,ColMajor,false>
349 {
350  template<typename Lhs, typename Rhs, typename Dest>
351  static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha)
352  {
353  // TODO if rhs is large enough it might be beneficial to make sure that dest is sequentially stored in memory, otherwise use a temp
354  typename nested_eval<Rhs,1>::type actual_rhs(rhs);
355  const Index size = rhs.rows();
356  for(Index k=0; k<size; ++k)
357  dest += (alpha*actual_rhs.coeff(k)) * lhs.col(k);
358  }
359 };
360 
361 template<> struct gemv_dense_selector<OnTheRight,RowMajor,false>
362 {
363  template<typename Lhs, typename Rhs, typename Dest>
364  static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha)
365  {
366  typename nested_eval<Rhs,Lhs::RowsAtCompileTime>::type actual_rhs(rhs);
367  const Index rows = dest.rows();
368  for(Index i=0; i<rows; ++i)
369  dest.coeffRef(i) += alpha * (lhs.row(i).cwiseProduct(actual_rhs.transpose())).sum();
370  }
371 };
372 
373 } // end namespace internal
374 
375 /***************************************************************************
376 * Implementation of matrix base methods
377 ***************************************************************************/
378 
385 #ifndef __CUDACC__
386 
387 template<typename Derived>
388 template<typename OtherDerived>
389 inline const Product<Derived, OtherDerived>
391 {
392  // A note regarding the function declaration: In MSVC, this function will sometimes
393  // not be inlined since DenseStorage is an unwindable object for dynamic
394  // matrices and product types are holding a member to store the result.
395  // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
396  enum {
397  ProductIsValid = Derived::ColsAtCompileTime==Dynamic
398  || OtherDerived::RowsAtCompileTime==Dynamic
399  || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
400  AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
401  SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
402  };
403  // note to the lost user:
404  // * for a dot product use: v1.dot(v2)
405  // * for a coeff-wise product use: v1.cwiseProduct(v2)
406  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
407  INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
408  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
409  INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
410  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
411 #ifdef EIGEN_DEBUG_PRODUCT
412  internal::product_type<Derived,OtherDerived>::debug();
413 #endif
414 
415  return Product<Derived, OtherDerived>(derived(), other.derived());
416 }
417 
418 #endif // __CUDACC__
419 
431 template<typename Derived>
432 template<typename OtherDerived>
435 {
436  enum {
437  ProductIsValid = Derived::ColsAtCompileTime==Dynamic
438  || OtherDerived::RowsAtCompileTime==Dynamic
439  || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
440  AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
441  SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
442  };
443  // note to the lost user:
444  // * for a dot product use: v1.dot(v2)
445  // * for a coeff-wise product use: v1.cwiseProduct(v2)
446  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
447  INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
448  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
449  INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
450  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
451 
452  return Product<Derived,OtherDerived,LazyProduct>(derived(), other.derived());
453 }
454 
455 } // end namespace Eigen
456 
457 #endif // EIGEN_PRODUCT_H
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:107
Definition: Constants.h:333
const Product< Derived, OtherDerived, LazyProduct > lazyProduct(const MatrixBase< OtherDerived > &other) const
Definition: GeneralProduct.h:434
Definition: LDLT.h:16
Definition: Constants.h:320
Definition: Constants.h:335
Definition: Constants.h:322
const ScalarMultipleReturnType operator*(const Scalar &scalar) const
Definition: MatrixBase.h:57
Definition: Eigen_Colamd.h:54
Definition: Constants.h:235
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48