diff options
Diffstat (limited to 'Eigen/src/Core/GeneralProduct.h')
| -rw-r--r-- | Eigen/src/Core/GeneralProduct.h | 286 |
1 files changed, 265 insertions, 21 deletions
diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h index 229d12c3f..734815986 100644 --- a/Eigen/src/Core/GeneralProduct.h +++ b/Eigen/src/Core/GeneralProduct.h @@ -13,6 +13,7 @@ namespace Eigen { +#ifndef EIGEN_TEST_EVALUATORS /** \class GeneralProduct * \ingroup Core_Module * @@ -34,6 +35,8 @@ namespace Eigen { */ template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value> class GeneralProduct; +#endif // EIGEN_TEST_EVALUATORS + enum { Large = 2, @@ -81,7 +84,8 @@ private: public: enum { - value = selector::ret + value = selector::ret, + ret = selector::ret }; #ifdef EIGEN_DEBUG_PRODUCT static void debug() @@ -97,6 +101,31 @@ public: #endif }; +// template<typename Lhs, typename Rhs> struct product_tag +// { +// private: +// +// typedef typename remove_all<Lhs>::type _Lhs; +// typedef typename remove_all<Rhs>::type _Rhs; +// enum { +// Rows = _Lhs::RowsAtCompileTime, +// Cols = _Rhs::ColsAtCompileTime, +// Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime, _Rhs::RowsAtCompileTime) +// }; +// +// enum { +// rows_select = Rows==1 ? int(Rows) : int(Large), +// cols_select = Cols==1 ? int(Cols) : int(Large), +// depth_select = Depth==1 ? int(Depth) : int(Large) +// }; +// typedef product_type_selector<rows_select, cols_select, depth_select> selector; +// +// public: +// enum { +// ret = selector::ret +// }; +// +// }; /* The following allows to select the kind of product at compile time * based on the three dimensions of the product. @@ -127,6 +156,7 @@ template<> struct product_type_selector<Large,Large,Small> { enum } // end namespace internal +#ifndef EIGEN_TEST_EVALUATORS /** \class ProductReturnType * \ingroup Core_Module * @@ -174,6 +204,7 @@ struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> template<typename Lhs, typename Rhs> struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> {}; +#endif /*********************************************************************** * Implementation of Inner Vector Vector Product @@ -185,6 +216,7 @@ struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedPr // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix // product ends up to a row-vector times col-vector product... To tackle this use // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x); +#ifndef EIGEN_TEST_EVALUATORS namespace internal { @@ -215,11 +247,12 @@ class GeneralProduct<Lhs, Rhs, InnerProduct> return Base::coeff(0,0); } }; - +#endif // EIGEN_TEST_EVALUATORS /*********************************************************************** * Implementation of Outer Vector Vector Product ***********************************************************************/ +#ifndef EIGEN_TEST_EVALUATORS namespace internal { // Column major @@ -299,6 +332,8 @@ class GeneralProduct<Lhs, Rhs, OuterProduct> } }; +#endif // EIGEN_TEST_EVALUATORS + /*********************************************************************** * Implementation of General Matrix Vector Product ***********************************************************************/ @@ -312,16 +347,22 @@ class GeneralProduct<Lhs, Rhs, OuterProduct> */ namespace internal { +#ifndef EIGEN_TEST_EVALUATORS template<typename Lhs, typename Rhs> struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> > : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> > {}; - template<int Side, int StorageOrder, bool BlasCompatible> struct gemv_selector; +#endif +#ifdef EIGEN_ENABLE_EVALUATORS +template<int Side, int StorageOrder, bool BlasCompatible> +struct gemv_dense_sense_selector; +#endif } // end namespace internal +#ifndef EIGEN_TEST_EVALUATORS template<typename Lhs, typename Rhs> class GeneralProduct<Lhs, Rhs, GemvProduct> : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> @@ -348,24 +389,10 @@ class GeneralProduct<Lhs, Rhs, GemvProduct> bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha); } }; +#endif namespace internal { -// The vector is on the left => transposition -template<int StorageOrder, bool BlasCompatible> -struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible> -{ - template<typename ProductType, typename Dest> - static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) - { - Transpose<Dest> destT(dest); - enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor }; - gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible> - ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct> - (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha); - } -}; - template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if; template<typename Scalar,int Size,int MaxSize> @@ -402,6 +429,23 @@ struct gemv_static_vector_if<Scalar,Size,MaxSize,true> #endif }; +#ifndef EIGEN_TEST_EVALUATORS + +// The vector is on the left => transposition +template<int StorageOrder, bool BlasCompatible> +struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible> +{ + template<typename ProductType, typename Dest> + static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) + { + Transpose<Dest> destT(dest); + enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor }; + gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible> + ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct> + (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha); + } +}; + template<> struct gemv_selector<OnTheRight,ColMajor,true> { template<typename ProductType, typename Dest> @@ -552,6 +596,179 @@ template<> struct gemv_selector<OnTheRight,RowMajor,false> } }; +#endif // EIGEN_TEST_EVALUATORS + +#ifdef EIGEN_ENABLE_EVALUATORS + +// The vector is on the left => transposition +template<int StorageOrder, bool BlasCompatible> +struct gemv_dense_sense_selector<OnTheLeft,StorageOrder,BlasCompatible> +{ + template<typename Lhs, typename Rhs, typename Dest> + static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha) + { + Transpose<Dest> destT(dest); + enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor }; + gemv_dense_sense_selector<OnTheRight,OtherStorageOrder,BlasCompatible> + ::run(rhs.transpose(), lhs.transpose(), destT, alpha); + } +}; + +template<> struct gemv_dense_sense_selector<OnTheRight,ColMajor,true> +{ + template<typename Lhs, typename Rhs, typename Dest> + static inline void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha) + { + typedef typename Dest::Index Index; + typedef typename Lhs::Scalar LhsScalar; + typedef typename Rhs::Scalar RhsScalar; + typedef typename Dest::Scalar ResScalar; + typedef typename Dest::RealScalar RealScalar; + + typedef internal::blas_traits<Lhs> LhsBlasTraits; + typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType; + typedef internal::blas_traits<Rhs> RhsBlasTraits; + typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; + + typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest; + + ActualLhsType actualLhs = LhsBlasTraits::extract(lhs); + ActualRhsType actualRhs = RhsBlasTraits::extract(rhs); + + ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs) + * RhsBlasTraits::extractScalarFactor(rhs); + + enum { + // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 + // on, the other hand it is good for the cache to pack the vector anyways... + EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1, + ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex), + MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal + }; + + gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest; + + bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); + bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; + + RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha); + + ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(), + evalToDest ? dest.data() : static_dest.data()); + + if(!evalToDest) + { + #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN + int size = dest.size(); + EIGEN_DENSE_STORAGE_CTOR_PLUGIN + #endif + if(!alphaIsCompatible) + { + MappedDest(actualDestPtr, dest.size()).setZero(); + compatibleAlpha = RhsScalar(1); + } + else + MappedDest(actualDestPtr, dest.size()) = dest; + } + + general_matrix_vector_product + <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( + actualLhs.rows(), actualLhs.cols(), + actualLhs.data(), actualLhs.outerStride(), + actualRhs.data(), actualRhs.innerStride(), + actualDestPtr, 1, + compatibleAlpha); + + if (!evalToDest) + { + if(!alphaIsCompatible) + dest += actualAlpha * MappedDest(actualDestPtr, dest.size()); + else + dest = MappedDest(actualDestPtr, dest.size()); + } + } +}; + +template<> struct gemv_dense_sense_selector<OnTheRight,RowMajor,true> +{ + template<typename Lhs, typename Rhs, typename Dest> + static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha) + { + typedef typename Dest::Index Index; + typedef typename Lhs::Scalar LhsScalar; + typedef typename Rhs::Scalar RhsScalar; + typedef typename Dest::Scalar ResScalar; + + typedef internal::blas_traits<Lhs> LhsBlasTraits; + typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType; + typedef internal::blas_traits<Rhs> RhsBlasTraits; + typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; + typedef typename internal::remove_all<ActualRhsType>::type ActualRhsTypeCleaned; + + typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(lhs); + typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(rhs); + + ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(lhs) + * RhsBlasTraits::extractScalarFactor(rhs); + + enum { + // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 + // on, the other hand it is good for the cache to pack the vector anyways... + DirectlyUseRhs = ActualRhsTypeCleaned::InnerStrideAtCompileTime==1 + }; + + gemv_static_vector_if<RhsScalar,ActualRhsTypeCleaned::SizeAtCompileTime,ActualRhsTypeCleaned::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs; + + ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(), + DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data()); + + if(!DirectlyUseRhs) + { + #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN + int size = actualRhs.size(); + EIGEN_DENSE_STORAGE_CTOR_PLUGIN + #endif + Map<typename ActualRhsTypeCleaned::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs; + } + + general_matrix_vector_product + <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( + actualLhs.rows(), actualLhs.cols(), + actualLhs.data(), actualLhs.outerStride(), + actualRhsPtr, 1, + dest.data(), dest.innerStride(), + actualAlpha); + } +}; + +template<> struct gemv_dense_sense_selector<OnTheRight,ColMajor,false> +{ + template<typename Lhs, typename Rhs, typename Dest> + static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha) + { + typedef typename Dest::Index Index; + // TODO makes sure dest is sequentially stored in memory, otherwise use a temp + const Index size = rhs.rows(); + for(Index k=0; k<size; ++k) + dest += (alpha*rhs.coeff(k)) * lhs.col(k); + } +}; + +template<> struct gemv_dense_sense_selector<OnTheRight,RowMajor,false> +{ + template<typename Lhs, typename Rhs, typename Dest> + static void run(const Lhs &lhs, const Rhs &rhs, Dest& dest, const typename Dest::Scalar& alpha) + { + typedef typename Dest::Index Index; + // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp + const Index rows = dest.rows(); + for(Index i=0; i<rows; ++i) + dest.coeffRef(i) += alpha * (lhs.row(i).cwiseProduct(rhs.transpose())).sum(); + } +}; + +#endif // EIGEN_ENABLE_EVALUATORS + } // end namespace internal /*************************************************************************** @@ -597,7 +814,7 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const return Product<Derived, OtherDerived>(derived(), other.derived()); } -#else +#else // EIGEN_TEST_EVALUATORS template<typename Derived> template<typename OtherDerived> inline const typename ProductReturnType<Derived, OtherDerived>::Type @@ -627,9 +844,10 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const #endif return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } -#endif +#endif // EIGEN_TEST_EVALUATORS + +#endif // __CUDACC__ -#endif /** \returns an expression of the matrix product of \c *this and \a other without implicit evaluation. * * The returned product will behave like any other expressions: the coefficients of the product will be @@ -641,6 +859,31 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const * * \sa operator*(const MatrixBase&) */ +#ifdef EIGEN_TEST_EVALUATORS +template<typename Derived> +template<typename OtherDerived> +const Product<Derived,OtherDerived,LazyProduct> +MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &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.cwiseProduct(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 Product<Derived,OtherDerived,LazyProduct>(derived(), other.derived()); +} +#else // EIGEN_TEST_EVALUATORS template<typename Derived> template<typename OtherDerived> const typename LazyProductReturnType<Derived,OtherDerived>::Type @@ -664,6 +907,7 @@ MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen |