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 /** \class GeneralProduct
17   * \ingroup Core_Module
18   *
19   * \brief Expression of the product of two general matrices or vectors
20   *
21   * \param LhsNested the type used to store the left-hand side
22   * \param RhsNested the type used to store the right-hand side
23   * \param ProductMode the type of the product
24   *
25   * This class represents an expression of the product of two general matrices.
26   * We call a general matrix, a dense matrix with full storage. For instance,
27   * This excludes triangular, selfadjoint, and sparse matrices.
28   * It is the return type of the operator* between general matrices. Its template
29   * arguments are determined automatically by ProductReturnType. Therefore,
30   * GeneralProduct should never be used direclty. To determine the result type of a
31   * function which involves a matrix product, use ProductReturnType::Type.
32   *
33   * \sa ProductReturnType, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
34   */
35 template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value>
36 class GeneralProduct;
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
debugproduct_type88   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 
131 /** \class ProductReturnType
132   * \ingroup Core_Module
133   *
134   * \brief Helper class to get the correct and optimized returned type of operator*
135   *
136   * \param Lhs the type of the left-hand side
137   * \param Rhs the type of the right-hand side
138   * \param ProductMode the type of the product (determined automatically by internal::product_mode)
139   *
140   * This class defines the typename Type representing the optimized product expression
141   * between two matrix expressions. In practice, using ProductReturnType<Lhs,Rhs>::Type
142   * is the recommended way to define the result type of a function returning an expression
143   * which involve a matrix product. The class Product should never be
144   * used directly.
145   *
146   * \sa class Product, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
147   */
148 template<typename Lhs, typename Rhs, int ProductType>
149 struct ProductReturnType
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       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
209         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
210 
211       Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
212     }
213 
214     /** Convertion to scalar */
215     operator const typename Base::Scalar() const {
216       return Base::coeff(0,0);
217     }
218 };
219 
220 /***********************************************************************
221 *  Implementation of Outer Vector Vector Product
222 ***********************************************************************/
223 
224 namespace internal {
225 
226 // Column major
227 template<typename ProductType, typename Dest, typename Func>
228 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const false_type&)
229 {
230   typedef typename Dest::Index Index;
231   // FIXME make sure lhs is sequentially stored
232   // FIXME not very good if rhs is real and lhs complex while alpha is real too
233   const Index cols = dest.cols();
234   for (Index j=0; j<cols; ++j)
235     func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
236 }
237 
238 // Row major
239 template<typename ProductType, typename Dest, typename Func>
240 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const true_type&) {
241   typedef typename Dest::Index Index;
242   // FIXME make sure rhs is sequentially stored
243   // FIXME not very good if lhs is real and rhs complex while alpha is real too
244   const Index rows = dest.rows();
245   for (Index i=0; i<rows; ++i)
246     func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
247 }
248 
249 template<typename Lhs, typename Rhs>
250 struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> >
251  : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> >
252 {};
253 
254 }
255 
256 template<typename Lhs, typename Rhs>
257 class GeneralProduct<Lhs, Rhs, OuterProduct>
258   : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs>
259 {
260     template<typename T> struct IsRowMajor : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
261 
262   public:
263     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
264 
265     GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
266     {
267       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
268         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
269     }
270 
271     struct set  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived()  = src; } };
272     struct add  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
273     struct sub  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
274     struct adds {
275       Scalar m_scale;
276       adds(const Scalar& s) : m_scale(s) {}
277       template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const {
278         dst.const_cast_derived() += m_scale * src;
279       }
280     };
281 
282     template<typename Dest>
283     inline void evalTo(Dest& dest) const {
284       internal::outer_product_selector_run(*this, dest, set(), IsRowMajor<Dest>());
285     }
286 
287     template<typename Dest>
288     inline void addTo(Dest& dest) const {
289       internal::outer_product_selector_run(*this, dest, add(), IsRowMajor<Dest>());
290     }
291 
292     template<typename Dest>
293     inline void subTo(Dest& dest) const {
294       internal::outer_product_selector_run(*this, dest, sub(), IsRowMajor<Dest>());
295     }
296 
297     template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
298     {
299       internal::outer_product_selector_run(*this, dest, adds(alpha), IsRowMajor<Dest>());
300     }
301 };
302 
303 /***********************************************************************
304 *  Implementation of General Matrix Vector Product
305 ***********************************************************************/
306 
307 /*  According to the shape/flags of the matrix we have to distinghish 3 different cases:
308  *   1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
309  *   2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
310  *   3 - all other cases are handled using a simple loop along the outer-storage direction.
311  *  Therefore we need a lower level meta selector.
312  *  Furthermore, if the matrix is the rhs, then the product has to be transposed.
313  */
314 namespace internal {
315 
316 template<typename Lhs, typename Rhs>
317 struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> >
318  : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> >
319 {};
320 
321 template<int Side, int StorageOrder, bool BlasCompatible>
322 struct gemv_selector;
323 
324 } // end namespace internal
325 
326 template<typename Lhs, typename Rhs>
327 class GeneralProduct<Lhs, Rhs, GemvProduct>
328   : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs>
329 {
330   public:
331     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
332 
333     typedef typename Lhs::Scalar LhsScalar;
334     typedef typename Rhs::Scalar RhsScalar;
335 
336     GeneralProduct(const Lhs& a_lhs, const Rhs& a_rhs) : Base(a_lhs,a_rhs)
337     {
338 //       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value),
339 //         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
340     }
341 
342     enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
343     typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType;
344 
345     template<typename Dest> void scaleAndAddTo(Dest& dst, const Scalar& alpha) const
346     {
347       eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols());
348       internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
349                        bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha);
350     }
351 };
352 
353 namespace internal {
354 
355 // The vector is on the left => transposition
356 template<int StorageOrder, bool BlasCompatible>
357 struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible>
358 {
359   template<typename ProductType, typename Dest>
360   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
361   {
362     Transpose<Dest> destT(dest);
363     enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
364     gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible>
365       ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct>
366         (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha);
367   }
368 };
369 
370 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if;
371 
372 template<typename Scalar,int Size,int MaxSize>
373 struct gemv_static_vector_if<Scalar,Size,MaxSize,false>
374 {
375   EIGEN_STRONG_INLINE  Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; }
376 };
377 
378 template<typename Scalar,int Size>
379 struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
380 {
381   EIGEN_STRONG_INLINE Scalar* data() { return 0; }
382 };
383 
384 template<typename Scalar,int Size,int MaxSize>
385 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
386 {
387   #if EIGEN_ALIGN_STATICALLY
388   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
389   EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
390   #else
391   // Some architectures cannot align on the stack,
392   // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
393   enum {
394     ForceAlignment  = internal::packet_traits<Scalar>::Vectorizable,
395     PacketSize      = internal::packet_traits<Scalar>::size
396   };
397   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
398   EIGEN_STRONG_INLINE Scalar* data() {
399     return ForceAlignment
400             ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
401             : m_data.array;
402   }
403   #endif
404 };
405 
406 template<> struct gemv_selector<OnTheRight,ColMajor,true>
407 {
408   template<typename ProductType, typename Dest>
409   static inline void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
410   {
411     typedef typename ProductType::Index Index;
412     typedef typename ProductType::LhsScalar   LhsScalar;
413     typedef typename ProductType::RhsScalar   RhsScalar;
414     typedef typename ProductType::Scalar      ResScalar;
415     typedef typename ProductType::RealScalar  RealScalar;
416     typedef typename ProductType::ActualLhsType ActualLhsType;
417     typedef typename ProductType::ActualRhsType ActualRhsType;
418     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
419     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
420     typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
421 
422     ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs());
423     ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs());
424 
425     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
426                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
427 
428     enum {
429       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
430       // on, the other hand it is good for the cache to pack the vector anyways...
431       EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
432       ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
433       MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
434     };
435 
436     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
437 
438     bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
439     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
440 
441     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
442 
443     ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
444                                                   evalToDest ? dest.data() : static_dest.data());
445 
446     if(!evalToDest)
447     {
448       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
449       int size = dest.size();
450       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
451       #endif
452       if(!alphaIsCompatible)
453       {
454         MappedDest(actualDestPtr, dest.size()).setZero();
455         compatibleAlpha = RhsScalar(1);
456       }
457       else
458         MappedDest(actualDestPtr, dest.size()) = dest;
459     }
460 
461     general_matrix_vector_product
462       <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
463         actualLhs.rows(), actualLhs.cols(),
464         actualLhs.data(), actualLhs.outerStride(),
465         actualRhs.data(), actualRhs.innerStride(),
466         actualDestPtr, 1,
467         compatibleAlpha);
468 
469     if (!evalToDest)
470     {
471       if(!alphaIsCompatible)
472         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
473       else
474         dest = MappedDest(actualDestPtr, dest.size());
475     }
476   }
477 };
478 
479 template<> struct gemv_selector<OnTheRight,RowMajor,true>
480 {
481   template<typename ProductType, typename Dest>
482   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
483   {
484     typedef typename ProductType::LhsScalar LhsScalar;
485     typedef typename ProductType::RhsScalar RhsScalar;
486     typedef typename ProductType::Scalar    ResScalar;
487     typedef typename ProductType::Index Index;
488     typedef typename ProductType::ActualLhsType ActualLhsType;
489     typedef typename ProductType::ActualRhsType ActualRhsType;
490     typedef typename ProductType::_ActualRhsType _ActualRhsType;
491     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
492     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
493 
494     typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
495     typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
496 
497     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
498                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
499 
500     enum {
501       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
502       // on, the other hand it is good for the cache to pack the vector anyways...
503       DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1
504     };
505 
506     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
507 
508     ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
509         DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
510 
511     if(!DirectlyUseRhs)
512     {
513       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
514       int size = actualRhs.size();
515       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
516       #endif
517       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
518     }
519 
520     general_matrix_vector_product
521       <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
522         actualLhs.rows(), actualLhs.cols(),
523         actualLhs.data(), actualLhs.outerStride(),
524         actualRhsPtr, 1,
525         dest.data(), dest.innerStride(),
526         actualAlpha);
527   }
528 };
529 
530 template<> struct gemv_selector<OnTheRight,ColMajor,false>
531 {
532   template<typename ProductType, typename Dest>
533   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
534   {
535     typedef typename Dest::Index Index;
536     // TODO makes sure dest is sequentially stored in memory, otherwise use a temp
537     const Index size = prod.rhs().rows();
538     for(Index k=0; k<size; ++k)
539       dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k);
540   }
541 };
542 
543 template<> struct gemv_selector<OnTheRight,RowMajor,false>
544 {
545   template<typename ProductType, typename Dest>
546   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
547   {
548     typedef typename Dest::Index Index;
549     // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp
550     const Index rows = prod.rows();
551     for(Index i=0; i<rows; ++i)
552       dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum();
553   }
554 };
555 
556 } // end namespace internal
557 
558 /***************************************************************************
559 * Implementation of matrix base methods
560 ***************************************************************************/
561 
562 /** \returns the matrix product of \c *this and \a other.
563   *
564   * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*().
565   *
566   * \sa lazyProduct(), operator*=(const MatrixBase&), Cwise::operator*()
567   */
568 template<typename Derived>
569 template<typename OtherDerived>
570 inline const typename ProductReturnType<Derived, OtherDerived>::Type
571 MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
572 {
573   // A note regarding the function declaration: In MSVC, this function will sometimes
574   // not be inlined since DenseStorage is an unwindable object for dynamic
575   // matrices and product types are holding a member to store the result.
576   // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
577   enum {
578     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
579                    || OtherDerived::RowsAtCompileTime==Dynamic
580                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
581     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
582     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
583   };
584   // note to the lost user:
585   //    * for a dot product use: v1.dot(v2)
586   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
587   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
588     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
589   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
590     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
591   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
592 #ifdef EIGEN_DEBUG_PRODUCT
593   internal::product_type<Derived,OtherDerived>::debug();
594 #endif
595   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
596 }
597 
598 /** \returns an expression of the matrix product of \c *this and \a other without implicit evaluation.
599   *
600   * The returned product will behave like any other expressions: the coefficients of the product will be
601   * computed once at a time as requested. This might be useful in some extremely rare cases when only
602   * a small and no coherent fraction of the result's coefficients have to be computed.
603   *
604   * \warning This version of the matrix product can be much much slower. So use it only if you know
605   * what you are doing and that you measured a true speed improvement.
606   *
607   * \sa operator*(const MatrixBase&)
608   */
609 template<typename Derived>
610 template<typename OtherDerived>
611 const typename LazyProductReturnType<Derived,OtherDerived>::Type
612 MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const
613 {
614   enum {
615     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
616                    || OtherDerived::RowsAtCompileTime==Dynamic
617                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
618     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
619     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
620   };
621   // note to the lost user:
622   //    * for a dot product use: v1.dot(v2)
623   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
624   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
625     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
626   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
627     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
628   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
629 
630   return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
631 }
632 
633 } // end namespace Eigen
634 
635 #endif // EIGEN_PRODUCT_H
636