// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud // Copyright (C) 2006-2008 Benoit Jacob // // 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_REDUX_H #define EIGEN_REDUX_H namespace Eigen { namespace internal { // TODO // * implement other kind of vectorization // * factorize code /*************************************************************************** * Part 1 : the logic deciding a strategy for vectorization and unrolling ***************************************************************************/ template struct redux_traits { public: enum { PacketSize = packet_traits::size, InnerMaxSize = int(Derived::IsRowMajor) ? Derived::MaxColsAtCompileTime : Derived::MaxRowsAtCompileTime }; enum { MightVectorize = (int(Derived::Flags)&ActualPacketAccessBit) && (functor_traits::PacketAccess), MayLinearVectorize = MightVectorize && (int(Derived::Flags)&LinearAccessBit), MaySliceVectorize = MightVectorize && int(InnerMaxSize)>=3*PacketSize }; public: enum { Traversal = int(MayLinearVectorize) ? int(LinearVectorizedTraversal) : int(MaySliceVectorize) ? int(SliceVectorizedTraversal) : int(DefaultTraversal) }; public: enum { Cost = ( Derived::SizeAtCompileTime == Dynamic || Derived::CoeffReadCost == Dynamic || (Derived::SizeAtCompileTime!=1 && functor_traits::Cost == Dynamic) ) ? Dynamic : Derived::SizeAtCompileTime * Derived::CoeffReadCost + (Derived::SizeAtCompileTime-1) * functor_traits::Cost, UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize)) }; public: enum { Unrolling = Cost != Dynamic && Cost <= UnrollingLimit ? CompleteUnrolling : NoUnrolling }; }; /*************************************************************************** * Part 2 : unrollers ***************************************************************************/ /*** no vectorization ***/ template struct redux_novec_unroller { enum { HalfLength = Length/2 }; typedef typename Derived::Scalar Scalar; static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func& func) { return func(redux_novec_unroller::run(mat,func), redux_novec_unroller::run(mat,func)); } }; template struct redux_novec_unroller { enum { outer = Start / Derived::InnerSizeAtCompileTime, inner = Start % Derived::InnerSizeAtCompileTime }; typedef typename Derived::Scalar Scalar; static EIGEN_STRONG_INLINE Scalar run(const Derived &mat, const Func&) { return mat.coeffByOuterInner(outer, inner); } }; // This is actually dead code and will never be called. It is required // to prevent false warnings regarding failed inlining though // for 0 length run() will never be called at all. template struct redux_novec_unroller { typedef typename Derived::Scalar Scalar; static EIGEN_STRONG_INLINE Scalar run(const Derived&, const Func&) { return Scalar(); } }; /*** vectorization ***/ template struct redux_vec_unroller { enum { PacketSize = packet_traits::size, HalfLength = Length/2 }; typedef typename Derived::Scalar Scalar; typedef typename packet_traits::type PacketScalar; static EIGEN_STRONG_INLINE PacketScalar run(const Derived &mat, const Func& func) { return func.packetOp( redux_vec_unroller::run(mat,func), redux_vec_unroller::run(mat,func) ); } }; template struct redux_vec_unroller { enum { index = Start * packet_traits::size, outer = index / int(Derived::InnerSizeAtCompileTime), inner = index % int(Derived::InnerSizeAtCompileTime), alignment = (Derived::Flags & AlignedBit) ? Aligned : Unaligned }; typedef typename Derived::Scalar Scalar; typedef typename packet_traits::type PacketScalar; static EIGEN_STRONG_INLINE PacketScalar run(const Derived &mat, const Func&) { return mat.template packetByOuterInner(outer, inner); } }; /*************************************************************************** * Part 3 : implementation of all cases ***************************************************************************/ template::Traversal, int Unrolling = redux_traits::Unrolling > struct redux_impl; template struct redux_impl { typedef typename Derived::Scalar Scalar; typedef typename Derived::Index Index; static EIGEN_STRONG_INLINE Scalar run(const Derived& mat, const Func& func) { eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix"); Scalar res; res = mat.coeffByOuterInner(0, 0); for(Index i = 1; i < mat.innerSize(); ++i) res = func(res, mat.coeffByOuterInner(0, i)); for(Index i = 1; i < mat.outerSize(); ++i) for(Index j = 0; j < mat.innerSize(); ++j) res = func(res, mat.coeffByOuterInner(i, j)); return res; } }; template struct redux_impl : public redux_novec_unroller {}; template struct redux_impl { typedef typename Derived::Scalar Scalar; typedef typename packet_traits::type PacketScalar; typedef typename Derived::Index Index; static Scalar run(const Derived& mat, const Func& func) { const Index size = mat.size(); eigen_assert(size && "you are using an empty matrix"); const Index packetSize = packet_traits::size; const Index alignedStart = internal::first_aligned(mat); enum { alignment = bool(Derived::Flags & DirectAccessBit) || bool(Derived::Flags & AlignedBit) ? Aligned : Unaligned }; const Index alignedSize2 = ((size-alignedStart)/(2*packetSize))*(2*packetSize); const Index alignedSize = ((size-alignedStart)/(packetSize))*(packetSize); const Index alignedEnd2 = alignedStart + alignedSize2; const Index alignedEnd = alignedStart + alignedSize; Scalar res; if(alignedSize) { PacketScalar packet_res0 = mat.template packet(alignedStart); if(alignedSize>packetSize) // we have at least two packets to partly unroll the loop { PacketScalar packet_res1 = mat.template packet(alignedStart+packetSize); for(Index index = alignedStart + 2*packetSize; index < alignedEnd2; index += 2*packetSize) { packet_res0 = func.packetOp(packet_res0, mat.template packet(index)); packet_res1 = func.packetOp(packet_res1, mat.template packet(index+packetSize)); } packet_res0 = func.packetOp(packet_res0,packet_res1); if(alignedEnd>alignedEnd2) packet_res0 = func.packetOp(packet_res0, mat.template packet(alignedEnd2)); } res = func.predux(packet_res0); for(Index index = 0; index < alignedStart; ++index) res = func(res,mat.coeff(index)); for(Index index = alignedEnd; index < size; ++index) res = func(res,mat.coeff(index)); } else // too small to vectorize anything. // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize. { res = mat.coeff(0); for(Index index = 1; index < size; ++index) res = func(res,mat.coeff(index)); } return res; } }; template struct redux_impl { typedef typename Derived::Scalar Scalar; typedef typename packet_traits::type PacketScalar; typedef typename Derived::Index Index; static Scalar run(const Derived& mat, const Func& func) { eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix"); const Index innerSize = mat.innerSize(); const Index outerSize = mat.outerSize(); enum { packetSize = packet_traits::size }; const Index packetedInnerSize = ((innerSize)/packetSize)*packetSize; Scalar res; if(packetedInnerSize) { PacketScalar packet_res = mat.template packet(0,0); for(Index j=0; j(j,i)); res = func.predux(packet_res); for(Index j=0; j::run(mat, func); } return res; } }; template struct redux_impl { typedef typename Derived::Scalar Scalar; typedef typename packet_traits::type PacketScalar; enum { PacketSize = packet_traits::size, Size = Derived::SizeAtCompileTime, VectorizedSize = (Size / PacketSize) * PacketSize }; static EIGEN_STRONG_INLINE Scalar run(const Derived& mat, const Func& func) { eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix"); Scalar res = func.predux(redux_vec_unroller::run(mat,func)); if (VectorizedSize != Size) res = func(res,redux_novec_unroller::run(mat,func)); return res; } }; } // end namespace internal /*************************************************************************** * Part 4 : public API ***************************************************************************/ /** \returns the result of a full redux operation on the whole matrix or vector using \a func * * The template parameter \a BinaryOp is the type of the functor \a func which must be * an associative operator. Both current STL and TR1 functor styles are handled. * * \sa DenseBase::sum(), DenseBase::minCoeff(), DenseBase::maxCoeff(), MatrixBase::colwise(), MatrixBase::rowwise() */ template template EIGEN_STRONG_INLINE typename internal::result_of::Scalar)>::type DenseBase::redux(const Func& func) const { typedef typename internal::remove_all::type ThisNested; return internal::redux_impl ::run(derived(), func); } /** \returns the minimum of all coefficients of \c *this. * \warning the result is undefined if \c *this contains NaN. */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar DenseBase::minCoeff() const { return this->redux(Eigen::internal::scalar_min_op()); } /** \returns the maximum of all coefficients of \c *this. * \warning the result is undefined if \c *this contains NaN. */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar DenseBase::maxCoeff() const { return this->redux(Eigen::internal::scalar_max_op()); } /** \returns the sum of all coefficients of *this * * \sa trace(), prod(), mean() */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar DenseBase::sum() const { if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0)) return Scalar(0); return this->redux(Eigen::internal::scalar_sum_op()); } /** \returns the mean of all coefficients of *this * * \sa trace(), prod(), sum() */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar DenseBase::mean() const { return Scalar(this->redux(Eigen::internal::scalar_sum_op())) / Scalar(this->size()); } /** \returns the product of all coefficients of *this * * Example: \include MatrixBase_prod.cpp * Output: \verbinclude MatrixBase_prod.out * * \sa sum(), mean(), trace() */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar DenseBase::prod() const { if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0)) return Scalar(1); return this->redux(Eigen::internal::scalar_product_op()); } /** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal. * * \c *this can be any matrix, not necessarily square. * * \sa diagonal(), sum() */ template EIGEN_STRONG_INLINE typename internal::traits::Scalar MatrixBase::trace() const { return derived().diagonal().sum(); } } // end namespace Eigen #endif // EIGEN_REDUX_H