// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2006-2008 Benoit Jacob // Copyright (C) 2008-2010 Gael Guennebaud // Copyright (C) 2011 Jitse Niesen // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_PRODUCTEVALUATORS_H #define EIGEN_PRODUCTEVALUATORS_H namespace Eigen { namespace internal { // We can evaluate the product either all at once, like GeneralProduct and its evalTo() function, or // traverse the matrix coefficient by coefficient, like CoeffBasedProduct. Use the existing logic // in ProductReturnType to decide. template struct product_evaluator_dispatcher; template struct evaluator_impl > : product_evaluator_dispatcher, typename ProductReturnType::Type> { typedef Product XprType; typedef product_evaluator_dispatcher::Type> Base; evaluator_impl(const XprType& xpr) : Base(xpr) { } }; template struct product_evaluator_traits_dispatcher; template struct evaluator_traits > : product_evaluator_traits_dispatcher, typename ProductReturnType::Type> { static const int AssumeAliasing = 1; }; // Case 1: Evaluate all at once // // We can view the GeneralProduct class as a part of the product evaluator. // Four sub-cases: InnerProduct, OuterProduct, GemmProduct and GemvProduct. // InnerProduct is special because GeneralProduct does not have an evalTo() method in this case. template struct product_evaluator_traits_dispatcher, GeneralProduct > { static const int HasEvalTo = 0; }; template struct product_evaluator_dispatcher, GeneralProduct > : public evaluator::PlainObject>::type { typedef Product XprType; typedef typename XprType::PlainObject PlainObject; typedef typename evaluator::type evaluator_base; // TODO: Computation is too early (?) product_evaluator_dispatcher(const XprType& xpr) : evaluator_base(m_result) { m_result.coeffRef(0,0) = (xpr.lhs().transpose().cwiseProduct(xpr.rhs())).sum(); } protected: PlainObject m_result; }; // For the other three subcases, simply call the evalTo() method of GeneralProduct // TODO: GeneralProduct should take evaluators, not expression objects. template struct product_evaluator_traits_dispatcher, GeneralProduct > { static const int HasEvalTo = 1; }; template struct product_evaluator_dispatcher, GeneralProduct > { typedef Product XprType; typedef typename XprType::PlainObject PlainObject; typedef typename evaluator::type evaluator_base; product_evaluator_dispatcher(const XprType& xpr) : m_xpr(xpr) { } template void evalTo(DstEvaluatorType /* not used */, DstXprType& dst) const { dst.resize(m_xpr.rows(), m_xpr.cols()); GeneralProduct(m_xpr.lhs(), m_xpr.rhs()).evalTo(dst); } protected: const XprType& m_xpr; }; // Case 2: Evaluate coeff by coeff // // This is mostly taken from CoeffBasedProduct.h // The main difference is that we add an extra argument to the etor_product_*_impl::run() function // for the inner dimension of the product, because evaluator object do not know their size. template struct etor_product_coeff_impl; template struct etor_product_packet_impl; template struct product_evaluator_traits_dispatcher, CoeffBasedProduct > { static const int HasEvalTo = 0; }; template struct product_evaluator_dispatcher, CoeffBasedProduct > : evaluator_impl_base > { typedef Product XprType; typedef CoeffBasedProduct CoeffBasedProductType; product_evaluator_dispatcher(const XprType& xpr) : m_lhsImpl(xpr.lhs()), m_rhsImpl(xpr.rhs()), m_innerDim(xpr.lhs().cols()) { } typedef typename XprType::Index Index; typedef typename XprType::Scalar Scalar; typedef typename XprType::CoeffReturnType CoeffReturnType; typedef typename XprType::PacketScalar PacketScalar; typedef typename XprType::PacketReturnType PacketReturnType; // Everything below here is taken from CoeffBasedProduct.h enum { RowsAtCompileTime = traits::RowsAtCompileTime, PacketSize = packet_traits::size, InnerSize = traits::InnerSize, CoeffReadCost = traits::CoeffReadCost, Unroll = CoeffReadCost != Dynamic && CoeffReadCost <= EIGEN_UNROLLING_LIMIT, CanVectorizeInner = traits::CanVectorizeInner }; typedef typename evaluator::type LhsEtorType; typedef typename evaluator::type RhsEtorType; typedef etor_product_coeff_impl CoeffImpl; const CoeffReturnType coeff(Index row, Index col) const { Scalar res; CoeffImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); return res; } /* Allow index-based non-packet access. It is impossible though to allow index-based packed access, * which is why we don't set the LinearAccessBit. */ const CoeffReturnType coeff(Index index) const { Scalar res; const Index row = RowsAtCompileTime == 1 ? 0 : index; const Index col = RowsAtCompileTime == 1 ? index : 0; CoeffImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); return res; } template const PacketReturnType packet(Index row, Index col) const { PacketScalar res; typedef etor_product_packet_impl PacketImpl; PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res); return res; } protected: typename evaluator::type m_lhsImpl; typename evaluator::type m_rhsImpl; // TODO: Get rid of m_innerDim if known at compile time Index m_innerDim; }; /*************************************************************************** * Normal product .coeff() implementation (with meta-unrolling) ***************************************************************************/ /************************************** *** Scalar path - no vectorization *** **************************************/ template struct etor_product_coeff_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar &res) { etor_product_coeff_impl::run(row, col, lhs, rhs, innerDim, res); res += lhs.coeff(row, UnrollingIndex) * rhs.coeff(UnrollingIndex, col); } }; template struct etor_product_coeff_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, RetScalar &res) { res = lhs.coeff(row, 0) * rhs.coeff(0, col); } }; template struct etor_product_coeff_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar& res) { eigen_assert(innerDim>0 && "you are using a non initialized matrix"); res = lhs.coeff(row, 0) * rhs.coeff(0, col); for(Index i = 1; i < innerDim; ++i) res += lhs.coeff(row, i) * rhs.coeff(i, col); } }; /******************************************* *** Scalar path with inner vectorization *** *******************************************/ template struct etor_product_coeff_vectorized_unroller { typedef typename Lhs::Index Index; enum { PacketSize = packet_traits::size }; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, typename Lhs::PacketScalar &pres) { etor_product_coeff_vectorized_unroller::run(row, col, lhs, rhs, innerDim, pres); pres = padd(pres, pmul( lhs.template packet(row, UnrollingIndex) , rhs.template packet(UnrollingIndex, col) )); } }; template struct etor_product_coeff_vectorized_unroller<0, Lhs, Rhs, Packet> { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::PacketScalar &pres) { pres = pmul(lhs.template packet(row, 0) , rhs.template packet(0, col)); } }; template struct etor_product_coeff_impl { typedef typename Lhs::PacketScalar Packet; typedef typename Lhs::Index Index; enum { PacketSize = packet_traits::size }; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, RetScalar &res) { Packet pres; etor_product_coeff_vectorized_unroller::run(row, col, lhs, rhs, innerDim, pres); etor_product_coeff_impl::run(row, col, lhs, rhs, innerDim, res); res = predux(pres); } }; template struct etor_product_coeff_vectorized_dyn_selector { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) { res = lhs.row(row).transpose().cwiseProduct(rhs.col(col)).sum(); } }; // NOTE the 3 following specializations are because taking .col(0) on a vector is a bit slower // NOTE maybe they are now useless since we have a specialization for Block template struct etor_product_coeff_vectorized_dyn_selector { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index /*row*/, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) { res = lhs.transpose().cwiseProduct(rhs.col(col)).sum(); } }; template struct etor_product_coeff_vectorized_dyn_selector { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index /*col*/, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) { res = lhs.row(row).transpose().cwiseProduct(rhs).sum(); } }; template struct etor_product_coeff_vectorized_dyn_selector { typedef typename Lhs::Index Index; EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, typename Lhs::Scalar &res) { res = lhs.transpose().cwiseProduct(rhs).sum(); } }; template struct etor_product_coeff_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, typename Lhs::Scalar &res) { etor_product_coeff_vectorized_dyn_selector::run(row, col, lhs, rhs, innerDim, res); } }; /******************* *** Packet path *** *******************/ template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res) { etor_product_packet_impl::run(row, col, lhs, rhs, innerDim, res); res = pmadd(pset1(lhs.coeff(row, UnrollingIndex)), rhs.template packet(UnrollingIndex, col), res); } }; template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res) { etor_product_packet_impl::run(row, col, lhs, rhs, innerDim, res); res = pmadd(lhs.template packet(row, UnrollingIndex), pset1(rhs.coeff(UnrollingIndex, col)), res); } }; template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res) { res = pmul(pset1(lhs.coeff(row, 0)),rhs.template packet(0, col)); } }; template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res) { res = pmul(lhs.template packet(row, 0), pset1(rhs.coeff(0, col))); } }; template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res) { eigen_assert(innerDim>0 && "you are using a non initialized matrix"); res = pmul(pset1(lhs.coeff(row, 0)),rhs.template packet(0, col)); for(Index i = 1; i < innerDim; ++i) res = pmadd(pset1(lhs.coeff(row, i)), rhs.template packet(i, col), res); } }; template struct etor_product_packet_impl { typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res) { eigen_assert(innerDim>0 && "you are using a non initialized matrix"); res = pmul(lhs.template packet(row, 0), pset1(rhs.coeff(0, col))); for(Index i = 1; i < innerDim; ++i) res = pmadd(lhs.template packet(row, i), pset1(rhs.coeff(i, col)), res); } }; } // end namespace internal } // end namespace Eigen #endif // EIGEN_PRODUCT_EVALUATORS_H