vcglib: Product.h Source File

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.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 /***************************
00030 *** Forward declarations ***
00031 ***************************/
00032 
00033 template<int VectorizationMode, int Index, typename Lhs, typename Rhs, typename RetScalar>
00034 struct ei_product_coeff_impl;
00035 
00036 template<int StorageOrder, int Index, typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00037 struct ei_product_packet_impl;
00038 
00055 template<typename Lhs, typename Rhs, int ProductMode>
00056 struct ProductReturnType
00057 {
00058   typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
00059   typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
00060 
00061   typedef Product<LhsNested, RhsNested, ProductMode> Type;
00062 };
00063 
00064 // cache friendly specialization
00065 // note that there is a DiagonalProduct specialization in DiagonalProduct.h
00066 template<typename Lhs, typename Rhs>
00067 struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
00068 {
00069   typedef const Lhs& LhsNested;
00070   typedef const Rhs& RhsNested;
00071 
00072   typedef Product<LhsNested, RhsNested, CacheFriendlyProduct> Type;
00073 };
00074 
00075 /*  Helper class to determine the type of the product, can be either:
00076  *    - NormalProduct
00077  *    - CacheFriendlyProduct
00078  *    - DiagonalProduct
00079  */
00080 template<typename Lhs, typename Rhs> struct ei_product_mode
00081 {
00082   enum{
00083 
00084     value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal)
00085           ? DiagonalProduct
00086           : Lhs::MaxColsAtCompileTime == Dynamic
00087             && ( Lhs::MaxRowsAtCompileTime == Dynamic
00088               || Rhs::MaxColsAtCompileTime == Dynamic )
00089             && (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit)  && (!(Lhs::Flags&DirectAccessBit))))
00090             && (!(Lhs::IsVectorAtCompileTime && (!(Rhs::Flags&RowMajorBit)) && (!(Rhs::Flags&DirectAccessBit))))
00091             && (ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)
00092           ? CacheFriendlyProduct
00093           : NormalProduct };
00094 };
00095 
00112 template<typename LhsNested, typename RhsNested, int ProductMode>
00113 struct ei_traits<Product<LhsNested, RhsNested, ProductMode> >
00114 {
00115   // clean the nested types:
00116   typedef typename ei_cleantype<LhsNested>::type _LhsNested;
00117   typedef typename ei_cleantype<RhsNested>::type _RhsNested;
00118   typedef typename ei_scalar_product_traits<typename _LhsNested::Scalar, typename _RhsNested::Scalar>::ReturnType Scalar;
00119 
00120   enum {
00121     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
00122     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
00123     LhsFlags = _LhsNested::Flags,
00124     RhsFlags = _RhsNested::Flags,
00125 
00126     RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
00127     ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
00128     InnerSize = EIGEN_SIZE_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
00129 
00130     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
00131     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
00132 
00133     LhsRowMajor = LhsFlags & RowMajorBit,
00134     RhsRowMajor = RhsFlags & RowMajorBit,
00135 
00136     CanVectorizeRhs = RhsRowMajor && (RhsFlags & PacketAccessBit)
00137                     && (ColsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
00138 
00139     CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit)
00140                     && (RowsAtCompileTime % ei_packet_traits<Scalar>::size == 0),
00141 
00142     EvalToRowMajor = RhsRowMajor && (ProductMode==(int)CacheFriendlyProduct ? LhsRowMajor : (!CanVectorizeLhs)),
00143 
00144     RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit)|DirectAccessBit),
00145 
00146     Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
00147           | EvalBeforeAssigningBit
00148           | EvalBeforeNestingBit
00149           | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0)
00150           | (LhsFlags & RhsFlags & AlignedBit),
00151 
00152     CoeffReadCost = InnerSize == Dynamic ? Dynamic
00153                   : InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
00154                     + (InnerSize - 1) * NumTraits<Scalar>::AddCost,
00155 
00156     /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside
00157      * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner
00158      * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
00159      * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
00160      */
00161     CanVectorizeInner = LhsRowMajor && (!RhsRowMajor) && (LhsFlags & RhsFlags & ActualPacketAccessBit)
00162                       && (InnerSize % ei_packet_traits<Scalar>::size == 0)
00163   };
00164 };
00165 
00166 template<typename LhsNested, typename RhsNested, int ProductMode> class Product : ei_no_assignment_operator,
00167   public MatrixBase<Product<LhsNested, RhsNested, ProductMode> >
00168 {
00169   public:
00170 
00171     EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
00172 
00173   private:
00174 
00175     typedef typename ei_traits<Product>::_LhsNested _LhsNested;
00176     typedef typename ei_traits<Product>::_RhsNested _RhsNested;
00177 
00178     enum {
00179       PacketSize = ei_packet_traits<Scalar>::size,
00180       InnerSize  = ei_traits<Product>::InnerSize,
00181       Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
00182       CanVectorizeInner = ei_traits<Product>::CanVectorizeInner
00183     };
00184 
00185     typedef ei_product_coeff_impl<CanVectorizeInner ? InnerVectorization : NoVectorization,
00186                                   Unroll ? InnerSize-1 : Dynamic,
00187                                   _LhsNested, _RhsNested, Scalar> ScalarCoeffImpl;
00188 
00189   public:
00190     typedef typename ei_nested<_LhsNested,_RhsNested::ColsAtCompileTime>::type LhsNestedX;
00191     typedef typename ei_nested<_RhsNested,_LhsNested::RowsAtCompileTime>::type RhsNestedX;
00192     typedef typename ei_cleantype<LhsNestedX>::type _LhsNestedX;
00193     typedef typename ei_cleantype<RhsNestedX>::type _RhsNestedX;
00194   
00195     
00196     enum {
00197       LhsNestedFlags = _LhsNestedX::Flags,
00198       RhsNestedFlags = _RhsNestedX::Flags
00199     };
00200     
00201 
00202     template<typename Lhs, typename Rhs>
00203     inline Product(const Lhs& lhs, const Rhs& rhs)
00204       : m_lhs(lhs), m_rhs(rhs)
00205     {
00206       // we don't allow taking products of matrices of different real types, as that wouldn't be vectorizable.
00207       // We still allow to mix T and complex<T>.
00208       EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret),
00209         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
00210       ei_assert(lhs.cols() == rhs.rows()
00211         && "invalid matrix product"
00212         && "if you wanted a coeff-wise or a dot product use the respective explicit functions");
00213     }
00214 
00218     template<typename DestDerived>
00219     void _cacheFriendlyEvalAndAdd(DestDerived& res) const;
00220 
00224     EIGEN_STRONG_INLINE bool _useCacheFriendlyProduct() const
00225     {
00226       return  m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
00227               && (  rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
00228                  || cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD);
00229     }
00230 
00231     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
00232     EIGEN_STRONG_INLINE int cols() const { return m_rhs.cols(); }
00233 
00234     EIGEN_STRONG_INLINE const Scalar coeff(int row, int col) const
00235     {
00236       Scalar res;
00237       ScalarCoeffImpl::run(row, col, m_lhs, m_rhs, res);
00238       return res;
00239     }
00240 
00241     /* Allow index-based non-packet access. It is impossible though to allow index-based packed access,
00242      * which is why we don't set the LinearAccessBit.
00243      */
00244     EIGEN_STRONG_INLINE const Scalar coeff(int index) const
00245     {
00246       Scalar res;
00247       const int row = RowsAtCompileTime == 1 ? 0 : index;
00248       const int col = RowsAtCompileTime == 1 ? index : 0;
00249       ScalarCoeffImpl::run(row, col, m_lhs, m_rhs, res);
00250       return res;
00251     }
00252 
00253     template<int LoadMode>
00254     EIGEN_STRONG_INLINE const PacketScalar packet(int row, int col) const
00255     {
00256       PacketScalar res;
00257       ei_product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor,
00258                                    Unroll ? InnerSize-1 : Dynamic,
00259                                    _LhsNested, _RhsNested, PacketScalar, LoadMode>
00260         ::run(row, col, m_lhs, m_rhs, res);
00261       return res;
00262     }
00263 
00264     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
00265     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
00266 
00267   protected:
00268     const LhsNested m_lhs;
00269     const RhsNested m_rhs;
00270 };
00271 
00278 template<typename Derived>
00279 template<typename OtherDerived>
00280 inline const typename ProductReturnType<Derived,OtherDerived>::Type
00281 MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
00282 {
00283   enum {
00284     ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
00285                    || OtherDerived::RowsAtCompileTime==Dynamic
00286                    || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
00287     AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
00288     SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
00289   };
00290   // note to the lost user:
00291   //    * for a dot product use: v1.dot(v2)
00292   //    * for a coeff-wise product use: v1.cwise()*v2
00293   EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00294     INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00295   EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00296     INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00297   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00298   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00299 }
00300 
00305 template<typename Derived>
00306 template<typename OtherDerived>
00307 inline Derived &
00308 MatrixBase<Derived>::operator*=(const MatrixBase<OtherDerived> &other)
00309 {
00310   return derived() = derived() * other.derived();
00311 }
00312 
00313 /***************************************************************************
00314 * Normal product .coeff() implementation (with meta-unrolling)
00315 ***************************************************************************/
00316 
00317 /**************************************
00318 *** Scalar path  - no vectorization ***
00319 **************************************/
00320 
00321 template<int Index, typename Lhs, typename Rhs, typename RetScalar>
00322 struct ei_product_coeff_impl<NoVectorization, Index, Lhs, Rhs, RetScalar>
00323 {
00324   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar &res)
00325   {
00326     ei_product_coeff_impl<NoVectorization, Index-1, Lhs, Rhs, RetScalar>::run(row, col, lhs, rhs, res);
00327     res += lhs.coeff(row, Index) * rhs.coeff(Index, col);
00328   }
00329 };
00330 
00331 template<typename Lhs, typename Rhs, typename RetScalar>
00332 struct ei_product_coeff_impl<NoVectorization, 0, Lhs, Rhs, RetScalar>
00333 {
00334   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar &res)
00335   {
00336     res = lhs.coeff(row, 0) * rhs.coeff(0, col);
00337   }
00338 };
00339 
00340 template<typename Lhs, typename Rhs, typename RetScalar>
00341 struct ei_product_coeff_impl<NoVectorization, Dynamic, Lhs, Rhs, RetScalar>
00342 {
00343   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar& res)
00344   {
00345     ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
00346     res = lhs.coeff(row, 0) * rhs.coeff(0, col);
00347       for(int i = 1; i < lhs.cols(); ++i)
00348         res += lhs.coeff(row, i) * rhs.coeff(i, col);
00349   }
00350 };
00351 
00352 // prevent buggy user code from causing an infinite recursion
00353 template<typename Lhs, typename Rhs, typename RetScalar>
00354 struct ei_product_coeff_impl<NoVectorization, -1, Lhs, Rhs, RetScalar>
00355 {
00356   EIGEN_STRONG_INLINE static void run(int, int, const Lhs&, const Rhs&, RetScalar&) {}
00357 };
00358 
00359 /*******************************************
00360 *** Scalar path with inner vectorization ***
00361 *******************************************/
00362 
00363 template<int Index, typename Lhs, typename Rhs, typename PacketScalar>
00364 struct ei_product_coeff_vectorized_unroller
00365 {
00366   enum { PacketSize = ei_packet_traits<typename Lhs::Scalar>::size };
00367   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::PacketScalar &pres)
00368   {
00369     ei_product_coeff_vectorized_unroller<Index-PacketSize, Lhs, Rhs, PacketScalar>::run(row, col, lhs, rhs, pres);
00370     pres = ei_padd(pres, ei_pmul( lhs.template packet<Aligned>(row, Index) , rhs.template packet<Aligned>(Index, col) ));
00371   }
00372 };
00373 
00374 template<typename Lhs, typename Rhs, typename PacketScalar>
00375 struct ei_product_coeff_vectorized_unroller<0, Lhs, Rhs, PacketScalar>
00376 {
00377   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::PacketScalar &pres)
00378   {
00379     pres = ei_pmul(lhs.template packet<Aligned>(row, 0) , rhs.template packet<Aligned>(0, col));
00380   }
00381 };
00382 
00383 template<int Index, typename Lhs, typename Rhs, typename RetScalar>
00384 struct ei_product_coeff_impl<InnerVectorization, Index, Lhs, Rhs, RetScalar>
00385 {
00386   typedef typename Lhs::PacketScalar PacketScalar;
00387   enum { PacketSize = ei_packet_traits<typename Lhs::Scalar>::size };
00388   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar &res)
00389   {
00390     PacketScalar pres;
00391     ei_product_coeff_vectorized_unroller<Index+1-PacketSize, Lhs, Rhs, PacketScalar>::run(row, col, lhs, rhs, pres);
00392     ei_product_coeff_impl<NoVectorization,Index,Lhs,Rhs,RetScalar>::run(row, col, lhs, rhs, res);
00393     res = ei_predux(pres);
00394   }
00395 };
00396 
00397 template<typename Lhs, typename Rhs, int LhsRows = Lhs::RowsAtCompileTime, int RhsCols = Rhs::ColsAtCompileTime>
00398 struct ei_product_coeff_vectorized_dyn_selector
00399 {
00400   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res)
00401   {
00402     res = ei_dot_impl<
00403       Block<Lhs, 1, ei_traits<Lhs>::ColsAtCompileTime>,
00404       Block<Rhs, ei_traits<Rhs>::RowsAtCompileTime, 1>,
00405       LinearVectorization, NoUnrolling>::run(lhs.row(row), rhs.col(col));
00406   }
00407 };
00408 
00409 // NOTE the 3 following specializations are because taking .col(0) on a vector is a bit slower
00410 // NOTE maybe they are now useless since we have a specialization for Block<Matrix>
00411 template<typename Lhs, typename Rhs, int RhsCols>
00412 struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,RhsCols>
00413 {
00414   EIGEN_STRONG_INLINE static void run(int /*row*/, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res)
00415   {
00416     res = ei_dot_impl<
00417       Lhs,
00418       Block<Rhs, ei_traits<Rhs>::RowsAtCompileTime, 1>,
00419       LinearVectorization, NoUnrolling>::run(lhs, rhs.col(col));
00420   }
00421 };
00422 
00423 template<typename Lhs, typename Rhs, int LhsRows>
00424 struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,LhsRows,1>
00425 {
00426   EIGEN_STRONG_INLINE static void run(int row, int /*col*/, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res)
00427   {
00428     res = ei_dot_impl<
00429       Block<Lhs, 1, ei_traits<Lhs>::ColsAtCompileTime>,
00430       Rhs,
00431       LinearVectorization, NoUnrolling>::run(lhs.row(row), rhs);
00432   }
00433 };
00434 
00435 template<typename Lhs, typename Rhs>
00436 struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,1>
00437 {
00438   EIGEN_STRONG_INLINE static void run(int /*row*/, int /*col*/, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res)
00439   {
00440     res = ei_dot_impl<
00441       Lhs,
00442       Rhs,
00443       LinearVectorization, NoUnrolling>::run(lhs, rhs);
00444   }
00445 };
00446 
00447 template<typename Lhs, typename Rhs, typename RetScalar>
00448 struct ei_product_coeff_impl<InnerVectorization, Dynamic, Lhs, Rhs, RetScalar>
00449 {
00450   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res)
00451   {
00452     ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs>::run(row, col, lhs, rhs, res);
00453   }
00454 };
00455 
00456 /*******************
00457 *** Packet path  ***
00458 *******************/
00459 
00460 template<int Index, typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00461 struct ei_product_packet_impl<RowMajor, Index, Lhs, Rhs, PacketScalar, LoadMode>
00462 {
00463   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar &res)
00464   {
00465     ei_product_packet_impl<RowMajor, Index-1, Lhs, Rhs, PacketScalar, LoadMode>::run(row, col, lhs, rhs, res);
00466     res =  ei_pmadd(ei_pset1(lhs.coeff(row, Index)), rhs.template packet<LoadMode>(Index, col), res);
00467   }
00468 };
00469 
00470 template<int Index, typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00471 struct ei_product_packet_impl<ColMajor, Index, Lhs, Rhs, PacketScalar, LoadMode>
00472 {
00473   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar &res)
00474   {
00475     ei_product_packet_impl<ColMajor, Index-1, Lhs, Rhs, PacketScalar, LoadMode>::run(row, col, lhs, rhs, res);
00476     res =  ei_pmadd(lhs.template packet<LoadMode>(row, Index), ei_pset1(rhs.coeff(Index, col)), res);
00477   }
00478 };
00479 
00480 template<typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00481 struct ei_product_packet_impl<RowMajor, 0, Lhs, Rhs, PacketScalar, LoadMode>
00482 {
00483   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar &res)
00484   {
00485     res = ei_pmul(ei_pset1(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col));
00486   }
00487 };
00488 
00489 template<typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00490 struct ei_product_packet_impl<ColMajor, 0, Lhs, Rhs, PacketScalar, LoadMode>
00491 {
00492   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar &res)
00493   {
00494     res = ei_pmul(lhs.template packet<LoadMode>(row, 0), ei_pset1(rhs.coeff(0, col)));
00495   }
00496 };
00497 
00498 template<typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00499 struct ei_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMode>
00500 {
00501   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar& res)
00502   {
00503     ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
00504     res = ei_pmul(ei_pset1(lhs.coeff(row, 0)),rhs.template packet<LoadMode>(0, col));
00505       for(int i = 1; i < lhs.cols(); ++i)
00506         res =  ei_pmadd(ei_pset1(lhs.coeff(row, i)), rhs.template packet<LoadMode>(i, col), res);
00507   }
00508 };
00509 
00510 template<typename Lhs, typename Rhs, typename PacketScalar, int LoadMode>
00511 struct ei_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMode>
00512 {
00513   EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, PacketScalar& res)
00514   {
00515     ei_assert(lhs.cols()>0 && "you are using a non initialized matrix");
00516     res = ei_pmul(lhs.template packet<LoadMode>(row, 0), ei_pset1(rhs.coeff(0, col)));
00517       for(int i = 1; i < lhs.cols(); ++i)
00518         res =  ei_pmadd(lhs.template packet<LoadMode>(row, i), ei_pset1(rhs.coeff(i, col)), res);
00519   }
00520 };
00521 
00522 /***************************************************************************
00523 * Cache friendly product callers and specific nested evaluation strategies
00524 ***************************************************************************/
00525 
00526 template<typename Scalar, typename RhsType>
00527 static void ei_cache_friendly_product_colmajor_times_vector(
00528   int size, const Scalar* lhs, int lhsStride, const RhsType& rhs, Scalar* res);
00529 
00530 template<typename Scalar, typename ResType>
00531 static void ei_cache_friendly_product_rowmajor_times_vector(
00532   const Scalar* lhs, int lhsStride, const Scalar* rhs, int rhsSize, ResType& res);
00533 
00534 template<typename ProductType,
00535   int LhsRows  = ei_traits<ProductType>::RowsAtCompileTime,
00536   int LhsOrder = int(ProductType::LhsNestedFlags)&RowMajorBit ? RowMajor : ColMajor,
00537   int LhsHasDirectAccess = int(ProductType::LhsNestedFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess,
00538   int RhsCols  = ProductType::ColsAtCompileTime,
00539   int RhsOrder = int(ProductType::RhsNestedFlags)&RowMajorBit ? RowMajor : ColMajor,
00540   int RhsHasDirectAccess = int(ProductType::RhsNestedFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess>
00541 struct ei_cache_friendly_product_selector
00542 {
00543   template<typename DestDerived>
00544   inline static void run(DestDerived& res, const ProductType& product)
00545   {
00546     product._cacheFriendlyEvalAndAdd(res);
00547   }
00548 };
00549 
00550 // optimized colmajor * vector path
00551 template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess>
00552 struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,NoDirectAccess,1,RhsOrder,RhsAccess>
00553 {
00554   template<typename DestDerived>
00555   inline static void run(DestDerived& res, const ProductType& product)
00556   {
00557     typename ProductType::LhsNestedX lhs(product.lhs());
00558     typename ProductType::RhsNestedX rhs(product.rhs());
00559     
00560     const int size = rhs.rows();
00561     for (int k=0; k<size; ++k)
00562         res += rhs.coeff(k) * lhs.col(k);
00563   }
00564 };
00565 
00566 // optimized cache friendly colmajor * vector path for matrix with direct access flag
00567 // NOTE this path could also be enabled for expressions if we add runtime align queries
00568 template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess>
00569 struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirectAccess,1,RhsOrder,RhsAccess>
00570 {
00571   typedef typename ProductType::Scalar Scalar;
00572 
00573   template<typename DestDerived>
00574   inline static void run(DestDerived& res, const ProductType& product)
00575   {
00576     typename ProductType::LhsNestedX lhs(product.lhs());
00577     typename ProductType::RhsNestedX rhs(product.rhs());
00578     
00579     enum {
00580       EvalToRes = (ei_packet_traits<Scalar>::size==1)
00581                 ||((DestDerived::Flags&ActualPacketAccessBit) && (!(DestDerived::Flags & RowMajorBit))) };
00582     Scalar* EIGEN_RESTRICT _res;
00583     if (EvalToRes)
00584        _res = &res.coeffRef(0);
00585     else
00586     {
00587       _res = ei_aligned_stack_new(Scalar,res.size());
00588       Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1,ColMajor> >(_res, res.size()) = res;
00589     }
00590     ei_cache_friendly_product_colmajor_times_vector(res.size(),
00591       &lhs.const_cast_derived().coeffRef(0,0), lhs.stride(),
00592       rhs, _res);
00593 
00594     if (!EvalToRes)
00595     {
00596       res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size());
00597       ei_aligned_stack_delete(Scalar, _res, res.size());
00598     }
00599   }
00600 };
00601 
00602 // optimized vector * rowmajor path
00603 template<typename ProductType, int LhsOrder, int LhsAccess, int RhsCols>
00604 struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCols,RowMajor,NoDirectAccess>
00605 {
00606   template<typename DestDerived>
00607   inline static void run(DestDerived& res, const ProductType& product)
00608   {
00609     typename ProductType::LhsNestedX lhs(product.lhs());
00610     typename ProductType::RhsNestedX rhs(product.rhs());
00611     
00612     const int cols = lhs.cols();
00613     for (int j=0; j<cols; ++j)
00614       res += lhs.coeff(j) * rhs.row(j);
00615   }
00616 };
00617 
00618 // optimized cache friendly vector * rowmajor path for matrix with direct access flag
00619 // NOTE this path coul also be enabled for expressions if we add runtime align queries
00620 template<typename ProductType, int LhsOrder, int LhsAccess, int RhsCols>
00621 struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCols,RowMajor,HasDirectAccess>
00622 {
00623   typedef typename ProductType::Scalar Scalar;
00624 
00625   template<typename DestDerived>
00626   inline static void run(DestDerived& res, const ProductType& product)
00627   {
00628     typename ProductType::LhsNestedX lhs(product.lhs());
00629     typename ProductType::RhsNestedX rhs(product.rhs());
00630     
00631     enum {
00632       EvalToRes = (ei_packet_traits<Scalar>::size==1)
00633                 ||((DestDerived::Flags & ActualPacketAccessBit) && (DestDerived::Flags & RowMajorBit)) };
00634     Scalar* EIGEN_RESTRICT _res;
00635     if (EvalToRes)
00636        _res = &res.coeffRef(0);
00637     else
00638     {
00639       _res = ei_aligned_stack_new(Scalar, res.size());
00640       Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size()) = res;
00641     }
00642     ei_cache_friendly_product_colmajor_times_vector(res.size(),
00643       &rhs.const_cast_derived().coeffRef(0,0), rhs.stride(),
00644       lhs.transpose(), _res);
00645 
00646     if (!EvalToRes)
00647     {
00648       res = Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size());
00649       ei_aligned_stack_delete(Scalar, _res, res.size());
00650     }
00651   }
00652 };
00653 
00654 // optimized rowmajor - vector product
00655 template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess>
00656 struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirectAccess,1,RhsOrder,RhsAccess>
00657 {
00658   typedef typename ProductType::Scalar Scalar;
00659   typedef typename ei_traits<ProductType>::_RhsNested Rhs;
00660   enum {
00661       UseRhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (ProductType::RhsNestedFlags&ActualPacketAccessBit))
00662                      && (ProductType::RhsNestedFlags&DirectAccessBit)
00663                      && (!(ProductType::RhsNestedFlags & RowMajorBit)) };
00664 
00665   template<typename DestDerived>
00666   inline static void run(DestDerived& res, const ProductType& product)
00667   {
00668     typename ProductType::LhsNestedX lhs(product.lhs());
00669     typename ProductType::RhsNestedX rhs(product.rhs());
00670   
00671     Scalar* EIGEN_RESTRICT _rhs;
00672     if (UseRhsDirectly)
00673        _rhs = &rhs.const_cast_derived().coeffRef(0);
00674     else
00675     {
00676       _rhs = ei_aligned_stack_new(Scalar, rhs.size());
00677       Map<Matrix<Scalar,Rhs::SizeAtCompileTime,1,ColMajor> >(_rhs, rhs.size()) = rhs;
00678     }
00679     ei_cache_friendly_product_rowmajor_times_vector(&lhs.const_cast_derived().coeffRef(0,0), lhs.stride(),
00680                                                     _rhs, rhs.size(), res);
00681 
00682     if (!UseRhsDirectly) ei_aligned_stack_delete(Scalar, _rhs, rhs.size());
00683   }
00684 };
00685 
00686 // optimized vector - colmajor product
00687 template<typename ProductType, int LhsOrder, int LhsAccess, int RhsCols>
00688 struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCols,ColMajor,HasDirectAccess>
00689 {
00690   typedef typename ProductType::Scalar Scalar;
00691   typedef typename ei_traits<ProductType>::_LhsNested Lhs;
00692   enum {
00693       UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (ProductType::LhsNestedFlags&ActualPacketAccessBit))
00694                      && (ProductType::LhsNestedFlags&DirectAccessBit)
00695                      && (ProductType::LhsNestedFlags & RowMajorBit) };
00696 
00697   template<typename DestDerived>
00698   inline static void run(DestDerived& res, const ProductType& product)
00699   {
00700     typename ProductType::LhsNestedX lhs(product.lhs());
00701     typename ProductType::RhsNestedX rhs(product.rhs());
00702   
00703     Scalar* EIGEN_RESTRICT _lhs;
00704     if (UseLhsDirectly)
00705        _lhs = &lhs.const_cast_derived().coeffRef(0);
00706     else
00707     {
00708       _lhs = ei_aligned_stack_new(Scalar, lhs.size());
00709       Map<Matrix<Scalar,1,Lhs::SizeAtCompileTime,ColMajor> >(_lhs, lhs.size()) = lhs;
00710     }
00711     ei_cache_friendly_product_rowmajor_times_vector(&rhs.const_cast_derived().coeffRef(0,0), rhs.stride(),
00712                                                     _lhs, lhs.size(), res);
00713 
00714     if(!UseLhsDirectly) ei_aligned_stack_delete(Scalar, _lhs, lhs.size());
00715   }
00716 };
00717 
00718 // discard this case which has to be handled by the default path
00719 // (we keep it to be sure to hit a compilation error if this is not the case)
00720 template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess>
00721 struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,NoDirectAccess,1,RhsOrder,RhsAccess>
00722 {};
00723 
00724 // discard this case which has to be handled by the default path
00725 // (we keep it to be sure to hit a compilation error if this is not the case)
00726 template<typename ProductType, int LhsOrder, int LhsAccess, int RhsCols>
00727 struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCols,ColMajor,NoDirectAccess>
00728 {};
00729 
00730 
00732 template<typename Derived>
00733 template<typename Lhs,typename Rhs>
00734 inline Derived&
00735 MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other)
00736 {
00737   if (other._expression()._useCacheFriendlyProduct())
00738     ei_cache_friendly_product_selector<Product<Lhs,Rhs,CacheFriendlyProduct> >::run(const_cast_derived(), other._expression());
00739   else
00740   {
00741     typedef typename ei_cleantype<Lhs>::type _Lhs;
00742     typedef typename ei_cleantype<Rhs>::type _Rhs;
00743   
00744     typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
00745     typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
00746 
00747     Product<LhsNested,RhsNested,NormalProduct> prod(other._expression().lhs(),other._expression().rhs());
00748     
00749     lazyAssign(derived() + prod);
00750   }
00751   return derived();
00752 }
00753 
00754 template<typename Derived>
00755 template<typename Lhs, typename Rhs>
00756 inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFriendlyProduct>& product)
00757 {
00758   if (product._useCacheFriendlyProduct())
00759   {
00760     setZero();
00761     ei_cache_friendly_product_selector<Product<Lhs,Rhs,CacheFriendlyProduct> >::run(const_cast_derived(), product);
00762   }
00763   else
00764   {
00765     typedef typename ei_cleantype<Lhs>::type _Lhs;
00766     typedef typename ei_cleantype<Rhs>::type _Rhs;
00767     
00768     typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
00769     typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
00770 
00771     typedef Product<LhsNested,RhsNested,NormalProduct> NormalProduct;
00772     NormalProduct normal_prod(product.lhs(),product.rhs());
00773 
00774     lazyAssign<NormalProduct>(normal_prod);
00775   }
00776   return derived();
00777 }
00778 
00779 template<typename T,int StorageOrder> struct ei_product_copy_rhs
00780 {
00781   typedef typename ei_meta_if<
00782          (ei_traits<T>::Flags & RowMajorBit)
00783       || (!(ei_traits<T>::Flags & DirectAccessBit)),
00784       typename ei_plain_matrix_type_column_major<T>::type,
00785       const T&
00786     >::ret type;
00787 };
00788 
00789 template<typename T> struct ei_product_copy_rhs<T,RowMajorBit>
00790 {
00791   typedef typename ei_meta_if<
00792       (!(ei_traits<T>::Flags & DirectAccessBit)),
00793       typename ei_plain_matrix_type<T>::type,
00794       const T&
00795     >::ret type;
00796 };
00797 
00798 template<typename T,int StorageOrder> struct ei_product_copy_lhs
00799 {
00800   typedef typename ei_meta_if<
00801       (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
00802       typename ei_plain_matrix_type_row_major<T>::type,
00803       const T&
00804     >::ret type;
00805 };
00806 
00807 template<typename T> struct ei_product_copy_lhs<T,RowMajorBit>
00808 {
00809   typedef typename ei_meta_if<
00810          ((ei_traits<T>::Flags & RowMajorBit)==0)
00811       || (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
00812       typename ei_plain_matrix_type_row_major<T>::type,
00813       const T&
00814     >::ret type;
00815 };
00816 
00817 template<typename Lhs, typename Rhs, int ProductMode>
00818 template<typename DestDerived>
00819 inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const
00820 {
00821   typedef typename ei_product_copy_lhs<_LhsNested,DestDerived::Flags&RowMajorBit>::type LhsCopy;
00822   typedef typename ei_unref<LhsCopy>::type _LhsCopy;
00823   typedef typename ei_product_copy_rhs<_RhsNested,DestDerived::Flags&RowMajorBit>::type RhsCopy;
00824   typedef typename ei_unref<RhsCopy>::type _RhsCopy;
00825   LhsCopy lhs(m_lhs);
00826   RhsCopy rhs(m_rhs);
00827   ei_cache_friendly_product<Scalar>(
00828     rows(), cols(), lhs.cols(),
00829     _LhsCopy::Flags&RowMajorBit, (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(),
00830     _RhsCopy::Flags&RowMajorBit, (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(),
00831     DestDerived::Flags&RowMajorBit, (Scalar*)&(res.coeffRef(0,0)), res.stride()
00832   );
00833 
00834 }
00835 
00836 #endif // EIGEN_PRODUCT_H