Eigen  3.2.9
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
35 template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value>
37 
38 enum {
39  Large = 2,
40  Small = 3
41 };
42 
43 namespace internal {
44 
45 template<int Rows, int Cols, int Depth> struct product_type_selector;
46 
47 template<int Size, int MaxSize> struct product_size_category
48 {
49  enum { is_large = MaxSize == Dynamic ||
50  Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD,
51  value = is_large ? Large
52  : Size == 1 ? 1
53  : Small
54  };
55 };
56 
57 template<typename Lhs, typename Rhs> struct product_type
58 {
59  typedef typename remove_all<Lhs>::type _Lhs;
60  typedef typename remove_all<Rhs>::type _Rhs;
61  enum {
62  MaxRows = _Lhs::MaxRowsAtCompileTime,
63  Rows = _Lhs::RowsAtCompileTime,
64  MaxCols = _Rhs::MaxColsAtCompileTime,
65  Cols = _Rhs::ColsAtCompileTime,
66  MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime,
67  _Rhs::MaxRowsAtCompileTime),
68  Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime,
69  _Rhs::RowsAtCompileTime),
70  LargeThreshold = EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
71  };
72 
73  // the splitting into different lines of code here, introducing the _select enums and the typedef below,
74  // is to work around an internal compiler error with gcc 4.1 and 4.2.
75 private:
76  enum {
77  rows_select = product_size_category<Rows,MaxRows>::value,
78  cols_select = product_size_category<Cols,MaxCols>::value,
79  depth_select = product_size_category<Depth,MaxDepth>::value
80  };
81  typedef product_type_selector<rows_select, cols_select, depth_select> selector;
82 
83 public:
84  enum {
85  value = selector::ret
86  };
87 #ifdef EIGEN_DEBUG_PRODUCT
88  static void debug()
89  {
90  EIGEN_DEBUG_VAR(Rows);
91  EIGEN_DEBUG_VAR(Cols);
92  EIGEN_DEBUG_VAR(Depth);
93  EIGEN_DEBUG_VAR(rows_select);
94  EIGEN_DEBUG_VAR(cols_select);
95  EIGEN_DEBUG_VAR(depth_select);
96  EIGEN_DEBUG_VAR(value);
97  }
98 #endif
99 };
100 
101 
102 /* The following allows to select the kind of product at compile time
103  * based on the three dimensions of the product.
104  * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
105 // FIXME I'm not sure the current mapping is the ideal one.
106 template<int M, int N> struct product_type_selector<M,N,1> { enum { ret = OuterProduct }; };
107 template<int Depth> struct product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; };
108 template<> struct product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; };
109 template<> struct product_type_selector<Small,1, Small> { enum { ret = CoeffBasedProductMode }; };
110 template<> struct product_type_selector<1, Small,Small> { enum { ret = CoeffBasedProductMode }; };
111 template<> struct product_type_selector<Small,Small,Small> { enum { ret = CoeffBasedProductMode }; };
112 template<> struct product_type_selector<Small, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; };
113 template<> struct product_type_selector<Small, Large, 1> { enum { ret = LazyCoeffBasedProductMode }; };
114 template<> struct product_type_selector<Large, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; };
115 template<> struct product_type_selector<1, Large,Small> { enum { ret = CoeffBasedProductMode }; };
116 template<> struct product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; };
117 template<> struct product_type_selector<1, Small,Large> { enum { ret = CoeffBasedProductMode }; };
118 template<> struct product_type_selector<Large,1, Small> { enum { ret = CoeffBasedProductMode }; };
119 template<> struct product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; };
120 template<> struct product_type_selector<Small,1, Large> { enum { ret = CoeffBasedProductMode }; };
121 template<> struct product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; };
122 template<> struct product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; };
123 template<> struct product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; };
124 template<> struct product_type_selector<Large,Large,Large> { enum { ret = GemmProduct }; };
125 template<> struct product_type_selector<Large,Small,Small> { enum { ret = GemmProduct }; };
126 template<> struct product_type_selector<Small,Large,Small> { enum { ret = GemmProduct }; };
127 template<> struct product_type_selector<Large,Large,Small> { enum { ret = GemmProduct }; };
128 
129 } // end namespace internal
130 
148 template<typename Lhs, typename Rhs, int ProductType>
150 {
151  // TODO use the nested type to reduce instanciations ????
152 // typedef typename internal::nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
153 // typedef typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
154 
155  typedef GeneralProduct<Lhs/*Nested*/, Rhs/*Nested*/, ProductType> Type;
156 };
157 
158 template<typename Lhs, typename Rhs>
159 struct ProductReturnType<Lhs,Rhs,CoeffBasedProductMode>
160 {
161  typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
162  typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
163  typedef CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit> Type;
164 };
165 
166 template<typename Lhs, typename Rhs>
167 struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
168 {
169  typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
170  typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
171  typedef CoeffBasedProduct<LhsNested, RhsNested, NestByRefBit> Type;
172 };
173 
174 // this is a workaround for sun CC
175 template<typename Lhs, typename Rhs>
176 struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
177 {};
178 
179 /***********************************************************************
180 * Implementation of Inner Vector Vector Product
181 ***********************************************************************/
182 
183 // FIXME : maybe the "inner product" could return a Scalar
184 // instead of a 1x1 matrix ??
185 // Pro: more natural for the user
186 // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
187 // product ends up to a row-vector times col-vector product... To tackle this use
188 // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
189 
190 namespace internal {
191 
192 template<typename Lhs, typename Rhs>
193 struct traits<GeneralProduct<Lhs,Rhs,InnerProduct> >
194  : traits<Matrix<typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> >
195 {};
196 
197 }
198 
199 template<typename Lhs, typename Rhs>
200 class GeneralProduct<Lhs, Rhs, InnerProduct>
201  : internal::no_assignment_operator,
202  public Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1>
203 {
204  typedef Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> Base;
205  public:
206  GeneralProduct(const Lhs& lhs, const Rhs& rhs)
207  {
208  Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
209  }
210 
212  operator const typename Base::Scalar() const {
213  return Base::coeff(0,0);
214  }
215 };
216 
217 /***********************************************************************
218 * Implementation of Outer Vector Vector Product
219 ***********************************************************************/
220 
221 namespace internal {
222 
223 // Column major
224 template<typename ProductType, typename Dest, typename Func>
225 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const false_type&)
226 {
227  typedef typename Dest::Index Index;
228  // FIXME make sure lhs is sequentially stored
229  // FIXME not very good if rhs is real and lhs complex while alpha is real too
230  const Index cols = dest.cols();
231  for (Index j=0; j<cols; ++j)
232  func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
233 }
234 
235 // Row major
236 template<typename ProductType, typename Dest, typename Func>
237 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const true_type&) {
238  typedef typename Dest::Index Index;
239  // FIXME make sure rhs is sequentially stored
240  // FIXME not very good if lhs is real and rhs complex while alpha is real too
241  const Index rows = dest.rows();
242  for (Index i=0; i<rows; ++i)
243  func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
244 }
245 
246 template<typename Lhs, typename Rhs>
247 struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> >
248  : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> >
249 {};
250 
251 }
252 
253 template<typename Lhs, typename Rhs>
254 class GeneralProduct<Lhs, Rhs, OuterProduct>
255  : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs>
256 {
257  template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
258 
259  public:
260  EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
261 
262  GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
263  {
264  }
265 
266  struct set { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } };
267  struct add { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
268  struct sub { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
269  struct adds {
270  Scalar m_scale;
271  adds(const Scalar& s) : m_scale(s) {}
272  template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const {
273  dst.const_cast_derived() += m_scale * src;
274  }
275  };
276 
277  template<typename Dest>
278  inline void evalTo(Dest& dest) const {
279  internal::outer_product_selector_run(*this, dest, set(), is_row_major<Dest>());
280  }
281 
282  template<typename Dest>
283  inline void addTo(Dest& dest) const {
284  internal::outer_product_selector_run(*this, dest, add(), is_row_major<Dest>());
285  }
286 
287  template<typename Dest>
288  inline void subTo(Dest& dest) const {
289  internal::outer_product_selector_run(*this, dest, sub(), is_row_major<Dest>());
290  }
291 
292  template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
293  {
294  internal::outer_product_selector_run(*this, dest, adds(alpha), is_row_major<Dest>());
295  }
296 };
297 
298 /***********************************************************************
299 * Implementation of General Matrix Vector Product
300 ***********************************************************************/
301 
302 /* According to the shape/flags of the matrix we have to distinghish 3 different cases:
303  * 1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
304  * 2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
305  * 3 - all other cases are handled using a simple loop along the outer-storage direction.
306  * Therefore we need a lower level meta selector.
307  * Furthermore, if the matrix is the rhs, then the product has to be transposed.
308  */
309 namespace internal {
310 
311 template<typename Lhs, typename Rhs>
312 struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> >
313  : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> >
314 {};
315 
316 template<int Side, int StorageOrder, bool BlasCompatible>
317 struct gemv_selector;
318 
319 } // end namespace internal
320 
321 template<typename Lhs, typename Rhs>
322 class GeneralProduct<Lhs, Rhs, GemvProduct>
323  : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs>
324 {
325  public:
326  EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
327 
328  typedef typename Lhs::Scalar LhsScalar;
329  typedef typename Rhs::Scalar RhsScalar;
330 
331  GeneralProduct(const Lhs& a_lhs, const Rhs& a_rhs) : Base(a_lhs,a_rhs)
332  {
333 // EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value),
334 // YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
335  }
336 
337  enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
338  typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType;
339 
340  template<typename Dest> void scaleAndAddTo(Dest& dst, const Scalar& alpha) const
341  {
342  eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols());
343  internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
344  bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha);
345  }
346 };
347 
348 namespace internal {
349 
350 // The vector is on the left => transposition
351 template<int StorageOrder, bool BlasCompatible>
352 struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible>
353 {
354  template<typename ProductType, typename Dest>
355  static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
356  {
357  Transpose<Dest> destT(dest);
358  enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
359  gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible>
360  ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct>
361  (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha);
362  }
363 };
364 
365 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if;
366 
367 template<typename Scalar,int Size,int MaxSize>
368 struct gemv_static_vector_if<Scalar,Size,MaxSize,false>
369 {
370  EIGEN_STRONG_INLINE Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; }
371 };
372 
373 template<typename Scalar,int Size>
374 struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
375 {
376  EIGEN_STRONG_INLINE Scalar* data() { return 0; }
377 };
378 
379 template<typename Scalar,int Size,int MaxSize>
380 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
381 {
382  #if EIGEN_ALIGN_STATICALLY
383  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
384  EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
385  #else
386  // Some architectures cannot align on the stack,
387  // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
388  enum {
389  ForceAlignment = internal::packet_traits<Scalar>::Vectorizable,
390  PacketSize = internal::packet_traits<Scalar>::size
391  };
392  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
393  EIGEN_STRONG_INLINE Scalar* data() {
394  return ForceAlignment
395  ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
396  : m_data.array;
397  }
398  #endif
399 };
400 
401 template<> struct gemv_selector<OnTheRight,ColMajor,true>
402 {
403  template<typename ProductType, typename Dest>
404  static inline void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
405  {
406  typedef typename ProductType::Index Index;
407  typedef typename ProductType::LhsScalar LhsScalar;
408  typedef typename ProductType::RhsScalar RhsScalar;
409  typedef typename ProductType::Scalar ResScalar;
410  typedef typename ProductType::RealScalar RealScalar;
411  typedef typename ProductType::ActualLhsType ActualLhsType;
412  typedef typename ProductType::ActualRhsType ActualRhsType;
413  typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
414  typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
415  typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
416 
417  ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs());
418  ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs());
419 
420  ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
421  * RhsBlasTraits::extractScalarFactor(prod.rhs());
422 
423  // make sure Dest is a compile-time vector type (bug 1166)
424  typedef typename conditional<Dest::IsVectorAtCompileTime, Dest, typename Dest::ColXpr>::type ActualDest;
425 
426  enum {
427  // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
428  // on, the other hand it is good for the cache to pack the vector anyways...
429  EvalToDestAtCompileTime = (ActualDest::InnerStrideAtCompileTime==1),
430  ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
431  MightCannotUseDest = (ActualDest::InnerStrideAtCompileTime!=1) || ComplexByReal
432  };
433 
434  gemv_static_vector_if<ResScalar,ActualDest::SizeAtCompileTime,ActualDest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
435 
436  bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
437  bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
438 
439  RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
440 
441  ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
442  evalToDest ? dest.data() : static_dest.data());
443 
444  if(!evalToDest)
445  {
446  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
447  int size = dest.size();
448  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
449  #endif
450  if(!alphaIsCompatible)
451  {
452  MappedDest(actualDestPtr, dest.size()).setZero();
453  compatibleAlpha = RhsScalar(1);
454  }
455  else
456  MappedDest(actualDestPtr, dest.size()) = dest;
457  }
458 
459  general_matrix_vector_product
460  <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
461  actualLhs.rows(), actualLhs.cols(),
462  actualLhs.data(), actualLhs.outerStride(),
463  actualRhs.data(), actualRhs.innerStride(),
464  actualDestPtr, 1,
465  compatibleAlpha);
466 
467  if (!evalToDest)
468  {
469  if(!alphaIsCompatible)
470  dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
471  else
472  dest = MappedDest(actualDestPtr, dest.size());
473  }
474  }
475 };
476 
477 template<> struct gemv_selector<OnTheRight,RowMajor,true>
478 {
479  template<typename ProductType, typename Dest>
480  static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
481  {
482  typedef typename ProductType::LhsScalar LhsScalar;
483  typedef typename ProductType::RhsScalar RhsScalar;
484  typedef typename ProductType::Scalar ResScalar;
485  typedef typename ProductType::Index Index;
486  typedef typename ProductType::ActualLhsType ActualLhsType;
487  typedef typename ProductType::ActualRhsType ActualRhsType;
488  typedef typename ProductType::_ActualRhsType _ActualRhsType;
489  typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
490  typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
491 
492  typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
493  typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
494 
495  ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
496  * RhsBlasTraits::extractScalarFactor(prod.rhs());
497 
498  enum {
499  // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
500  // on, the other hand it is good for the cache to pack the vector anyways...
501  DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1
502  };
503 
504  gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
505 
506  ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
507  DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
508 
509  if(!DirectlyUseRhs)
510  {
511  #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
512  int size = actualRhs.size();
513  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
514  #endif
515  Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
516  }
517 
518  general_matrix_vector_product
519  <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
520  actualLhs.rows(), actualLhs.cols(),
521  actualLhs.data(), actualLhs.outerStride(),
522  actualRhsPtr, 1,
523  dest.data(), dest.col(0).innerStride(), //NOTE if dest is not a vector at compile-time, then dest.innerStride() might be wrong. (bug 1166)
524  actualAlpha);
525  }
526 };
527 
528 template<> struct gemv_selector<OnTheRight,ColMajor,false>
529 {
530  template<typename ProductType, typename Dest>
531  static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
532  {
533  typedef typename Dest::Index Index;
534  // TODO makes sure dest is sequentially stored in memory, otherwise use a temp
535  const Index size = prod.rhs().rows();
536  for(Index k=0; k<size; ++k)
537  dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k);
538  }
539 };
540 
541 template<> struct gemv_selector<OnTheRight,RowMajor,false>
542 {
543  template<typename ProductType, typename Dest>
544  static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
545  {
546  typedef typename Dest::Index Index;
547  // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp
548  const Index rows = prod.rows();
549  for(Index i=0; i<rows; ++i)
550  dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum();
551  }
552 };
553 
554 } // end namespace internal
555 
556 /***************************************************************************
557 * Implementation of matrix base methods
558 ***************************************************************************/
559 
566 template<typename Derived>
567 template<typename OtherDerived>
568 inline const typename ProductReturnType<Derived, OtherDerived>::Type
570 {
571  // A note regarding the function declaration: In MSVC, this function will sometimes
572  // not be inlined since DenseStorage is an unwindable object for dynamic
573  // matrices and product types are holding a member to store the result.
574  // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
575  enum {
576  ProductIsValid = Derived::ColsAtCompileTime==Dynamic
577  || OtherDerived::RowsAtCompileTime==Dynamic
578  || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
579  AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
580  SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
581  };
582  // note to the lost user:
583  // * for a dot product use: v1.dot(v2)
584  // * for a coeff-wise product use: v1.cwiseProduct(v2)
585  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
586  INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
587  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
588  INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
589  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
590 #ifdef EIGEN_DEBUG_PRODUCT
591  internal::product_type<Derived,OtherDerived>::debug();
592 #endif
593  return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
594 }
595 
607 template<typename Derived>
608 template<typename OtherDerived>
611 {
612  enum {
613  ProductIsValid = Derived::ColsAtCompileTime==Dynamic
614  || OtherDerived::RowsAtCompileTime==Dynamic
615  || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
616  AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
617  SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
618  };
619  // note to the lost user:
620  // * for a dot product use: v1.dot(v2)
621  // * for a coeff-wise product use: v1.cwiseProduct(v2)
622  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
623  INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
624  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
625  INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
626  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
627 
628  return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
629 }
630 
631 } // end namespace Eigen
632 
633 #endif // EIGEN_PRODUCT_H
Expression of the product of two general matrices or vectors.
Definition: GeneralProduct.h:36
Definition: Constants.h:194
Definition: Constants.h:279
const int Dynamic
Definition: Constants.h:21
Definition: Constants.h:266
const ScalarMultipleReturnType operator*(const Scalar &scalar) const
Definition: MatrixBase.h:50
Definition: Constants.h:277
Definition: Constants.h:264
const LazyProductReturnType< Derived, OtherDerived >::Type lazyProduct(const MatrixBase< OtherDerived > &other) const
Definition: GeneralProduct.h:610
const unsigned int RowMajorBit
Definition: Constants.h:53
Helper class to get the correct and optimized returned type of operator*.
Definition: GeneralProduct.h:149
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48