From 88147e0a91733bde2c095cbe51c20b6d8a32ea88 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 5 Aug 2009 15:23:35 +0200 Subject: big refactoring in Product.h: - all specialized products now inherits ProductBase - the default product evaluated by Assign is still here, but it is currently enabled for small fixed sizes only - => this significantly speed up compilation for large matrices - I left the OuterProduct specialization empty as an exercise... --- Eigen/src/Core/Matrix.h | 15 - Eigen/src/Core/MatrixBase.h | 32 +- Eigen/src/Core/Product.h | 895 ++++++++++++++---------------- Eigen/src/Core/ProductBase.h | 36 +- Eigen/src/Core/util/ForwardDeclarations.h | 9 +- 5 files changed, 464 insertions(+), 523 deletions(-) (limited to 'Eigen') diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 848236bac..f58424ba2 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -339,13 +339,6 @@ class Matrix return Base::operator=(func); } - template - EIGEN_STRONG_INLINE Matrix& operator=(const ProductBase& other) - { - resize(other.rows(), other.cols()); - return Base::operator=(other); - } - using Base::operator +=; using Base::operator -=; using Base::operator *=; @@ -452,14 +445,6 @@ class Matrix other.evalTo(*this); } - template - EIGEN_STRONG_INLINE Matrix(const ProductBase& other) - { - _check_template_params(); - resize(other.rows(), other.cols()); - other.evalTo(*this); - } - /** Destructor */ inline ~Matrix() {} diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f94764731..6ec7ddbb7 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -318,17 +318,6 @@ template class MatrixBase Derived& operator-=(const AnyMatrixBase &other) { other.derived().subToDense(derived()); return derived(); } - - template - Derived& operator=(const ProductBase &other); - - template - Derived& operator+=(const ProductBase &other); - - template - Derived& operator-=(const ProductBase &other); - - template Derived& operator=(const ReturnByValue& func); @@ -337,14 +326,21 @@ template class MatrixBase template Derived& lazyAssign(const MatrixBase& other); - /** Overloaded for cache friendly product evaluation */ - template - Derived& lazyAssign(const Product& product); - /** Overloaded for cache friendly product evaluation */ template Derived& lazyAssign(const Flagged& other) { return lazyAssign(other._expression()); } + + template + Derived& lazyAssign(const ProductBase& other); + + template + Derived& operator+=(const Flagged, 0, + EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); + + template + Derived& operator-=(const Flagged, 0, + EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); #endif // not EIGEN_PARSED_BY_DOXYGEN CommaInitializer operator<< (const Scalar& s); @@ -412,12 +408,6 @@ template class MatrixBase template Derived& operator-=(const MatrixBase& other); - template - Derived& operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); - - template - Derived& operator-=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); - Derived& operator*=(const Scalar& other); Derived& operator/=(const Scalar& other); diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 1a32eb5de..151639993 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -26,15 +26,73 @@ #ifndef EIGEN_PRODUCT_H #define EIGEN_PRODUCT_H -/*************************** -*** Forward declarations *** -***************************/ +/** \class GeneralProduct + * + * \brief Expression of the product of two general matrices or vectors + * + * \param LhsNested the type used to store the left-hand side + * \param RhsNested the type used to store the right-hand side + * \param ProductMode the type of the product + * + * This class represents an expression of the product of two general matrices. + * We call a general matrix, a dense matrix with full storage. For instance, + * This excludes triangular, selfadjoint, and sparse matrices. + * It is the return type of the operator* between general matrices. Its template + * arguments are determined automatically by ProductReturnType. Therefore, + * GeneralProduct should never be used direclty. To determine the result type of a + * function which involves a matrix product, use ProductReturnType::Type. + * + * \sa ProductReturnType, MatrixBase::operator*(const MatrixBase&) + */ +template::value> +class GeneralProduct; -template -struct ei_product_coeff_impl; +template struct ei_product_type_selector; -template -struct ei_product_packet_impl; +enum { + Large = Dynamic, + Small = -Dynamic +}; + +enum { OuterProduct, InnerProduct, UnrolledProduct, GemvProduct, GemmProduct }; + +template struct ei_product_type +{ + enum { + Rows = Lhs::RowsAtCompileTime, + Cols = Rhs::ColsAtCompileTime, + Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime), + + value = ei_product_type_selector<(Rows>8 ? Large : Rows==1 ? 1 : Small), + (Cols>8 ? Large : Cols==1 ? 1 : Small), + (Depth>8 ? Large : Depth==1 ? 1 : Small)>::ret + }; +}; + +template struct ei_product_type_selector { enum { ret = OuterProduct }; }; +template struct ei_product_type_selector<1,1,Depth> { enum { ret = InnerProduct }; }; +template<> struct ei_product_type_selector<1,1,1> { enum { ret = InnerProduct }; }; +template<> struct ei_product_type_selector { enum { ret = UnrolledProduct }; }; +template<> struct ei_product_type_selector<1,Small,Small> { enum { ret = UnrolledProduct }; }; +template<> struct ei_product_type_selector { enum { ret = UnrolledProduct }; }; + +// template<> struct ei_product_type_selector { enum { ret = GemvProduct }; }; +// template<> struct ei_product_type_selector<1,Small,Small> { enum { ret = GemvProduct }; }; +// template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; + +template<> struct ei_product_type_selector<1,Large,Small> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1,Large,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1,Small,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; +template<> struct ei_product_type_selector { enum { ret = GemmProduct }; }; /** \class ProductReturnType * @@ -52,133 +110,365 @@ struct ei_product_packet_impl; * * \sa class Product, MatrixBase::operator*(const MatrixBase&) */ -template +template struct ProductReturnType +{ + // TODO use the nested type to reduce instanciations ???? +// typedef typename ei_nested::type LhsNested; +// typedef typename ei_nested::type RhsNested; + + typedef GeneralProduct Type; +}; + +template +struct ProductReturnType { typedef typename ei_nested::type LhsNested; typedef typename ei_nested::type RhsNested; + typedef GeneralProduct Type; +}; + - typedef Product Type; +/*********************************************************************** +* Implementation of General Matrix Matrix Product +***********************************************************************/ + +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +class GeneralProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) + + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + + template void addTo(Dest& dst, Scalar alpha) const + { + ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); + + const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs); + const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs) + * RhsBlasTraits::extractScalarFactor(m_rhs); + + ei_general_matrix_matrix_product< + Scalar, + (_ActualLhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(LhsBlasTraits::NeedToConjugate), + (_ActualRhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(RhsBlasTraits::NeedToConjugate), + (Dest::Flags&RowMajorBit)?RowMajor:ColMajor> + ::run( + this->rows(), this->cols(), lhs.cols(), + (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(), + (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(), + (Scalar*)&(dst.coeffRef(0,0)), dst.stride(), + actualAlpha); + } }; -// cache friendly specialization +/*********************************************************************** +* Implementation of Inner Vector Vector Product +***********************************************************************/ + template -struct ProductReturnType +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +class GeneralProduct + : public ProductBase, Lhs, Rhs> { - typedef typename ei_nested::type LhsNested; - typedef typename ei_nested::type - >::type RhsNested; + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) + + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} - typedef Product Type; + template void addTo(Dest& dst, Scalar alpha) const + { + ei_assert(dst.rows()==1 && dst.cols()==1); + dst.coeffRef(0,0) += (m_lhs.cwise()*m_rhs).sum(); + } }; -/* Helper class to determine the type of the product, can be either: - * - NormalProduct - * - CacheFriendlyProduct +/*********************************************************************** +* Implementation of Outer Vector Vector Product +***********************************************************************/ + +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +class GeneralProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) + + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + + template void addTo(Dest& dst, Scalar alpha) const + { + // TODO + } +}; + +/*********************************************************************** +* Implementation of General Matrix Vector Product +***********************************************************************/ + +/* According to the shape/flags of the matrix we have to distinghish 3 different cases: + * 1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine + * 2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine + * 3 - all other cases are handled using a simple loop along the outer-storage direction. + * Therefore we need a lower level meta selector. + * Furthermore, if the matrix is the rhs, then the product has to be transposed. */ -template struct ei_product_mode -{ - enum{ - // workaround sun studio: - LhsIsVectorAtCompileTime = ei_traits::ColsAtCompileTime==1 || ei_traits::ColsAtCompileTime==1, - value = ei_traits::MaxColsAtCompileTime == Dynamic - && ( ei_traits::MaxRowsAtCompileTime == Dynamic - || ei_traits::MaxColsAtCompileTime == Dynamic ) - && (!(Rhs::IsVectorAtCompileTime && (ei_traits::Flags&RowMajorBit) && (!(ei_traits::Flags&DirectAccessBit)))) - && (!(LhsIsVectorAtCompileTime && (!(ei_traits::Flags&RowMajorBit)) && (!(ei_traits::Flags&DirectAccessBit)))) - && (ei_is_same_type::Scalar, typename ei_traits::Scalar>::ret) - ? CacheFriendlyProduct - : NormalProduct }; -}; - -/** \class Product - * - * \brief Expression of the product of two matrices - * - * \param LhsNested the type used to store the left-hand side - * \param RhsNested the type used to store the right-hand side - * \param ProductMode the type of the product - * - * This class represents an expression of the product of two matrices. - * It is the return type of the operator* between matrices. Its template - * arguments are determined automatically by ProductReturnType. Therefore, - * Product should never be used direclty. To determine the result type of a - * function which involves a matrix product, use ProductReturnType::Type. - * - * \sa ProductReturnType, MatrixBase::operator*(const MatrixBase&) - */ -template -struct ei_traits > +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +struct ei_gemv_selector; + +template +class GeneralProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) + + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + + enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; + typedef typename ei_meta_if::ret MatrixType; + + template void addTo(Dest& dst, Scalar alpha) const + { + ei_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols()); + ei_gemv_selector::ActualAccess>::run(*this, dst, alpha); + } +}; + +// The vector is on the left => transposition +template +struct ei_gemv_selector +{ + template + static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) + { + Transpose destT(dest); + ei_gemv_selector + ::run(GeneralProduct,Transpose > + (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha); + } +}; + +template<> struct ei_gemv_selector +{ + template + static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) + { + typedef typename ProductType::Scalar Scalar; + typedef typename ProductType::ActualLhsType ActualLhsType; + typedef typename ProductType::ActualRhsType ActualRhsType; + typedef typename ProductType::LhsBlasTraits LhsBlasTraits; + typedef typename ProductType::RhsBlasTraits RhsBlasTraits; + + ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs()); + ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs()); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) + * RhsBlasTraits::extractScalarFactor(prod.rhs()); + + enum { + EvalToDest = (ei_packet_traits::size==1) + ||((Dest::Flags&ActualPacketAccessBit) && (!(Dest::Flags & RowMajorBit))) + }; + Scalar* EIGEN_RESTRICT actualDest; + if (EvalToDest) + actualDest = &dest.coeffRef(0); + else + { + actualDest = ei_aligned_stack_new(Scalar,dest.size()); + Map >(actualDest, dest.size()) = dest; + } + + ei_cache_friendly_product_colmajor_times_vector + ( + dest.size(), + &actualLhs.const_cast_derived().coeffRef(0,0), actualLhs.stride(), + actualRhs, actualDest, actualAlpha); + + if (!EvalToDest) + { + dest = Map >(actualDest, dest.size()); + ei_aligned_stack_delete(Scalar, actualDest, dest.size()); + } + } +}; + +template<> struct ei_gemv_selector +{ + template + static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) + { + typedef typename ProductType::Scalar Scalar; + typedef typename ProductType::ActualLhsType ActualLhsType; + typedef typename ProductType::ActualRhsType ActualRhsType; + typedef typename ProductType::_ActualRhsType _ActualRhsType; + typedef typename ProductType::LhsBlasTraits LhsBlasTraits; + typedef typename ProductType::RhsBlasTraits RhsBlasTraits; + + ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs()); + ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs()); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) + * RhsBlasTraits::extractScalarFactor(prod.rhs()); + + enum { + DirectlyUseRhs = ((ei_packet_traits::size==1) || (_ActualRhsType::Flags&ActualPacketAccessBit)) + && (!(_ActualRhsType::Flags & RowMajorBit)) + }; + + Scalar* EIGEN_RESTRICT rhs_data; + if (DirectlyUseRhs) + rhs_data = &actualRhs.const_cast_derived().coeffRef(0); + else + { + rhs_data = ei_aligned_stack_new(Scalar, actualRhs.size()); + Map >(rhs_data, actualRhs.size()) = actualRhs; + } + + ei_cache_friendly_product_rowmajor_times_vector + ( + &actualLhs.const_cast_derived().coeffRef(0,0), actualLhs.stride(), + rhs_data, prod.rhs().size(), dest, actualAlpha); + + if (!DirectlyUseRhs) ei_aligned_stack_delete(Scalar, rhs_data, prod.rhs().size()); + } +}; + +template<> struct ei_gemv_selector +{ + template + static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) + { + // TODO makes sure dest is sequentially stored in memory, otherwise use a temp + const int size = prod.rhs().rows(); + for(int k=0; k struct ei_gemv_selector +{ + template + static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) + { + // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp + const int rows = prod.rows(); + for(int i=0; i +struct ei_product_coeff_impl; + +template +struct ei_product_packet_impl; + +template +struct ei_traits > { - // clean the nested types: typedef typename ei_cleantype::type _LhsNested; typedef typename ei_cleantype::type _RhsNested; typedef typename ei_scalar_product_traits::ReturnType Scalar; - + enum { - LhsCoeffReadCost = _LhsNested::CoeffReadCost, - RhsCoeffReadCost = _RhsNested::CoeffReadCost, - LhsFlags = _LhsNested::Flags, - RhsFlags = _RhsNested::Flags, + LhsCoeffReadCost = _LhsNested::CoeffReadCost, + RhsCoeffReadCost = _RhsNested::CoeffReadCost, + LhsFlags = _LhsNested::Flags, + RhsFlags = _RhsNested::Flags, - RowsAtCompileTime = _LhsNested::RowsAtCompileTime, - ColsAtCompileTime = _RhsNested::ColsAtCompileTime, - InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), + RowsAtCompileTime = _LhsNested::RowsAtCompileTime, + ColsAtCompileTime = _RhsNested::ColsAtCompileTime, + InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), - MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, - MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, + MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, + MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, - LhsRowMajor = LhsFlags & RowMajorBit, - RhsRowMajor = RhsFlags & RowMajorBit, + LhsRowMajor = LhsFlags & RowMajorBit, + RhsRowMajor = RhsFlags & RowMajorBit, - CanVectorizeRhs = RhsRowMajor && (RhsFlags & PacketAccessBit) - && (ColsAtCompileTime == Dynamic || (ColsAtCompileTime % ei_packet_traits::size) == 0), + CanVectorizeRhs = RhsRowMajor && (RhsFlags & PacketAccessBit) + && (ColsAtCompileTime == Dynamic || (ColsAtCompileTime % ei_packet_traits::size) == 0), - CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) - && (RowsAtCompileTime == Dynamic || (RowsAtCompileTime % ei_packet_traits::size) == 0), + CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) + && (RowsAtCompileTime == Dynamic || (RowsAtCompileTime % ei_packet_traits::size) == 0), - EvalToRowMajor = RhsRowMajor && (ProductMode==(int)CacheFriendlyProduct ? LhsRowMajor : (!CanVectorizeLhs)), + EvalToRowMajor = RhsRowMajor && (!CanVectorizeLhs), - RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), + RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), - Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) - | EvalBeforeAssigningBit - | EvalBeforeNestingBit - | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0) - | (LhsFlags & RhsFlags & AlignedBit), + Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) + | EvalBeforeAssigningBit + | EvalBeforeNestingBit + | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0) + | (LhsFlags & RhsFlags & AlignedBit), - CoeffReadCost = InnerSize == Dynamic ? Dynamic - : InnerSize * (NumTraits::MulCost + LhsCoeffReadCost + RhsCoeffReadCost) - + (InnerSize - 1) * NumTraits::AddCost, + CoeffReadCost = InnerSize == Dynamic ? Dynamic + : InnerSize * (NumTraits::MulCost + LhsCoeffReadCost + RhsCoeffReadCost) + + (InnerSize - 1) * NumTraits::AddCost, - /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside - * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner - * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect - * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI. - */ - CanVectorizeInner = LhsRowMajor && (!RhsRowMajor) && (LhsFlags & RhsFlags & ActualPacketAccessBit) - && (InnerSize % ei_packet_traits::size == 0) - }; + /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside + * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner + * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect + * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI. + */ + CanVectorizeInner = LhsRowMajor && (!RhsRowMajor) && (LhsFlags & RhsFlags & ActualPacketAccessBit) + && (InnerSize % ei_packet_traits::size == 0) + }; }; -template class Product : ei_no_assignment_operator, - public MatrixBase > +template class GeneralProduct + : ei_no_assignment_operator, + public MatrixBase > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(Product) + EIGEN_GENERIC_PUBLIC_INTERFACE(GeneralProduct) private: - typedef typename ei_traits::_LhsNested _LhsNested; - typedef typename ei_traits::_RhsNested _RhsNested; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; enum { PacketSize = ei_packet_traits::size, - InnerSize = ei_traits::InnerSize, + InnerSize = ei_traits::InnerSize, Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT, - CanVectorizeInner = ei_traits::CanVectorizeInner + CanVectorizeInner = ei_traits::CanVectorizeInner }; typedef ei_product_coeff_impl class Product public: template - inline Product(const Lhs& lhs, const Rhs& rhs) + inline GeneralProduct(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) { // we don't allow taking products of matrices of different real types, as that wouldn't be vectorizable. @@ -200,23 +490,6 @@ template class Product && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); } - /** \internal - * compute \a res += \c *this using the cache friendly product. - */ - template - void _cacheFriendlyEvalAndAdd(DestDerived& res, Scalar alpha) const; - - /** \internal - * \returns whether it is worth it to use the cache friendly product. - */ - EIGEN_STRONG_INLINE bool _useCacheFriendlyProduct() const - { - // TODO do something more accurate here (especially for mat-vec products) - return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - && ( rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD - || cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD); - } - EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } EIGEN_STRONG_INLINE int cols() const { return m_rhs.cols(); } @@ -250,54 +523,11 @@ template class Product return res; } - EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } - EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } - protected: const LhsNested m_lhs; const RhsNested m_rhs; }; -/** \returns the matrix product of \c *this and \a other. - * - * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*(). - * - * \sa lazy(), operator*=(const MatrixBase&), Cwise::operator*() - */ -template -template -inline const typename ProductReturnType::Type -MatrixBase::operator*(const MatrixBase &other) const -{ - enum { - ProductIsValid = Derived::ColsAtCompileTime==Dynamic - || OtherDerived::RowsAtCompileTime==Dynamic - || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), - AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, - SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) - }; - // note to the lost user: - // * for a dot product use: v1.dot(v2) - // * for a coeff-wise product use: v1.cwise()*v2 - EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), - INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) - EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), - INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) - EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) - return typename ProductReturnType::Type(derived(), other.derived()); -} - -/** replaces \c *this by \c *this * \a other. - * - * \returns a reference to \c *this - */ -template -template -inline Derived & -MatrixBase::operator*=(const AnyMatrixBase &other) -{ - return derived() = derived() * other.derived(); -} /*************************************************************************** * Normal product .coeff() implementation (with meta-unrolling) @@ -509,335 +739,50 @@ struct ei_product_packet_impl::RowsAtCompileTime, - int LhsOrder = int(ei_traits::LhsFlags)&RowMajorBit ? RowMajor : ColMajor, - int LhsHasDirectAccess = ei_blas_traits::_LhsNested>::ActualAccess, - int RhsCols = ei_traits::ColsAtCompileTime, - int RhsOrder = int(ei_traits::RhsFlags)&RowMajorBit ? RowMajor : ColMajor, - int RhsHasDirectAccess = ei_blas_traits::_RhsNested>::ActualAccess> -struct ei_cache_friendly_product_selector -{ - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - product._cacheFriendlyEvalAndAdd(res, alpha); - } -}; - -// optimized colmajor * vector path -template -struct ei_cache_friendly_product_selector -{ - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ei_assert(alpha==typename ProductType::Scalar(1)); - const int size = product.rhs().rows(); - for (int k=0; k -struct ei_cache_friendly_product_selector -{ - typedef typename ProductType::Scalar Scalar; - typedef ei_blas_traits::_LhsNested> LhsProductTraits; - typedef ei_blas_traits::_RhsNested> RhsProductTraits; - - typedef typename LhsProductTraits::ExtractType ActualLhsType; - typedef typename RhsProductTraits::ExtractType ActualRhsType; - - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ActualLhsType actualLhs = LhsProductTraits::extract(product.lhs()); - ActualRhsType actualRhs = RhsProductTraits::extract(product.rhs()); - - Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(product.lhs()) - * RhsProductTraits::extractScalarFactor(product.rhs()); - - enum { - EvalToRes = (ei_packet_traits::size==1) - ||((DestDerived::Flags&ActualPacketAccessBit) && (!(DestDerived::Flags & RowMajorBit))) }; - Scalar* EIGEN_RESTRICT _res; - if (EvalToRes) - _res = &res.coeffRef(0); - else - { - _res = ei_aligned_stack_new(Scalar,res.size()); - Map >(_res, res.size()) = res; - } -// std::cerr << "colmajor * vector " << EvalToRes << "\n"; - ei_cache_friendly_product_colmajor_times_vector - ( - res.size(), - &actualLhs.const_cast_derived().coeffRef(0,0), actualLhs.stride(), - actualRhs, _res, actualAlpha); - - if (!EvalToRes) - { - res = Map >(_res, res.size()); - ei_aligned_stack_delete(Scalar, _res, res.size()); - } - } -}; - -// optimized vector * rowmajor path -template -struct ei_cache_friendly_product_selector -{ - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ei_assert(alpha==typename ProductType::Scalar(1)); - const int cols = product.lhs().cols(); - for (int j=0; j -struct ei_cache_friendly_product_selector -{ - typedef typename ProductType::Scalar Scalar; - typedef ei_blas_traits::_LhsNested> LhsProductTraits; - typedef ei_blas_traits::_RhsNested> RhsProductTraits; - - typedef typename LhsProductTraits::ExtractType ActualLhsType; - typedef typename RhsProductTraits::ExtractType ActualRhsType; - - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ActualLhsType actualLhs = LhsProductTraits::extract(product.lhs()); - ActualRhsType actualRhs = RhsProductTraits::extract(product.rhs()); - - Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(product.lhs()) - * RhsProductTraits::extractScalarFactor(product.rhs()); - - enum { - EvalToRes = (ei_packet_traits::size==1) - ||((DestDerived::Flags & ActualPacketAccessBit) && (DestDerived::Flags & RowMajorBit)) }; - Scalar* EIGEN_RESTRICT _res; - if (EvalToRes) - _res = &res.coeffRef(0); - else - { - _res = ei_aligned_stack_new(Scalar, res.size()); - Map >(_res, res.size()) = res; - } - - ei_cache_friendly_product_colmajor_times_vector - (res.size(), - &actualRhs.const_cast_derived().coeffRef(0,0), actualRhs.stride(), - actualLhs.transpose(), _res, actualAlpha); - - if (!EvalToRes) - { - res = Map >(_res, res.size()); - ei_aligned_stack_delete(Scalar, _res, res.size()); - } - } -}; - -// optimized rowmajor - vector product -template -struct ei_cache_friendly_product_selector -{ - typedef typename ProductType::Scalar Scalar; - - typedef ei_blas_traits::_LhsNested> LhsProductTraits; - typedef ei_blas_traits::_RhsNested> RhsProductTraits; - - typedef typename LhsProductTraits::ExtractType ActualLhsType; - typedef typename RhsProductTraits::ExtractType ActualRhsType; - typedef typename ei_cleantype::type _ActualRhsType; - - enum { - UseRhsDirectly = ((ei_packet_traits::size==1) || (_ActualRhsType::Flags&ActualPacketAccessBit)) - && (!(_ActualRhsType::Flags & RowMajorBit)) }; - - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ActualLhsType actualLhs = LhsProductTraits::extract(product.lhs()); - ActualRhsType actualRhs = RhsProductTraits::extract(product.rhs()); - - Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(product.lhs()) - * RhsProductTraits::extractScalarFactor(product.rhs()); - - Scalar* EIGEN_RESTRICT _rhs; - if (UseRhsDirectly) - _rhs = &actualRhs.const_cast_derived().coeffRef(0); - else - { - _rhs = ei_aligned_stack_new(Scalar, actualRhs.size()); - Map >(_rhs, actualRhs.size()) = actualRhs; - } - - ei_cache_friendly_product_rowmajor_times_vector - ( - &actualLhs.const_cast_derived().coeffRef(0,0), actualLhs.stride(), - _rhs, product.rhs().size(), res, actualAlpha); - - if (!UseRhsDirectly) ei_aligned_stack_delete(Scalar, _rhs, product.rhs().size()); - } -}; - -// optimized vector - colmajor product -template -struct ei_cache_friendly_product_selector +/** \returns the matrix product of \c *this and \a other. + * + * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*(). + * + * \sa lazy(), operator*=(const MatrixBase&), Cwise::operator*() + */ +template +template +inline const typename ProductReturnType::Type +MatrixBase::operator*(const MatrixBase &other) const { - typedef typename ProductType::Scalar Scalar; - - typedef ei_blas_traits::_LhsNested> LhsProductTraits; - typedef ei_blas_traits::_RhsNested> RhsProductTraits; - - typedef typename LhsProductTraits::ExtractType ActualLhsType; - typedef typename RhsProductTraits::ExtractType ActualRhsType; - typedef typename ei_cleantype::type _ActualLhsType; - enum { - UseLhsDirectly = ((ei_packet_traits::size==1) || (_ActualLhsType::Flags&ActualPacketAccessBit)) - && (_ActualLhsType::Flags & RowMajorBit) }; - - template - inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) - { - ActualLhsType actualLhs = LhsProductTraits::extract(product.lhs()); - ActualRhsType actualRhs = RhsProductTraits::extract(product.rhs()); - - Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(product.lhs()) - * RhsProductTraits::extractScalarFactor(product.rhs()); - - Scalar* EIGEN_RESTRICT _lhs; - if (UseLhsDirectly) - _lhs = &actualLhs.const_cast_derived().coeffRef(0); - else - { - _lhs = ei_aligned_stack_new(Scalar, actualLhs.size()); - Map >(_lhs, actualLhs.size()) = actualLhs; - } - - ei_cache_friendly_product_rowmajor_times_vector - ( - &actualRhs.const_cast_derived().coeffRef(0,0), actualRhs.stride(), - _lhs, product.lhs().size(), res, actualAlpha); - - if(!UseLhsDirectly) ei_aligned_stack_delete(Scalar, _lhs, product.lhs().size()); - } -}; - -// discard this case which has to be handled by the default path -// (we keep it to be sure to hit a compilation error if this is not the case) -template -struct ei_cache_friendly_product_selector -{}; - -// discard this case which has to be handled by the default path -// (we keep it to be sure to hit a compilation error if this is not the case) -template -struct ei_cache_friendly_product_selector -{}; - - -/** \internal - * Overloaded to perform an efficient C += A*B */ -template -template -inline Derived& -MatrixBase::operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) -{//std::cerr << "operator+=\n"; - if (other._expression()._useCacheFriendlyProduct()) - ei_cache_friendly_product_selector >::run(const_cast_derived(), other._expression(), Scalar(1)); - else { //std::cerr << "no cf\n"; - lazyAssign(derived() + other._expression()); - } - return derived(); + ProductIsValid = Derived::ColsAtCompileTime==Dynamic + || OtherDerived::RowsAtCompileTime==Dynamic + || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), + AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, + SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) + }; + // note to the lost user: + // * for a dot product use: v1.dot(v2) + // * for a coeff-wise product use: v1.cwise()*v2 + EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), + INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) + EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), + INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) + EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) + return typename ProductReturnType::Type(derived(), other.derived()); } -/** \internal - * Overloaded to perform an efficient C -= A*B */ -template -template -inline Derived& -MatrixBase::operator-=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) -{ - if (other._expression()._useCacheFriendlyProduct()) - ei_cache_friendly_product_selector >::run(const_cast_derived(), other._expression(), Scalar(-1)); - else - lazyAssign(derived() - other._expression()); - return derived(); -} -/** \internal - * Overloaded to perform an efficient C = A*B */ -template -template -inline Derived& MatrixBase::lazyAssign(const Product& product) -{ - if (product._useCacheFriendlyProduct()) - { - setZero(); - ei_cache_friendly_product_selector >::run(const_cast_derived(), product, Scalar(1)); - } - else - { - lazyAssign(Product(product.lhs(),product.rhs())); - } - return derived(); -} -template -template -inline void Product::_cacheFriendlyEvalAndAdd(DestDerived& res, Scalar alpha) const +/** replaces \c *this by \c *this * \a other. + * + * \returns a reference to \c *this + */ +template +template +inline Derived & +MatrixBase::operator*=(const AnyMatrixBase &other) { - typedef ei_blas_traits<_LhsNested> LhsProductTraits; - typedef ei_blas_traits<_RhsNested> RhsProductTraits; - - typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType; - typedef typename RhsProductTraits::DirectLinearAccessType ActualRhsType; - - typedef typename ei_cleantype::type _ActualLhsType; - typedef typename ei_cleantype::type _ActualRhsType; - - const ActualLhsType lhs = LhsProductTraits::extract(m_lhs); - const ActualRhsType rhs = RhsProductTraits::extract(m_rhs); - - Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(m_lhs) - * RhsProductTraits::extractScalarFactor(m_rhs); - - ei_general_matrix_matrix_product< - Scalar, - (_ActualLhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(LhsProductTraits::NeedToConjugate), - (_ActualRhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(RhsProductTraits::NeedToConjugate), - (DestDerived::Flags&RowMajorBit)?RowMajor:ColMajor> - ::run( - rows(), cols(), lhs.cols(), - (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(), - (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(), - (Scalar*)&(res.coeffRef(0,0)), res.stride(), - actualAlpha); + return derived() = derived() * other.derived(); } #endif // EIGEN_PRODUCT_H diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index 1f146babf..bae2302b3 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -39,11 +39,11 @@ struct ei_traits > ColsAtCompileTime = ei_traits::ColsAtCompileTime, MaxRowsAtCompileTime = ei_traits::MaxRowsAtCompileTime, MaxColsAtCompileTime = ei_traits::MaxColsAtCompileTime, - Flags = EvalBeforeNestingBit, + Flags = EvalBeforeNestingBit | EvalBeforeAssigningBit, CoeffReadCost = 0 // FIXME why is it needed ? }; }; -* + // enforce evaluation before nesting template struct ei_nested, N, EvalType> @@ -90,7 +90,11 @@ class ProductBase : public MatrixBase ProductBase(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) - {} + { + ei_assert(lhs.cols() == rhs.rows() + && "invalid matrix product" + && "if you wanted a coeff-wise or a dot product use the respective explicit functions"); + } inline int rows() const { return m_lhs.rows(); } inline int cols() const { return m_rhs.cols(); } @@ -115,6 +119,14 @@ class ProductBase : public MatrixBase return res; } + const Flagged lazy() const + { + return *this; + } + + const _LhsNested& lhs() const { return m_lhs; } + const _RhsNested& rhs() const { return m_rhs; } + protected: const LhsNested m_lhs; @@ -129,25 +141,33 @@ class ProductBase : public MatrixBase void coeffRef(int); }; +/** \internal + * Overloaded to perform an efficient C = (A*B).lazy() */ template template -Derived& MatrixBase::operator=(const ProductBase& other) +Derived& MatrixBase::lazyAssign(const ProductBase& other) { other.evalTo(derived()); return derived(); } +/** \internal + * Overloaded to perform an efficient C += (A*B).lazy() */ template template -Derived& MatrixBase::operator+=(const ProductBase& other) +Derived& MatrixBase::operator+=(const Flagged, 0, + EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { - other.addTo(derived()); return derived(); + other._expression().addTo(derived()); return derived(); } +/** \internal + * Overloaded to perform an efficient C -= (A*B).lazy() */ template template -Derived& MatrixBase::operator-=(const ProductBase& other) +Derived& MatrixBase::operator-=(const Flagged, 0, + EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) { - other.subTo(derived()); return derived(); + other._expression().subTo(derived()); return derived(); } #endif // EIGEN_PRODUCTBASE_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index d755445c1..c9ce174d2 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -48,8 +48,7 @@ template class CwiseNullaryOp; template class CwiseUnaryOp; template class CwiseUnaryView; template class CwiseBinaryOp; -template class ProductBase; -template class Product; +template class ProductBase; template class DiagonalBase; template class DiagonalWrapper; @@ -69,8 +68,10 @@ template class ReturnByValue; template class BandMatrix; -template struct ei_product_mode; -template::value> struct ProductReturnType; +template struct ei_product_type; +template::value> +struct ProductReturnType; template struct ei_scalar_sum_op; template struct ei_scalar_difference_op; -- cgit v1.2.3