GeneralProduct.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_GENERAL_PRODUCT_H
00012 #define EIGEN_GENERAL_PRODUCT_H
00013 
00014 namespace Eigen { 
00015 
00016 /** \class GeneralProduct
00017   * \ingroup Core_Module
00018   *
00019   * \brief Expression of the product of two general matrices or vectors
00020   *
00021   * \param LhsNested the type used to store the left-hand side
00022   * \param RhsNested the type used to store the right-hand side
00023   * \param ProductMode the type of the product
00024   *
00025   * This class represents an expression of the product of two general matrices.
00026   * We call a general matrix, a dense matrix with full storage. For instance,
00027   * This excludes triangular, selfadjoint, and sparse matrices.
00028   * It is the return type of the operator* between general matrices. Its template
00029   * arguments are determined automatically by ProductReturnType. Therefore,
00030   * GeneralProduct should never be used direclty. To determine the result type of a
00031   * function which involves a matrix product, use ProductReturnType::Type.
00032   *
00033   * \sa ProductReturnType, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
00034   */
00035 template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value>
00036 class GeneralProduct;
00037 
00038 enum {
00039   Large = 2,
00040   Small = 3
00041 };
00042 
00043 namespace internal {
00044 
00045 template<int Rows, int Cols, int Depth> struct product_type_selector;
00046 
00047 template<int Size, int MaxSize> struct product_size_category
00048 {
00049   enum { is_large = MaxSize == Dynamic ||
00050                     Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD,
00051          value = is_large  ? Large
00052                : Size == 1 ? 1
00053                            : Small
00054   };
00055 };
00056 
00057 template<typename Lhs, typename Rhs> struct product_type
00058 {
00059   typedef typename remove_all<Lhs>::type _Lhs;
00060   typedef typename remove_all<Rhs>::type _Rhs;
00061   enum {
00062     MaxRows  = _Lhs::MaxRowsAtCompileTime,
00063     Rows  = _Lhs::RowsAtCompileTime,
00064     MaxCols  = _Rhs::MaxColsAtCompileTime,
00065     Cols  = _Rhs::ColsAtCompileTime,
00066     MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime,
00067                                            _Rhs::MaxRowsAtCompileTime),
00068     Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime,
00069                                         _Rhs::RowsAtCompileTime),
00070     LargeThreshold = EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
00071   };
00072 
00073   // the splitting into different lines of code here, introducing the _select enums and the typedef below,
00074   // is to work around an internal compiler error with gcc 4.1 and 4.2.
00075 private:
00076   enum {
00077     rows_select = product_size_category<Rows,MaxRows>::value,
00078     cols_select = product_size_category<Cols,MaxCols>::value,
00079     depth_select = product_size_category<Depth,MaxDepth>::value
00080   };
00081   typedef product_type_selector<rows_select, cols_select, depth_select> selector;
00082 
00083 public:
00084   enum {
00085     value = selector::ret
00086   };
00087 #ifdef EIGEN_DEBUG_PRODUCT
00088   static void debug()
00089   {
00090       EIGEN_DEBUG_VAR(Rows);
00091       EIGEN_DEBUG_VAR(Cols);
00092       EIGEN_DEBUG_VAR(Depth);
00093       EIGEN_DEBUG_VAR(rows_select);
00094       EIGEN_DEBUG_VAR(cols_select);
00095       EIGEN_DEBUG_VAR(depth_select);
00096       EIGEN_DEBUG_VAR(value);
00097   }
00098 #endif
00099 };
00100 
00101 
00102 /* The following allows to select the kind of product at compile time
00103  * based on the three dimensions of the product.
00104  * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
00105 // FIXME I'm not sure the current mapping is the ideal one.
00106 template<int M, int N>  struct product_type_selector<M,N,1>              { enum { ret = OuterProduct }; };
00107 template<int Depth>     struct product_type_selector<1,    1,    Depth>  { enum { ret = InnerProduct }; };
00108 template<>              struct product_type_selector<1,    1,    1>      { enum { ret = InnerProduct }; };
00109 template<>              struct product_type_selector<Small,1,    Small>  { enum { ret = CoeffBasedProductMode }; };
00110 template<>              struct product_type_selector<1,    Small,Small>  { enum { ret = CoeffBasedProductMode }; };
00111 template<>              struct product_type_selector<Small,Small,Small>  { enum { ret = CoeffBasedProductMode }; };
00112 template<>              struct product_type_selector<Small, Small, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00113 template<>              struct product_type_selector<Small, Large, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00114 template<>              struct product_type_selector<Large, Small, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00115 template<>              struct product_type_selector<1,    Large,Small>  { enum { ret = CoeffBasedProductMode }; };
00116 template<>              struct product_type_selector<1,    Large,Large>  { enum { ret = GemvProduct }; };
00117 template<>              struct product_type_selector<1,    Small,Large>  { enum { ret = CoeffBasedProductMode }; };
00118 template<>              struct product_type_selector<Large,1,    Small>  { enum { ret = CoeffBasedProductMode }; };
00119 template<>              struct product_type_selector<Large,1,    Large>  { enum { ret = GemvProduct }; };
00120 template<>              struct product_type_selector<Small,1,    Large>  { enum { ret = CoeffBasedProductMode }; };
00121 template<>              struct product_type_selector<Small,Small,Large>  { enum { ret = GemmProduct }; };
00122 template<>              struct product_type_selector<Large,Small,Large>  { enum { ret = GemmProduct }; };
00123 template<>              struct product_type_selector<Small,Large,Large>  { enum { ret = GemmProduct }; };
00124 template<>              struct product_type_selector<Large,Large,Large>  { enum { ret = GemmProduct }; };
00125 template<>              struct product_type_selector<Large,Small,Small>  { enum { ret = GemmProduct }; };
00126 template<>              struct product_type_selector<Small,Large,Small>  { enum { ret = GemmProduct }; };
00127 template<>              struct product_type_selector<Large,Large,Small>  { enum { ret = GemmProduct }; };
00128 
00129 } // end namespace internal
00130 
00131 /** \class ProductReturnType
00132   * \ingroup Core_Module
00133   *
00134   * \brief Helper class to get the correct and optimized returned type of operator*
00135   *
00136   * \param Lhs the type of the left-hand side
00137   * \param Rhs the type of the right-hand side
00138   * \param ProductMode the type of the product (determined automatically by internal::product_mode)
00139   *
00140   * This class defines the typename Type representing the optimized product expression
00141   * between two matrix expressions. In practice, using ProductReturnType<Lhs,Rhs>::Type
00142   * is the recommended way to define the result type of a function returning an expression
00143   * which involve a matrix product. The class Product should never be
00144   * used directly.
00145   *
00146   * \sa class Product, MatrixBase::operator*(const MatrixBase<OtherDerived>&)
00147   */
00148 template<typename Lhs, typename Rhs, int ProductType>
00149 struct ProductReturnType
00150 {
00151   // TODO use the nested type to reduce instanciations ????
00152 //   typedef typename internal::nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
00153 //   typedef typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
00154 
00155   typedef GeneralProduct<Lhs/*Nested*/, Rhs/*Nested*/, ProductType> Type;
00156 };
00157 
00158 template<typename Lhs, typename Rhs>
00159 struct ProductReturnType<Lhs,Rhs,CoeffBasedProductMode>
00160 {
00161   typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
00162   typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
00163   typedef CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit> Type;
00164 };
00165 
00166 template<typename Lhs, typename Rhs>
00167 struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
00168 {
00169   typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
00170   typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
00171   typedef CoeffBasedProduct<LhsNested, RhsNested, NestByRefBit> Type;
00172 };
00173 
00174 // this is a workaround for sun CC
00175 template<typename Lhs, typename Rhs>
00176 struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
00177 {};
00178 
00179 /***********************************************************************
00180 *  Implementation of Inner Vector Vector Product
00181 ***********************************************************************/
00182 
00183 // FIXME : maybe the "inner product" could return a Scalar
00184 // instead of a 1x1 matrix ??
00185 // Pro: more natural for the user
00186 // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
00187 // product ends up to a row-vector times col-vector product... To tackle this use
00188 // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
00189 
00190 namespace internal {
00191 
00192 template<typename Lhs, typename Rhs>
00193 struct traits<GeneralProduct<Lhs,Rhs,InnerProduct> >
00194  : traits<Matrix<typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> >
00195 {};
00196 
00197 }
00198 
00199 template<typename Lhs, typename Rhs>
00200 class GeneralProduct<Lhs, Rhs, InnerProduct>
00201   : internal::no_assignment_operator,
00202     public Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1>
00203 {
00204     typedef Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> Base;
00205   public:
00206     GeneralProduct(const Lhs& lhs, const Rhs& rhs)
00207     {
00208       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
00209         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00210 
00211       Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
00212     }
00213 
00214     /** Convertion to scalar */
00215     operator const typename Base::Scalar() const {
00216       return Base::coeff(0,0);
00217     }
00218 };
00219 
00220 /***********************************************************************
00221 *  Implementation of Outer Vector Vector Product
00222 ***********************************************************************/
00223 
00224 namespace internal {
00225 
00226 // Column major
00227 template<typename ProductType, typename Dest, typename Func>
00228 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const false_type&)
00229 {
00230   typedef typename Dest::Index Index;
00231   // FIXME make sure lhs is sequentially stored
00232   // FIXME not very good if rhs is real and lhs complex while alpha is real too
00233   const Index cols = dest.cols();
00234   for (Index j=0; j<cols; ++j)
00235     func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
00236 }
00237 
00238 // Row major
00239 template<typename ProductType, typename Dest, typename Func>
00240 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const true_type&) {
00241   typedef typename Dest::Index Index;
00242   // FIXME make sure rhs is sequentially stored
00243   // FIXME not very good if lhs is real and rhs complex while alpha is real too
00244   const Index rows = dest.rows();
00245   for (Index i=0; i<rows; ++i)
00246     func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
00247 }
00248 
00249 template<typename Lhs, typename Rhs>
00250 struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> >
00251  : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> >
00252 {};
00253 
00254 }
00255 
00256 template<typename Lhs, typename Rhs>
00257 class GeneralProduct<Lhs, Rhs, OuterProduct>
00258   : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs>
00259 {
00260     template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
00261     
00262   public:
00263     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
00264 
00265     GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00266     {
00267       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
00268         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00269     }
00270     
00271     struct set  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived()  = src; } };
00272     struct add  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
00273     struct sub  { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
00274     struct adds {
00275       Scalar m_scale;
00276       adds(const Scalar& s) : m_scale(s) {}
00277       template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const {
00278         dst.const_cast_derived() += m_scale * src;
00279       }
00280     };
00281     
00282     template<typename Dest>
00283     inline void evalTo(Dest& dest) const {
00284       internal::outer_product_selector_run(*this, dest, set(), is_row_major<Dest>());
00285     }
00286     
00287     template<typename Dest>
00288     inline void addTo(Dest& dest) const {
00289       internal::outer_product_selector_run(*this, dest, add(), is_row_major<Dest>());
00290     }
00291 
00292     template<typename Dest>
00293     inline void subTo(Dest& dest) const {
00294       internal::outer_product_selector_run(*this, dest, sub(), is_row_major<Dest>());
00295     }
00296 
00297     template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const
00298     {
00299       internal::outer_product_selector_run(*this, dest, adds(alpha), is_row_major<Dest>());
00300     }
00301 };
00302 
00303 /***********************************************************************
00304 *  Implementation of General Matrix Vector Product
00305 ***********************************************************************/
00306 
00307 /*  According to the shape/flags of the matrix we have to distinghish 3 different cases:
00308  *   1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
00309  *   2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
00310  *   3 - all other cases are handled using a simple loop along the outer-storage direction.
00311  *  Therefore we need a lower level meta selector.
00312  *  Furthermore, if the matrix is the rhs, then the product has to be transposed.
00313  */
00314 namespace internal {
00315 
00316 template<typename Lhs, typename Rhs>
00317 struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> >
00318  : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> >
00319 {};
00320 
00321 template<int Side, int StorageOrder, bool BlasCompatible>
00322 struct gemv_selector;
00323 
00324 } // end namespace internal
00325 
00326 template<typename Lhs, typename Rhs>
00327 class GeneralProduct<Lhs, Rhs, GemvProduct>
00328   : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs>
00329 {
00330   public:
00331     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
00332 
00333     typedef typename Lhs::Scalar LhsScalar;
00334     typedef typename Rhs::Scalar RhsScalar;
00335 
00336     GeneralProduct(const Lhs& a_lhs, const Rhs& a_rhs) : Base(a_lhs,a_rhs)
00337     {
00338 //       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value),
00339 //         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00340     }
00341 
00342     enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
00343     typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType;
00344 
00345     template<typename Dest> void scaleAndAddTo(Dest& dst, const Scalar& alpha) const
00346     {
00347       eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols());
00348       internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
00349                        bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha);
00350     }
00351 };
00352 
00353 namespace internal {
00354 
00355 // The vector is on the left => transposition
00356 template<int StorageOrder, bool BlasCompatible>
00357 struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible>
00358 {
00359   template<typename ProductType, typename Dest>
00360   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
00361   {
00362     Transpose<Dest> destT(dest);
00363     enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
00364     gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible>
00365       ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct>
00366         (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha);
00367   }
00368 };
00369 
00370 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if;
00371 
00372 template<typename Scalar,int Size,int MaxSize>
00373 struct gemv_static_vector_if<Scalar,Size,MaxSize,false>
00374 {
00375   EIGEN_STRONG_INLINE  Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; }
00376 };
00377 
00378 template<typename Scalar,int Size>
00379 struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
00380 {
00381   EIGEN_STRONG_INLINE Scalar* data() { return 0; }
00382 };
00383 
00384 template<typename Scalar,int Size,int MaxSize>
00385 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
00386 {
00387   #if EIGEN_ALIGN_STATICALLY
00388   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
00389   EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
00390   #else
00391   // Some architectures cannot align on the stack,
00392   // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
00393   enum {
00394     ForceAlignment  = internal::packet_traits<Scalar>::Vectorizable,
00395     PacketSize      = internal::packet_traits<Scalar>::size
00396   };
00397   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
00398   EIGEN_STRONG_INLINE Scalar* data() {
00399     return ForceAlignment
00400             ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
00401             : m_data.array;
00402   }
00403   #endif
00404 };
00405 
00406 template<> struct gemv_selector<OnTheRight,ColMajor,true>
00407 {
00408   template<typename ProductType, typename Dest>
00409   static inline void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
00410   {
00411     typedef typename ProductType::Index Index;
00412     typedef typename ProductType::LhsScalar   LhsScalar;
00413     typedef typename ProductType::RhsScalar   RhsScalar;
00414     typedef typename ProductType::Scalar      ResScalar;
00415     typedef typename ProductType::RealScalar  RealScalar;
00416     typedef typename ProductType::ActualLhsType ActualLhsType;
00417     typedef typename ProductType::ActualRhsType ActualRhsType;
00418     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00419     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00420     typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
00421 
00422     ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs());
00423     ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs());
00424 
00425     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00426                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00427 
00428     // make sure Dest is a compile-time vector type (bug 1166)
00429     typedef typename conditional<Dest::IsVectorAtCompileTime, Dest, typename Dest::ColXpr>::type ActualDest;
00430 
00431     enum {
00432       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
00433       // on, the other hand it is good for the cache to pack the vector anyways...
00434       EvalToDestAtCompileTime = (ActualDest::InnerStrideAtCompileTime==1),
00435       ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
00436       MightCannotUseDest = (ActualDest::InnerStrideAtCompileTime!=1) || ComplexByReal
00437     };
00438 
00439     gemv_static_vector_if<ResScalar,ActualDest::SizeAtCompileTime,ActualDest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
00440 
00441     bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
00442     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
00443     
00444     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
00445 
00446     ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
00447                                                   evalToDest ? dest.data() : static_dest.data());
00448     
00449     if(!evalToDest)
00450     {
00451       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00452       int size = dest.size();
00453       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00454       #endif
00455       if(!alphaIsCompatible)
00456       {
00457         MappedDest(actualDestPtr, dest.size()).setZero();
00458         compatibleAlpha = RhsScalar(1);
00459       }
00460       else
00461         MappedDest(actualDestPtr, dest.size()) = dest;
00462     }
00463 
00464     general_matrix_vector_product
00465       <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
00466         actualLhs.rows(), actualLhs.cols(),
00467         actualLhs.data(), actualLhs.outerStride(),
00468         actualRhs.data(), actualRhs.innerStride(),
00469         actualDestPtr, 1,
00470         compatibleAlpha);
00471 
00472     if (!evalToDest)
00473     {
00474       if(!alphaIsCompatible)
00475         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
00476       else
00477         dest = MappedDest(actualDestPtr, dest.size());
00478     }
00479   }
00480 };
00481 
00482 template<> struct gemv_selector<OnTheRight,RowMajor,true>
00483 {
00484   template<typename ProductType, typename Dest>
00485   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
00486   {
00487     typedef typename ProductType::LhsScalar LhsScalar;
00488     typedef typename ProductType::RhsScalar RhsScalar;
00489     typedef typename ProductType::Scalar    ResScalar;
00490     typedef typename ProductType::Index Index;
00491     typedef typename ProductType::ActualLhsType ActualLhsType;
00492     typedef typename ProductType::ActualRhsType ActualRhsType;
00493     typedef typename ProductType::_ActualRhsType _ActualRhsType;
00494     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00495     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00496 
00497     typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
00498     typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
00499 
00500     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00501                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00502 
00503     enum {
00504       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
00505       // on, the other hand it is good for the cache to pack the vector anyways...
00506       DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1
00507     };
00508 
00509     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
00510 
00511     ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
00512         DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
00513 
00514     if(!DirectlyUseRhs)
00515     {
00516       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00517       int size = actualRhs.size();
00518       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00519       #endif
00520       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
00521     }
00522 
00523     general_matrix_vector_product
00524       <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
00525         actualLhs.rows(), actualLhs.cols(),
00526         actualLhs.data(), actualLhs.outerStride(),
00527         actualRhsPtr, 1,
00528         dest.data(), dest.col(0).innerStride(), //NOTE  if dest is not a vector at compile-time, then dest.innerStride() might be wrong. (bug 1166)
00529         actualAlpha);
00530   }
00531 };
00532 
00533 template<> struct gemv_selector<OnTheRight,ColMajor,false>
00534 {
00535   template<typename ProductType, typename Dest>
00536   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
00537   {
00538     typedef typename Dest::Index Index;
00539     // TODO makes sure dest is sequentially stored in memory, otherwise use a temp
00540     const Index size = prod.rhs().rows();
00541     for(Index k=0; k<size; ++k)
00542       dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k);
00543   }
00544 };
00545 
00546 template<> struct gemv_selector<OnTheRight,RowMajor,false>
00547 {
00548   template<typename ProductType, typename Dest>
00549   static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha)
00550   {
00551     typedef typename Dest::Index Index;
00552     // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp
00553     const Index rows = prod.rows();
00554     for(Index i=0; i<rows; ++i)
00555       dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum();
00556   }
00557 };
00558 
00559 } // end namespace internal
00560 
00561 /***************************************************************************
00562 * Implementation of matrix base methods
00563 ***************************************************************************/
00564 
00565 /** \returns the matrix product of \c *this and \a other.
00566   *
00567   * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*().
00568   *
00569   * \sa lazyProduct(), operator*=(const MatrixBase&), Cwise::operator*()
00570   */
00571 template<typename Derived>
00572 template<typename OtherDerived>
00573 inline const typename ProductReturnType<Derived, OtherDerived>::Type
00574 MatrixBase<Derived>::operator* (const MatrixBase<OtherDerived> &other) const
00575 {
00576   // A note regarding the function declaration: In MSVC, this function will sometimes
00577   // not be inlined since DenseStorage is an unwindable object for dynamic
00578   // matrices and product types are holding a member to store the result.
00579   // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
00580   enum {
00581     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
00582                    || OtherDerived::RowsAtCompileTime==Dynamic
00583                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
00584     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
00585     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
00586   };
00587   // note to the lost user:
00588   //    * for a dot product use: v1.dot(v2)
00589   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
00590   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00591     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00592   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00593     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00594   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00595 #ifdef EIGEN_DEBUG_PRODUCT
00596   internal::product_type<Derived,OtherDerived>::debug();
00597 #endif
00598   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00599 }
00600 
00601 /** \returns an expression of the matrix product of \c *this and \a other without implicit evaluation.
00602   *
00603   * The returned product will behave like any other expressions: the coefficients of the product will be
00604   * computed once at a time as requested. This might be useful in some extremely rare cases when only
00605   * a small and no coherent fraction of the result's coefficients have to be computed.
00606   *
00607   * \warning This version of the matrix product can be much much slower. So use it only if you know
00608   * what you are doing and that you measured a true speed improvement.
00609   *
00610   * \sa operator*(const MatrixBase&)
00611   */
00612 template<typename Derived>
00613 template<typename OtherDerived>
00614 const typename LazyProductReturnType<Derived,OtherDerived>::Type
00615 MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const
00616 {
00617   enum {
00618     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
00619                    || OtherDerived::RowsAtCompileTime==Dynamic
00620                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
00621     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
00622     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
00623   };
00624   // note to the lost user:
00625   //    * for a dot product use: v1.dot(v2)
00626   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
00627   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00628     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00629   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00630     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00631   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00632 
00633   return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00634 }
00635 
00636 } // end namespace Eigen
00637 
00638 #endif // EIGEN_PRODUCT_H