GeneralProduct.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_GENERAL_PRODUCT_H
00027 #define EIGEN_GENERAL_PRODUCT_H
00028 
00029 namespace Eigen { 
00030 
00050 template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value>
00051 class GeneralProduct;
00052 
00053 enum {
00054   Large = 2,
00055   Small = 3
00056 };
00057 
00058 namespace internal {
00059 
00060 template<int Rows, int Cols, int Depth> struct product_type_selector;
00061 
00062 template<int Size, int MaxSize> struct product_size_category
00063 {
00064   enum { is_large = MaxSize == Dynamic ||
00065                     Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD,
00066          value = is_large  ? Large
00067                : Size == 1 ? 1
00068                            : Small
00069   };
00070 };
00071 
00072 template<typename Lhs, typename Rhs> struct product_type
00073 {
00074   typedef typename remove_all<Lhs>::type _Lhs;
00075   typedef typename remove_all<Rhs>::type _Rhs;
00076   enum {
00077     MaxRows  = _Lhs::MaxRowsAtCompileTime,
00078     Rows  = _Lhs::RowsAtCompileTime,
00079     MaxCols  = _Rhs::MaxColsAtCompileTime,
00080     Cols  = _Rhs::ColsAtCompileTime,
00081     MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime,
00082                                            _Rhs::MaxRowsAtCompileTime),
00083     Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime,
00084                                         _Rhs::RowsAtCompileTime),
00085     LargeThreshold = EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
00086   };
00087 
00088   // the splitting into different lines of code here, introducing the _select enums and the typedef below,
00089   // is to work around an internal compiler error with gcc 4.1 and 4.2.
00090 private:
00091   enum {
00092     rows_select = product_size_category<Rows,MaxRows>::value,
00093     cols_select = product_size_category<Cols,MaxCols>::value,
00094     depth_select = product_size_category<Depth,MaxDepth>::value
00095   };
00096   typedef product_type_selector<rows_select, cols_select, depth_select> selector;
00097 
00098 public:
00099   enum {
00100     value = selector::ret
00101   };
00102 #ifdef EIGEN_DEBUG_PRODUCT
00103   static void debug()
00104   {
00105       EIGEN_DEBUG_VAR(Rows);
00106       EIGEN_DEBUG_VAR(Cols);
00107       EIGEN_DEBUG_VAR(Depth);
00108       EIGEN_DEBUG_VAR(rows_select);
00109       EIGEN_DEBUG_VAR(cols_select);
00110       EIGEN_DEBUG_VAR(depth_select);
00111       EIGEN_DEBUG_VAR(value);
00112   }
00113 #endif
00114 };
00115 
00116 
00117 /* The following allows to select the kind of product at compile time
00118  * based on the three dimensions of the product.
00119  * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
00120 // FIXME I'm not sure the current mapping is the ideal one.
00121 template<int M, int N>  struct product_type_selector<M,N,1>              { enum { ret = OuterProduct }; };
00122 template<int Depth>     struct product_type_selector<1,    1,    Depth>  { enum { ret = InnerProduct }; };
00123 template<>              struct product_type_selector<1,    1,    1>      { enum { ret = InnerProduct }; };
00124 template<>              struct product_type_selector<Small,1,    Small>  { enum { ret = CoeffBasedProductMode }; };
00125 template<>              struct product_type_selector<1,    Small,Small>  { enum { ret = CoeffBasedProductMode }; };
00126 template<>              struct product_type_selector<Small,Small,Small>  { enum { ret = CoeffBasedProductMode }; };
00127 template<>              struct product_type_selector<Small, Small, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00128 template<>              struct product_type_selector<Small, Large, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00129 template<>              struct product_type_selector<Large, Small, 1>    { enum { ret = LazyCoeffBasedProductMode }; };
00130 template<>              struct product_type_selector<1,    Large,Small>  { enum { ret = CoeffBasedProductMode }; };
00131 template<>              struct product_type_selector<1,    Large,Large>  { enum { ret = GemvProduct }; };
00132 template<>              struct product_type_selector<1,    Small,Large>  { enum { ret = CoeffBasedProductMode }; };
00133 template<>              struct product_type_selector<Large,1,    Small>  { enum { ret = CoeffBasedProductMode }; };
00134 template<>              struct product_type_selector<Large,1,    Large>  { enum { ret = GemvProduct }; };
00135 template<>              struct product_type_selector<Small,1,    Large>  { enum { ret = CoeffBasedProductMode }; };
00136 template<>              struct product_type_selector<Small,Small,Large>  { enum { ret = GemmProduct }; };
00137 template<>              struct product_type_selector<Large,Small,Large>  { enum { ret = GemmProduct }; };
00138 template<>              struct product_type_selector<Small,Large,Large>  { enum { ret = GemmProduct }; };
00139 template<>              struct product_type_selector<Large,Large,Large>  { enum { ret = GemmProduct }; };
00140 template<>              struct product_type_selector<Large,Small,Small>  { enum { ret = GemmProduct }; };
00141 template<>              struct product_type_selector<Small,Large,Small>  { enum { ret = GemmProduct }; };
00142 template<>              struct product_type_selector<Large,Large,Small>  { enum { ret = GemmProduct }; };
00143 
00144 } // end namespace internal
00145 
00163 template<typename Lhs, typename Rhs, int ProductType>
00164 struct ProductReturnType
00165 {
00166   // TODO use the nested type to reduce instanciations ????
00167 //   typedef typename internal::nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
00168 //   typedef typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
00169 
00170   typedef GeneralProduct<Lhs/*Nested*/, Rhs/*Nested*/, ProductType> Type;
00171 };
00172 
00173 template<typename Lhs, typename Rhs>
00174 struct ProductReturnType<Lhs,Rhs,CoeffBasedProductMode>
00175 {
00176   typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
00177   typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
00178   typedef CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit> Type;
00179 };
00180 
00181 template<typename Lhs, typename Rhs>
00182 struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
00183 {
00184   typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested;
00185   typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested;
00186   typedef CoeffBasedProduct<LhsNested, RhsNested, NestByRefBit> Type;
00187 };
00188 
00189 // this is a workaround for sun CC
00190 template<typename Lhs, typename Rhs>
00191 struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode>
00192 {};
00193 
00194 /***********************************************************************
00195 *  Implementation of Inner Vector Vector Product
00196 ***********************************************************************/
00197 
00198 // FIXME : maybe the "inner product" could return a Scalar
00199 // instead of a 1x1 matrix ??
00200 // Pro: more natural for the user
00201 // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix
00202 // product ends up to a row-vector times col-vector product... To tackle this use
00203 // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x);
00204 
00205 namespace internal {
00206 
00207 template<typename Lhs, typename Rhs>
00208 struct traits<GeneralProduct<Lhs,Rhs,InnerProduct> >
00209  : traits<Matrix<typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> >
00210 {};
00211 
00212 }
00213 
00214 template<typename Lhs, typename Rhs>
00215 class GeneralProduct<Lhs, Rhs, InnerProduct>
00216   : internal::no_assignment_operator,
00217     public Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1>
00218 {
00219     typedef Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> Base;
00220   public:
00221     GeneralProduct(const Lhs& lhs, const Rhs& rhs)
00222     {
00223       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
00224         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00225 
00226       Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
00227     }
00228 
00230     operator const typename Base::Scalar() const {
00231       return Base::coeff(0,0);
00232     }
00233 };
00234 
00235 /***********************************************************************
00236 *  Implementation of Outer Vector Vector Product
00237 ***********************************************************************/
00238 
00239 namespace internal {
00240 template<int StorageOrder> struct outer_product_selector;
00241 
00242 template<typename Lhs, typename Rhs>
00243 struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> >
00244  : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> >
00245 {};
00246 
00247 }
00248 
00249 template<typename Lhs, typename Rhs>
00250 class GeneralProduct<Lhs, Rhs, OuterProduct>
00251   : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs>
00252 {
00253   public:
00254     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
00255 
00256     GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00257     {
00258       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
00259         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00260     }
00261 
00262     template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
00263     {
00264       internal::outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha);
00265     }
00266 };
00267 
00268 namespace internal {
00269 
00270 template<> struct outer_product_selector<ColMajor> {
00271   template<typename ProductType, typename Dest>
00272   static EIGEN_DONT_INLINE void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) {
00273     typedef typename Dest::Index Index;
00274     // FIXME make sure lhs is sequentially stored
00275     // FIXME not very good if rhs is real and lhs complex while alpha is real too
00276     const Index cols = dest.cols();
00277     for (Index j=0; j<cols; ++j)
00278       dest.col(j) += (alpha * prod.rhs().coeff(j)) * prod.lhs();
00279   }
00280 };
00281 
00282 template<> struct outer_product_selector<RowMajor> {
00283   template<typename ProductType, typename Dest>
00284   static EIGEN_DONT_INLINE void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) {
00285     typedef typename Dest::Index Index;
00286     // FIXME make sure rhs is sequentially stored
00287     // FIXME not very good if lhs is real and rhs complex while alpha is real too
00288     const Index rows = dest.rows();
00289     for (Index i=0; i<rows; ++i)
00290       dest.row(i) += (alpha * prod.lhs().coeff(i)) * prod.rhs();
00291   }
00292 };
00293 
00294 } // end namespace internal
00295 
00296 /***********************************************************************
00297 *  Implementation of General Matrix Vector Product
00298 ***********************************************************************/
00299 
00300 /*  According to the shape/flags of the matrix we have to distinghish 3 different cases:
00301  *   1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
00302  *   2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
00303  *   3 - all other cases are handled using a simple loop along the outer-storage direction.
00304  *  Therefore we need a lower level meta selector.
00305  *  Furthermore, if the matrix is the rhs, then the product has to be transposed.
00306  */
00307 namespace internal {
00308 
00309 template<typename Lhs, typename Rhs>
00310 struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> >
00311  : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> >
00312 {};
00313 
00314 template<int Side, int StorageOrder, bool BlasCompatible>
00315 struct gemv_selector;
00316 
00317 } // end namespace internal
00318 
00319 template<typename Lhs, typename Rhs>
00320 class GeneralProduct<Lhs, Rhs, GemvProduct>
00321   : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs>
00322 {
00323   public:
00324     EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
00325 
00326     typedef typename Lhs::Scalar LhsScalar;
00327     typedef typename Rhs::Scalar RhsScalar;
00328 
00329     GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00330     {
00331 //       EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value),
00332 //         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00333     }
00334 
00335     enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
00336     typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType;
00337 
00338     template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00339     {
00340       eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols());
00341       internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
00342                        bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha);
00343     }
00344 };
00345 
00346 namespace internal {
00347 
00348 // The vector is on the left => transposition
00349 template<int StorageOrder, bool BlasCompatible>
00350 struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible>
00351 {
00352   template<typename ProductType, typename Dest>
00353   static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha)
00354   {
00355     Transpose<Dest> destT(dest);
00356     enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
00357     gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible>
00358       ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct>
00359         (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha);
00360   }
00361 };
00362 
00363 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if;
00364 
00365 template<typename Scalar,int Size,int MaxSize>
00366 struct gemv_static_vector_if<Scalar,Size,MaxSize,false>
00367 {
00368   EIGEN_STRONG_INLINE  Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; }
00369 };
00370 
00371 template<typename Scalar,int Size>
00372 struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
00373 {
00374   EIGEN_STRONG_INLINE Scalar* data() { return 0; }
00375 };
00376 
00377 template<typename Scalar,int Size,int MaxSize>
00378 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
00379 {
00380   #if EIGEN_ALIGN_STATICALLY
00381   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
00382   EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
00383   #else
00384   // Some architectures cannot align on the stack,
00385   // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
00386   enum {
00387     ForceAlignment  = internal::packet_traits<Scalar>::Vectorizable,
00388     PacketSize      = internal::packet_traits<Scalar>::size
00389   };
00390   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
00391   EIGEN_STRONG_INLINE Scalar* data() {
00392     return ForceAlignment
00393             ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
00394             : m_data.array;
00395   }
00396   #endif
00397 };
00398 
00399 template<> struct gemv_selector<OnTheRight,ColMajor,true>
00400 {
00401   template<typename ProductType, typename Dest>
00402   static inline void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha)
00403   {
00404     typedef typename ProductType::Index Index;
00405     typedef typename ProductType::LhsScalar   LhsScalar;
00406     typedef typename ProductType::RhsScalar   RhsScalar;
00407     typedef typename ProductType::Scalar      ResScalar;
00408     typedef typename ProductType::RealScalar  RealScalar;
00409     typedef typename ProductType::ActualLhsType ActualLhsType;
00410     typedef typename ProductType::ActualRhsType ActualRhsType;
00411     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00412     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00413     typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest;
00414 
00415     ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs());
00416     ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs());
00417 
00418     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00419                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00420 
00421     enum {
00422       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
00423       // on, the other hand it is good for the cache to pack the vector anyways...
00424       EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1,
00425       ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
00426       MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal
00427     };
00428 
00429     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
00430 
00431     bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0));
00432     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
00433     
00434     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
00435 
00436     ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
00437                                                   evalToDest ? dest.data() : static_dest.data());
00438     
00439     if(!evalToDest)
00440     {
00441       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00442       int size = dest.size();
00443       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00444       #endif
00445       if(!alphaIsCompatible)
00446       {
00447         MappedDest(actualDestPtr, dest.size()).setZero();
00448         compatibleAlpha = RhsScalar(1);
00449       }
00450       else
00451         MappedDest(actualDestPtr, dest.size()) = dest;
00452     }
00453 
00454     general_matrix_vector_product
00455       <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
00456         actualLhs.rows(), actualLhs.cols(),
00457         actualLhs.data(), actualLhs.outerStride(),
00458         actualRhs.data(), actualRhs.innerStride(),
00459         actualDestPtr, 1,
00460         compatibleAlpha);
00461 
00462     if (!evalToDest)
00463     {
00464       if(!alphaIsCompatible)
00465         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
00466       else
00467         dest = MappedDest(actualDestPtr, dest.size());
00468     }
00469   }
00470 };
00471 
00472 template<> struct gemv_selector<OnTheRight,RowMajor,true>
00473 {
00474   template<typename ProductType, typename Dest>
00475   static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha)
00476   {
00477     typedef typename ProductType::LhsScalar LhsScalar;
00478     typedef typename ProductType::RhsScalar RhsScalar;
00479     typedef typename ProductType::Scalar    ResScalar;
00480     typedef typename ProductType::Index Index;
00481     typedef typename ProductType::ActualLhsType ActualLhsType;
00482     typedef typename ProductType::ActualRhsType ActualRhsType;
00483     typedef typename ProductType::_ActualRhsType _ActualRhsType;
00484     typedef typename ProductType::LhsBlasTraits LhsBlasTraits;
00485     typedef typename ProductType::RhsBlasTraits RhsBlasTraits;
00486 
00487     typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
00488     typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
00489 
00490     ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs())
00491                                   * RhsBlasTraits::extractScalarFactor(prod.rhs());
00492 
00493     enum {
00494       // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
00495       // on, the other hand it is good for the cache to pack the vector anyways...
00496       DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1
00497     };
00498 
00499     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
00500 
00501     ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
00502         DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
00503 
00504     if(!DirectlyUseRhs)
00505     {
00506       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00507       int size = actualRhs.size();
00508       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
00509       #endif
00510       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
00511     }
00512 
00513     general_matrix_vector_product
00514       <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run(
00515         actualLhs.rows(), actualLhs.cols(),
00516         actualLhs.data(), actualLhs.outerStride(),
00517         actualRhsPtr, 1,
00518         dest.data(), dest.innerStride(),
00519         actualAlpha);
00520   }
00521 };
00522 
00523 template<> struct gemv_selector<OnTheRight,ColMajor,false>
00524 {
00525   template<typename ProductType, typename Dest>
00526   static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha)
00527   {
00528     typedef typename Dest::Index Index;
00529     // TODO makes sure dest is sequentially stored in memory, otherwise use a temp
00530     const Index size = prod.rhs().rows();
00531     for(Index k=0; k<size; ++k)
00532       dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k);
00533   }
00534 };
00535 
00536 template<> struct gemv_selector<OnTheRight,RowMajor,false>
00537 {
00538   template<typename ProductType, typename Dest>
00539   static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha)
00540   {
00541     typedef typename Dest::Index Index;
00542     // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp
00543     const Index rows = prod.rows();
00544     for(Index i=0; i<rows; ++i)
00545       dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum();
00546   }
00547 };
00548 
00549 } // end namespace internal
00550 
00551 /***************************************************************************
00552 * Implementation of matrix base methods
00553 ***************************************************************************/
00554 
00561 template<typename Derived>
00562 template<typename OtherDerived>
00563 inline const typename ProductReturnType<Derived, OtherDerived>::Type
00564 MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
00565 {
00566   // A note regarding the function declaration: In MSVC, this function will sometimes
00567   // not be inlined since DenseStorage is an unwindable object for dynamic
00568   // matrices and product types are holding a member to store the result.
00569   // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
00570   enum {
00571     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
00572                    || OtherDerived::RowsAtCompileTime==Dynamic
00573                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
00574     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
00575     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
00576   };
00577   // note to the lost user:
00578   //    * for a dot product use: v1.dot(v2)
00579   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
00580   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00581     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00582   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00583     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00584   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00585 #ifdef EIGEN_DEBUG_PRODUCT
00586   internal::product_type<Derived,OtherDerived>::debug();
00587 #endif
00588   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00589 }
00590 
00602 template<typename Derived>
00603 template<typename OtherDerived>
00604 const typename LazyProductReturnType<Derived,OtherDerived>::Type
00605 MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const
00606 {
00607   enum {
00608     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
00609                    || OtherDerived::RowsAtCompileTime==Dynamic
00610                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
00611     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
00612     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
00613   };
00614   // note to the lost user:
00615   //    * for a dot product use: v1.dot(v2)
00616   //    * for a coeff-wise product use: v1.cwiseProduct(v2)
00617   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00618     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00619   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00620     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00621   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00622 
00623   return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00624 }
00625 
00626 } // end namespace Eigen
00627 
00628 #endif // EIGEN_PRODUCT_H