diff --git a/Eigen/Core b/Eigen/Core index e08bde4ab..b9c7e6c0f 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -48,6 +48,7 @@ namespace Eigen { #include "src/Core/Dot.h" #include "src/Core/DiagonalMatrix.h" #include "src/Core/DiagonalCoeffs.h" +#include "src/Core/Sum.h" #include "src/Core/Redux.h" #include "src/Core/Visitor.h" #include "src/Core/Fuzzy.h" diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 1c292d104..63eda1e85 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -31,19 +31,6 @@ * Part 1 : the logic deciding a strategy for vectorization and unrolling ***************************************************************************/ -enum { - NoVectorization, - InnerVectorization, - LinearVectorization, - SliceVectorization -}; - -enum { - CompleteUnrolling, - InnerUnrolling, - NoUnrolling -}; - template struct ei_assign_traits { diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 80b0e4ac7..cda47fa25 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -173,10 +173,13 @@ struct ei_functor_traits > template struct ei_scalar_abs2_op EIGEN_EMPTY_STRUCT { typedef typename NumTraits::Real result_type; inline const result_type operator() (const Scalar& a) const { return ei_abs2(a); } + template + inline const PacketScalar packetOp(const PacketScalar& a) const + { return ei_pmul(a,a); } }; template struct ei_functor_traits > -{ enum { Cost = NumTraits::MulCost, PacketAccess = false }; }; +{ enum { Cost = NumTraits::MulCost, PacketAccess = NumTraits::IsComplex==false && int(ei_packet_traits::size)>1 }; }; /** \internal * \brief Template functor to compute the conjugate of a complex value @@ -223,7 +226,7 @@ struct ei_functor_traits > * * \sa class CwiseUnaryOp, MatrixBase::operator*, MatrixBase::operator/ */ -template::size)>1?true:false) > struct ei_scalar_multiple_op; +template::size)>1) > struct ei_scalar_multiple_op; template struct ei_scalar_multiple_op { diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index 4b93e20fd..bede8362f 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h @@ -96,30 +96,6 @@ MatrixBase::redux(const BinaryOp& func) const ::run(derived(), func); } -/** \returns the sum of all coefficients of *this - * - * \sa trace() - */ -template -inline typename ei_traits::Scalar -MatrixBase::sum() const -{ - return this->redux(Eigen::ei_scalar_sum_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 -inline typename ei_traits::Scalar -MatrixBase::trace() const -{ - return diagonal().sum(); -} - /** \returns the minimum of all coefficients of *this */ template diff --git a/Eigen/src/Core/Sum.h b/Eigen/src/Core/Sum.h new file mode 100644 index 000000000..05a8722dc --- /dev/null +++ b/Eigen/src/Core/Sum.h @@ -0,0 +1,282 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2008 Benoit Jacob +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_SUM_H +#define EIGEN_SUM_H + +/*************************************************************************** +* Part 1 : the logic deciding a strategy for vectorization and unrolling +***************************************************************************/ + +template +struct ei_sum_traits +{ +public: + enum { + Vectorization = (int(Derived::Flags)&PacketAccessBit) + && (int(Derived::Flags)&LinearAccessBit) + ? LinearVectorization + : NoVectorization + }; + +private: + enum { + PacketSize = ei_packet_traits::size, + Cost = Derived::SizeAtCompileTime * Derived::CoeffReadCost + + (Derived::SizeAtCompileTime-1) * NumTraits::AddCost, + UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize)) + }; + +public: + enum { + Unrolling = Cost <= UnrollingLimit + ? CompleteUnrolling + : NoUnrolling + }; +}; + +/*************************************************************************** +* Part 2 : unrollers +***************************************************************************/ + +/*** no vectorization ***/ + +template +struct ei_sum_novec_unroller +{ + enum { + HalfLength = Length/2 + }; + + typedef typename Derived::Scalar Scalar; + + inline static Scalar run(const Derived &mat) + { + return ei_sum_novec_unroller::run(mat) + + ei_sum_novec_unroller::run(mat); + } +}; + +template +struct ei_sum_novec_unroller +{ + enum { + col = Start / Derived::RowsAtCompileTime, + row = Start % Derived::RowsAtCompileTime + }; + + typedef typename Derived::Scalar Scalar; + + inline static Scalar run(const Derived &mat) + { + return mat.coeff(row, col); + } +}; + +/*** vectorization ***/ + +template::size)> +struct ei_sum_vec_unroller +{ + enum { + row = int(Derived::Flags)&RowMajorBit + ? Index / int(Derived::ColsAtCompileTime) + : Index % Derived::RowsAtCompileTime, + col = int(Derived::Flags)&RowMajorBit + ? Index % int(Derived::ColsAtCompileTime) + : Index / Derived::RowsAtCompileTime + }; + + typedef typename Derived::Scalar Scalar; + typedef typename ei_packet_traits::type PacketScalar; + + inline static PacketScalar run(const Derived &mat) + { + return ei_padd( + mat.template packet(row, col), + ei_sum_vec_unroller::size, Stop>::run(mat) + ); + } +}; + +template +struct ei_sum_vec_unroller +{ + enum { + row = int(Derived::Flags)&RowMajorBit + ? Index / int(Derived::ColsAtCompileTime) + : Index % Derived::RowsAtCompileTime, + col = int(Derived::Flags)&RowMajorBit + ? Index % int(Derived::ColsAtCompileTime) + : Index / Derived::RowsAtCompileTime + }; + + typedef typename Derived::Scalar Scalar; + typedef typename ei_packet_traits::type PacketScalar; + + inline static PacketScalar run(const Derived &mat) + { + return mat.template packet(row, col); + } +}; + +/*************************************************************************** +* Part 3 : implementation of all cases +***************************************************************************/ + +template::Vectorization, + int Unrolling = ei_sum_traits::Unrolling +> +struct ei_sum_impl; + +template +struct ei_sum_impl +{ + typedef typename Derived::Scalar Scalar; + static Scalar run(const Derived& mat) + { + Scalar res; + res = mat.coeff(0, 0); + for(int i = 1; i < mat.rows(); i++) + res += mat.coeff(i, 0); + for(int j = 1; j < mat.cols(); j++) + for(int i = 0; i < mat.rows(); i++) + res += mat.coeff(i, j); + return res; + } +}; + +template +struct ei_sum_impl + : public ei_sum_novec_unroller +{}; + +template +struct ei_sum_impl +{ + typedef typename Derived::Scalar Scalar; + typedef typename ei_packet_traits::type PacketScalar; + + static Scalar run(const Derived& mat) + { + const int size = mat.size(); + const int packetSize = ei_packet_traits::size; + const int alignedSize = (size/packetSize)*packetSize; + const bool rowMajor = Derived::Flags&RowMajorBit; + const int innerSize = rowMajor ? mat.cols() : mat.rows(); + const int outerSize = rowMajor ? mat.rows() : mat.cols(); + Scalar res; + + // do the vectorizable part of the sum + if(size >= packetSize) + { + PacketScalar packet_res; + packet_res = mat.template packet(0, 0); + int index; + for(index = packetSize; index(row, col)); + } + res = ei_predux(packet_res); + + // now we must do the rest without vectorization. + if(alignedSize == size) return res; + } + else // too small to vectorize anything. + // since this is dynamic-size hence inefficient anyway, don't try to optimize. + { + res = Scalar(0); + } + + const int k = alignedSize/innerSize; + + // do the remainder of the current row or col + for(int i = alignedSize%innerSize; i < innerSize; i++) + { + const int row = rowMajor ? k : i; + const int col = rowMajor ? i : k; + res += mat.coeff(row, col); + } + + // do the remaining rows or cols + for(int j = k+1; j < outerSize; j++) + for(int i = 0; i < innerSize; i++) + { + const int row = rowMajor ? i : j; + const int col = rowMajor ? j : i; + res += mat.coeff(row, col); + } + + return res; + } +}; + +template +struct ei_sum_impl +{ + typedef typename Derived::Scalar Scalar; + static Scalar run(const Derived& mat) + { + return ei_predux( + ei_sum_vec_unroller::run(mat) + ); + } +}; + +/*************************************************************************** +* Part 4 : implementation of MatrixBase methods +***************************************************************************/ + +/** \returns the sum of all coefficients of *this + * + * \sa trace() + */ +template +inline typename ei_traits::Scalar +MatrixBase::sum() const +{ + return ei_sum_impl::run(derived()); +} + +/** \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 +inline typename ei_traits::Scalar +MatrixBase::trace() const +{ + return diagonal().sum(); +} + + +#endif // EIGEN_SUM_H diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 0b65a57aa..46d0a87b8 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -142,5 +142,18 @@ enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal }; enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct }; +enum { + NoVectorization, + InnerVectorization, + LinearVectorization, + SliceVectorization +}; + +enum { + CompleteUnrolling, + InnerUnrolling, + NoUnrolling +}; + #endif // EIGEN_CONSTANTS_H