Main MRPT website > C++ reference
MRPT logo

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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011