Redux.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_REDUX_H
00027 #define EIGEN_REDUX_H
00028 
00029 namespace internal {
00030 
00031 // TODO
00032 //  * implement other kind of vectorization
00033 //  * factorize code
00034 
00035 /***************************************************************************
00036 * Part 1 : the logic deciding a strategy for vectorization and unrolling
00037 ***************************************************************************/
00038 
00039 template<typename Func, typename Derived>
00040 struct redux_traits
00041 {
00042 public:
00043   enum {
00044     PacketSize = packet_traits<typename Derived::Scalar>::size,
00045     InnerMaxSize = int(Derived::IsRowMajor)
00046                  ? Derived::MaxColsAtCompileTime
00047                  : Derived::MaxRowsAtCompileTime
00048   };
00049 
00050   enum {
00051     MightVectorize = (int(Derived::Flags)&ActualPacketAccessBit)
00052                   && (functor_traits<Func>::PacketAccess),
00053     MayLinearVectorize = MightVectorize && (int(Derived::Flags)&LinearAccessBit),
00054     MaySliceVectorize  = MightVectorize && int(InnerMaxSize)>=3*PacketSize
00055   };
00056 
00057 public:
00058   enum {
00059     Traversal = int(MayLinearVectorize) ? int(LinearVectorizedTraversal)
00060               : int(MaySliceVectorize)  ? int(SliceVectorizedTraversal)
00061                                         : int(DefaultTraversal)
00062   };
00063 
00064 public:
00065   enum {
00066     Cost = (  Derived::SizeAtCompileTime == Dynamic
00067            || Derived::CoeffReadCost == Dynamic
00068            || (Derived::SizeAtCompileTime!=1 && functor_traits<Func>::Cost == Dynamic)
00069            ) ? Dynamic
00070            : Derived::SizeAtCompileTime * Derived::CoeffReadCost
00071                + (Derived::SizeAtCompileTime-1) * functor_traits<Func>::Cost,
00072     UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize))
00073   };
00074 
00075 public:
00076   enum {
00077     Unrolling = Cost != Dynamic && Cost <= UnrollingLimit
00078               ? CompleteUnrolling
00079               : NoUnrolling
00080   };
00081 };
00082 
00083 /***************************************************************************
00084 * Part 2 : unrollers
00085 ***************************************************************************/
00086 
00087 /*** no vectorization ***/
00088 
00089 template<typename Func, typename Derived, int Start, int Length>
00090 struct redux_novec_unroller
00091 {
00092   enum {
00093     HalfLength = Length/2
00094   };
00095 
00096   typedef typename Derived::Scalar Scalar;
00097 
00098   EIGEN_STRONG_INLINE static Scalar run(const Derived &mat, const Func& func)
00099   {
00100     return func(redux_novec_unroller<Func, Derived, Start, HalfLength>::run(mat,func),
00101                 redux_novec_unroller<Func, Derived, Start+HalfLength, Length-HalfLength>::run(mat,func));
00102   }
00103 };
00104 
00105 template<typename Func, typename Derived, int Start>
00106 struct redux_novec_unroller<Func, Derived, Start, 1>
00107 {
00108   enum {
00109     outer = Start / Derived::InnerSizeAtCompileTime,
00110     inner = Start % Derived::InnerSizeAtCompileTime
00111   };
00112 
00113   typedef typename Derived::Scalar Scalar;
00114 
00115   EIGEN_STRONG_INLINE static Scalar run(const Derived &mat, const Func&)
00116   {
00117     return mat.coeffByOuterInner(outer, inner);
00118   }
00119 };
00120 
00121 // This is actually dead code and will never be called. It is required
00122 // to prevent false warnings regarding failed inlining though
00123 // for 0 length run() will never be called at all.
00124 template<typename Func, typename Derived, int Start>
00125 struct redux_novec_unroller<Func, Derived, Start, 0>
00126 {
00127   typedef typename Derived::Scalar Scalar;
00128   EIGEN_STRONG_INLINE static Scalar run(const Derived&, const Func&) { return Scalar(); }
00129 };
00130 
00131 /*** vectorization ***/
00132 
00133 template<typename Func, typename Derived, int Start, int Length>
00134 struct redux_vec_unroller
00135 {
00136   enum {
00137     PacketSize = packet_traits<typename Derived::Scalar>::size,
00138     HalfLength = Length/2
00139   };
00140 
00141   typedef typename Derived::Scalar Scalar;
00142   typedef typename packet_traits<Scalar>::type PacketScalar;
00143 
00144   EIGEN_STRONG_INLINE static PacketScalar run(const Derived &mat, const Func& func)
00145   {
00146     return func.packetOp(
00147             redux_vec_unroller<Func, Derived, Start, HalfLength>::run(mat,func),
00148             redux_vec_unroller<Func, Derived, Start+HalfLength, Length-HalfLength>::run(mat,func) );
00149   }
00150 };
00151 
00152 template<typename Func, typename Derived, int Start>
00153 struct redux_vec_unroller<Func, Derived, Start, 1>
00154 {
00155   enum {
00156     index = Start * packet_traits<typename Derived::Scalar>::size,
00157     outer = index / int(Derived::InnerSizeAtCompileTime),
00158     inner = index % int(Derived::InnerSizeAtCompileTime),
00159     alignment = (Derived::Flags & AlignedBit) ? Aligned : Unaligned
00160   };
00161 
00162   typedef typename Derived::Scalar Scalar;
00163   typedef typename packet_traits<Scalar>::type PacketScalar;
00164 
00165   EIGEN_STRONG_INLINE static PacketScalar run(const Derived &mat, const Func&)
00166   {
00167     return mat.template packetByOuterInner<alignment>(outer, inner);
00168   }
00169 };
00170 
00171 /***************************************************************************
00172 * Part 3 : implementation of all cases
00173 ***************************************************************************/
00174 
00175 template<typename Func, typename Derived,
00176          int Traversal = redux_traits<Func, Derived>::Traversal,
00177          int Unrolling = redux_traits<Func, Derived>::Unrolling
00178 >
00179 struct redux_impl;
00180 
00181 template<typename Func, typename Derived>
00182 struct redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>
00183 {
00184   typedef typename Derived::Scalar Scalar;
00185   typedef typename Derived::Index Index;
00186   static EIGEN_STRONG_INLINE Scalar run(const Derived& mat, const Func& func)
00187   {
00188     eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
00189     Scalar res;
00190     res = mat.coeffByOuterInner(0, 0);
00191     for(Index i = 1; i < mat.innerSize(); ++i)
00192       res = func(res, mat.coeffByOuterInner(0, i));
00193     for(Index i = 1; i < mat.outerSize(); ++i)
00194       for(Index j = 0; j < mat.innerSize(); ++j)
00195         res = func(res, mat.coeffByOuterInner(i, j));
00196     return res;
00197   }
00198 };
00199 
00200 template<typename Func, typename Derived>
00201 struct redux_impl<Func,Derived, DefaultTraversal, CompleteUnrolling>
00202   : public redux_novec_unroller<Func,Derived, 0, Derived::SizeAtCompileTime>
00203 {};
00204 
00205 template<typename Func, typename Derived>
00206 struct redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling>
00207 {
00208   typedef typename Derived::Scalar Scalar;
00209   typedef typename packet_traits<Scalar>::type PacketScalar;
00210   typedef typename Derived::Index Index;
00211 
00212   static Scalar run(const Derived& mat, const Func& func)
00213   {
00214     const Index size = mat.size();
00215     eigen_assert(size && "you are using an empty matrix");
00216     const Index packetSize = packet_traits<Scalar>::size;
00217     const Index alignedStart = first_aligned(mat);
00218     enum {
00219       alignment = bool(Derived::Flags & DirectAccessBit) || bool(Derived::Flags & AlignedBit)
00220                 ? Aligned : Unaligned
00221     };
00222     const Index alignedSize = ((size-alignedStart)/packetSize)*packetSize;
00223     const Index alignedEnd = alignedStart + alignedSize;
00224     Scalar res;
00225     if(alignedSize)
00226     {
00227       PacketScalar packet_res = mat.template packet<alignment>(alignedStart);
00228       for(Index index = alignedStart + packetSize; index < alignedEnd; index += packetSize)
00229         packet_res = func.packetOp(packet_res, mat.template packet<alignment>(index));
00230       res = func.predux(packet_res);
00231 
00232       for(Index index = 0; index < alignedStart; ++index)
00233         res = func(res,mat.coeff(index));
00234 
00235       for(Index index = alignedEnd; index < size; ++index)
00236         res = func(res,mat.coeff(index));
00237     }
00238     else // too small to vectorize anything.
00239          // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
00240     {
00241       res = mat.coeff(0);
00242       for(Index index = 1; index < size; ++index)
00243         res = func(res,mat.coeff(index));
00244     }
00245 
00246     return res;
00247   }
00248 };
00249 
00250 template<typename Func, typename Derived>
00251 struct redux_impl<Func, Derived, SliceVectorizedTraversal, NoUnrolling>
00252 {
00253   typedef typename Derived::Scalar Scalar;
00254   typedef typename packet_traits<Scalar>::type PacketScalar;
00255   typedef typename Derived::Index Index;
00256 
00257   static Scalar run(const Derived& mat, const Func& func)
00258   {
00259     eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
00260     const Index innerSize = mat.innerSize();
00261     const Index outerSize = mat.outerSize();
00262     enum {
00263       packetSize = packet_traits<Scalar>::size
00264     };
00265     const Index packetedInnerSize = ((innerSize)/packetSize)*packetSize;
00266     Scalar res;
00267     if(packetedInnerSize)
00268     {
00269       PacketScalar packet_res = mat.template packet<Unaligned>(0,0);
00270       for(Index j=0; j<outerSize; ++j)
00271         for(Index i=(j==0?packetSize:0); i<packetedInnerSize; i+=Index(packetSize))
00272           packet_res = func.packetOp(packet_res, mat.template packetByOuterInner<Unaligned>(j,i));
00273 
00274       res = func.predux(packet_res);
00275       for(Index j=0; j<outerSize; ++j)
00276         for(Index i=packetedInnerSize; i<innerSize; ++i)
00277           res = func(res, mat.coeffByOuterInner(j,i));
00278     }
00279     else // too small to vectorize anything.
00280          // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize.
00281     {
00282       res = redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>::run(mat, func);
00283     }
00284 
00285     return res;
00286   }
00287 };
00288 
00289 template<typename Func, typename Derived>
00290 struct redux_impl<Func, Derived, LinearVectorizedTraversal, CompleteUnrolling>
00291 {
00292   typedef typename Derived::Scalar Scalar;
00293   typedef typename packet_traits<Scalar>::type PacketScalar;
00294   enum {
00295     PacketSize = packet_traits<Scalar>::size,
00296     Size = Derived::SizeAtCompileTime,
00297     VectorizedSize = (Size / PacketSize) * PacketSize
00298   };
00299   EIGEN_STRONG_INLINE static Scalar run(const Derived& mat, const Func& func)
00300   {
00301     eigen_assert(mat.rows()>0 && mat.cols()>0 && "you are using an empty matrix");
00302     Scalar res = func.predux(redux_vec_unroller<Func, Derived, 0, Size / PacketSize>::run(mat,func));
00303     if (VectorizedSize != Size)
00304       res = func(res,redux_novec_unroller<Func, Derived, VectorizedSize, Size-VectorizedSize>::run(mat,func));
00305     return res;
00306   }
00307 };
00308 
00309 } // end namespace internal
00310 
00311 /***************************************************************************
00312 * Part 4 : public API
00313 ***************************************************************************/
00314 
00315 
00323 template<typename Derived>
00324 template<typename Func>
00325 EIGEN_STRONG_INLINE typename internal::result_of<Func(typename internal::traits<Derived>::Scalar)>::type
00326 DenseBase<Derived>::redux(const Func& func) const
00327 {
00328   typedef typename internal::remove_all<typename Derived::Nested>::type ThisNested;
00329   return internal::redux_impl<Func, ThisNested>
00330             ::run(derived(), func);
00331 }
00332 
00335 template<typename Derived>
00336 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00337 DenseBase<Derived>::minCoeff() const
00338 {
00339   return this->redux(Eigen::internal::scalar_min_op<Scalar>());
00340 }
00341 
00344 template<typename Derived>
00345 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00346 DenseBase<Derived>::maxCoeff() const
00347 {
00348   return this->redux(Eigen::internal::scalar_max_op<Scalar>());
00349 }
00350 
00355 template<typename Derived>
00356 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00357 DenseBase<Derived>::sum() const
00358 {
00359   if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
00360     return Scalar(0);
00361   return this->redux(Eigen::internal::scalar_sum_op<Scalar>());
00362 }
00363 
00368 template<typename Derived>
00369 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00370 DenseBase<Derived>::mean() const
00371 {
00372   return Scalar(this->redux(Eigen::internal::scalar_sum_op<Scalar>())) / Scalar(this->size());
00373 }
00374 
00382 template<typename Derived>
00383 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00384 DenseBase<Derived>::prod() const
00385 {
00386   if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
00387     return Scalar(1);
00388   return this->redux(Eigen::internal::scalar_product_op<Scalar>());
00389 }
00390 
00397 template<typename Derived>
00398 EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
00399 MatrixBase<Derived>::trace() const
00400 {
00401   return derived().diagonal().sum();
00402 }
00403 
00404 #endif // EIGEN_REDUX_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:18