File indexing completed on 2025-01-18 09:56:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #ifndef EIGEN_REDUX_H
0012 #define EIGEN_REDUX_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 template<typename Func, typename Evaluator>
0027 struct redux_traits
0028 {
0029 public:
0030 typedef typename find_best_packet<typename Evaluator::Scalar,Evaluator::SizeAtCompileTime>::type PacketType;
0031 enum {
0032 PacketSize = unpacket_traits<PacketType>::size,
0033 InnerMaxSize = int(Evaluator::IsRowMajor)
0034 ? Evaluator::MaxColsAtCompileTime
0035 : Evaluator::MaxRowsAtCompileTime,
0036 OuterMaxSize = int(Evaluator::IsRowMajor)
0037 ? Evaluator::MaxRowsAtCompileTime
0038 : Evaluator::MaxColsAtCompileTime,
0039 SliceVectorizedWork = int(InnerMaxSize)==Dynamic ? Dynamic
0040 : int(OuterMaxSize)==Dynamic ? (int(InnerMaxSize)>=int(PacketSize) ? Dynamic : 0)
0041 : (int(InnerMaxSize)/int(PacketSize)) * int(OuterMaxSize)
0042 };
0043
0044 enum {
0045 MightVectorize = (int(Evaluator::Flags)&ActualPacketAccessBit)
0046 && (functor_traits<Func>::PacketAccess),
0047 MayLinearVectorize = bool(MightVectorize) && (int(Evaluator::Flags)&LinearAccessBit),
0048 MaySliceVectorize = bool(MightVectorize) && (int(SliceVectorizedWork)==Dynamic || int(SliceVectorizedWork)>=3)
0049 };
0050
0051 public:
0052 enum {
0053 Traversal = int(MayLinearVectorize) ? int(LinearVectorizedTraversal)
0054 : int(MaySliceVectorize) ? int(SliceVectorizedTraversal)
0055 : int(DefaultTraversal)
0056 };
0057
0058 public:
0059 enum {
0060 Cost = Evaluator::SizeAtCompileTime == Dynamic ? HugeCost
0061 : int(Evaluator::SizeAtCompileTime) * int(Evaluator::CoeffReadCost) + (Evaluator::SizeAtCompileTime-1) * functor_traits<Func>::Cost,
0062 UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize))
0063 };
0064
0065 public:
0066 enum {
0067 Unrolling = Cost <= UnrollingLimit ? CompleteUnrolling : NoUnrolling
0068 };
0069
0070 #ifdef EIGEN_DEBUG_ASSIGN
0071 static void debug()
0072 {
0073 std::cerr << "Xpr: " << typeid(typename Evaluator::XprType).name() << std::endl;
0074 std::cerr.setf(std::ios::hex, std::ios::basefield);
0075 EIGEN_DEBUG_VAR(Evaluator::Flags)
0076 std::cerr.unsetf(std::ios::hex);
0077 EIGEN_DEBUG_VAR(InnerMaxSize)
0078 EIGEN_DEBUG_VAR(OuterMaxSize)
0079 EIGEN_DEBUG_VAR(SliceVectorizedWork)
0080 EIGEN_DEBUG_VAR(PacketSize)
0081 EIGEN_DEBUG_VAR(MightVectorize)
0082 EIGEN_DEBUG_VAR(MayLinearVectorize)
0083 EIGEN_DEBUG_VAR(MaySliceVectorize)
0084 std::cerr << "Traversal" << " = " << Traversal << " (" << demangle_traversal(Traversal) << ")" << std::endl;
0085 EIGEN_DEBUG_VAR(UnrollingLimit)
0086 std::cerr << "Unrolling" << " = " << Unrolling << " (" << demangle_unrolling(Unrolling) << ")" << std::endl;
0087 std::cerr << std::endl;
0088 }
0089 #endif
0090 };
0091
0092
0093
0094
0095
0096
0097
0098 template<typename Func, typename Evaluator, int Start, int Length>
0099 struct redux_novec_unroller
0100 {
0101 enum {
0102 HalfLength = Length/2
0103 };
0104
0105 typedef typename Evaluator::Scalar Scalar;
0106
0107 EIGEN_DEVICE_FUNC
0108 static EIGEN_STRONG_INLINE Scalar run(const Evaluator &eval, const Func& func)
0109 {
0110 return func(redux_novec_unroller<Func, Evaluator, Start, HalfLength>::run(eval,func),
0111 redux_novec_unroller<Func, Evaluator, Start+HalfLength, Length-HalfLength>::run(eval,func));
0112 }
0113 };
0114
0115 template<typename Func, typename Evaluator, int Start>
0116 struct redux_novec_unroller<Func, Evaluator, Start, 1>
0117 {
0118 enum {
0119 outer = Start / Evaluator::InnerSizeAtCompileTime,
0120 inner = Start % Evaluator::InnerSizeAtCompileTime
0121 };
0122
0123 typedef typename Evaluator::Scalar Scalar;
0124
0125 EIGEN_DEVICE_FUNC
0126 static EIGEN_STRONG_INLINE Scalar run(const Evaluator &eval, const Func&)
0127 {
0128 return eval.coeffByOuterInner(outer, inner);
0129 }
0130 };
0131
0132
0133
0134
0135 template<typename Func, typename Evaluator, int Start>
0136 struct redux_novec_unroller<Func, Evaluator, Start, 0>
0137 {
0138 typedef typename Evaluator::Scalar Scalar;
0139 EIGEN_DEVICE_FUNC
0140 static EIGEN_STRONG_INLINE Scalar run(const Evaluator&, const Func&) { return Scalar(); }
0141 };
0142
0143
0144
0145 template<typename Func, typename Evaluator, int Start, int Length>
0146 struct redux_vec_unroller
0147 {
0148 template<typename PacketType>
0149 EIGEN_DEVICE_FUNC
0150 static EIGEN_STRONG_INLINE PacketType run(const Evaluator &eval, const Func& func)
0151 {
0152 enum {
0153 PacketSize = unpacket_traits<PacketType>::size,
0154 HalfLength = Length/2
0155 };
0156
0157 return func.packetOp(
0158 redux_vec_unroller<Func, Evaluator, Start, HalfLength>::template run<PacketType>(eval,func),
0159 redux_vec_unroller<Func, Evaluator, Start+HalfLength, Length-HalfLength>::template run<PacketType>(eval,func) );
0160 }
0161 };
0162
0163 template<typename Func, typename Evaluator, int Start>
0164 struct redux_vec_unroller<Func, Evaluator, Start, 1>
0165 {
0166 template<typename PacketType>
0167 EIGEN_DEVICE_FUNC
0168 static EIGEN_STRONG_INLINE PacketType run(const Evaluator &eval, const Func&)
0169 {
0170 enum {
0171 PacketSize = unpacket_traits<PacketType>::size,
0172 index = Start * PacketSize,
0173 outer = index / int(Evaluator::InnerSizeAtCompileTime),
0174 inner = index % int(Evaluator::InnerSizeAtCompileTime),
0175 alignment = Evaluator::Alignment
0176 };
0177 return eval.template packetByOuterInner<alignment,PacketType>(outer, inner);
0178 }
0179 };
0180
0181
0182
0183
0184
0185 template<typename Func, typename Evaluator,
0186 int Traversal = redux_traits<Func, Evaluator>::Traversal,
0187 int Unrolling = redux_traits<Func, Evaluator>::Unrolling
0188 >
0189 struct redux_impl;
0190
0191 template<typename Func, typename Evaluator>
0192 struct redux_impl<Func, Evaluator, DefaultTraversal, NoUnrolling>
0193 {
0194 typedef typename Evaluator::Scalar Scalar;
0195
0196 template<typename XprType>
0197 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE
0198 Scalar run(const Evaluator &eval, const Func& func, const XprType& xpr)
0199 {
0200 eigen_assert(xpr.rows()>0 && xpr.cols()>0 && "you are using an empty matrix");
0201 Scalar res;
0202 res = eval.coeffByOuterInner(0, 0);
0203 for(Index i = 1; i < xpr.innerSize(); ++i)
0204 res = func(res, eval.coeffByOuterInner(0, i));
0205 for(Index i = 1; i < xpr.outerSize(); ++i)
0206 for(Index j = 0; j < xpr.innerSize(); ++j)
0207 res = func(res, eval.coeffByOuterInner(i, j));
0208 return res;
0209 }
0210 };
0211
0212 template<typename Func, typename Evaluator>
0213 struct redux_impl<Func,Evaluator, DefaultTraversal, CompleteUnrolling>
0214 : redux_novec_unroller<Func,Evaluator, 0, Evaluator::SizeAtCompileTime>
0215 {
0216 typedef redux_novec_unroller<Func,Evaluator, 0, Evaluator::SizeAtCompileTime> Base;
0217 typedef typename Evaluator::Scalar Scalar;
0218 template<typename XprType>
0219 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE
0220 Scalar run(const Evaluator &eval, const Func& func, const XprType& )
0221 {
0222 return Base::run(eval,func);
0223 }
0224 };
0225
0226 template<typename Func, typename Evaluator>
0227 struct redux_impl<Func, Evaluator, LinearVectorizedTraversal, NoUnrolling>
0228 {
0229 typedef typename Evaluator::Scalar Scalar;
0230 typedef typename redux_traits<Func, Evaluator>::PacketType PacketScalar;
0231
0232 template<typename XprType>
0233 static Scalar run(const Evaluator &eval, const Func& func, const XprType& xpr)
0234 {
0235 const Index size = xpr.size();
0236
0237 const Index packetSize = redux_traits<Func, Evaluator>::PacketSize;
0238 const int packetAlignment = unpacket_traits<PacketScalar>::alignment;
0239 enum {
0240 alignment0 = (bool(Evaluator::Flags & DirectAccessBit) && bool(packet_traits<Scalar>::AlignedOnScalar)) ? int(packetAlignment) : int(Unaligned),
0241 alignment = EIGEN_PLAIN_ENUM_MAX(alignment0, Evaluator::Alignment)
0242 };
0243 const Index alignedStart = internal::first_default_aligned(xpr);
0244 const Index alignedSize2 = ((size-alignedStart)/(2*packetSize))*(2*packetSize);
0245 const Index alignedSize = ((size-alignedStart)/(packetSize))*(packetSize);
0246 const Index alignedEnd2 = alignedStart + alignedSize2;
0247 const Index alignedEnd = alignedStart + alignedSize;
0248 Scalar res;
0249 if(alignedSize)
0250 {
0251 PacketScalar packet_res0 = eval.template packet<alignment,PacketScalar>(alignedStart);
0252 if(alignedSize>packetSize)
0253 {
0254 PacketScalar packet_res1 = eval.template packet<alignment,PacketScalar>(alignedStart+packetSize);
0255 for(Index index = alignedStart + 2*packetSize; index < alignedEnd2; index += 2*packetSize)
0256 {
0257 packet_res0 = func.packetOp(packet_res0, eval.template packet<alignment,PacketScalar>(index));
0258 packet_res1 = func.packetOp(packet_res1, eval.template packet<alignment,PacketScalar>(index+packetSize));
0259 }
0260
0261 packet_res0 = func.packetOp(packet_res0,packet_res1);
0262 if(alignedEnd>alignedEnd2)
0263 packet_res0 = func.packetOp(packet_res0, eval.template packet<alignment,PacketScalar>(alignedEnd2));
0264 }
0265 res = func.predux(packet_res0);
0266
0267 for(Index index = 0; index < alignedStart; ++index)
0268 res = func(res,eval.coeff(index));
0269
0270 for(Index index = alignedEnd; index < size; ++index)
0271 res = func(res,eval.coeff(index));
0272 }
0273 else
0274
0275 {
0276 res = eval.coeff(0);
0277 for(Index index = 1; index < size; ++index)
0278 res = func(res,eval.coeff(index));
0279 }
0280
0281 return res;
0282 }
0283 };
0284
0285
0286 template<typename Func, typename Evaluator, int Unrolling>
0287 struct redux_impl<Func, Evaluator, SliceVectorizedTraversal, Unrolling>
0288 {
0289 typedef typename Evaluator::Scalar Scalar;
0290 typedef typename redux_traits<Func, Evaluator>::PacketType PacketType;
0291
0292 template<typename XprType>
0293 EIGEN_DEVICE_FUNC static Scalar run(const Evaluator &eval, const Func& func, const XprType& xpr)
0294 {
0295 eigen_assert(xpr.rows()>0 && xpr.cols()>0 && "you are using an empty matrix");
0296 const Index innerSize = xpr.innerSize();
0297 const Index outerSize = xpr.outerSize();
0298 enum {
0299 packetSize = redux_traits<Func, Evaluator>::PacketSize
0300 };
0301 const Index packetedInnerSize = ((innerSize)/packetSize)*packetSize;
0302 Scalar res;
0303 if(packetedInnerSize)
0304 {
0305 PacketType packet_res = eval.template packet<Unaligned,PacketType>(0,0);
0306 for(Index j=0; j<outerSize; ++j)
0307 for(Index i=(j==0?packetSize:0); i<packetedInnerSize; i+=Index(packetSize))
0308 packet_res = func.packetOp(packet_res, eval.template packetByOuterInner<Unaligned,PacketType>(j,i));
0309
0310 res = func.predux(packet_res);
0311 for(Index j=0; j<outerSize; ++j)
0312 for(Index i=packetedInnerSize; i<innerSize; ++i)
0313 res = func(res, eval.coeffByOuterInner(j,i));
0314 }
0315 else
0316
0317 {
0318 res = redux_impl<Func, Evaluator, DefaultTraversal, NoUnrolling>::run(eval, func, xpr);
0319 }
0320
0321 return res;
0322 }
0323 };
0324
0325 template<typename Func, typename Evaluator>
0326 struct redux_impl<Func, Evaluator, LinearVectorizedTraversal, CompleteUnrolling>
0327 {
0328 typedef typename Evaluator::Scalar Scalar;
0329
0330 typedef typename redux_traits<Func, Evaluator>::PacketType PacketType;
0331 enum {
0332 PacketSize = redux_traits<Func, Evaluator>::PacketSize,
0333 Size = Evaluator::SizeAtCompileTime,
0334 VectorizedSize = (int(Size) / int(PacketSize)) * int(PacketSize)
0335 };
0336
0337 template<typename XprType>
0338 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE
0339 Scalar run(const Evaluator &eval, const Func& func, const XprType &xpr)
0340 {
0341 EIGEN_ONLY_USED_FOR_DEBUG(xpr)
0342 eigen_assert(xpr.rows()>0 && xpr.cols()>0 && "you are using an empty matrix");
0343 if (VectorizedSize > 0) {
0344 Scalar res = func.predux(redux_vec_unroller<Func, Evaluator, 0, Size / PacketSize>::template run<PacketType>(eval,func));
0345 if (VectorizedSize != Size)
0346 res = func(res,redux_novec_unroller<Func, Evaluator, VectorizedSize, Size-VectorizedSize>::run(eval,func));
0347 return res;
0348 }
0349 else {
0350 return redux_novec_unroller<Func, Evaluator, 0, Size>::run(eval,func);
0351 }
0352 }
0353 };
0354
0355
0356 template<typename _XprType>
0357 class redux_evaluator : public internal::evaluator<_XprType>
0358 {
0359 typedef internal::evaluator<_XprType> Base;
0360 public:
0361 typedef _XprType XprType;
0362 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0363 explicit redux_evaluator(const XprType &xpr) : Base(xpr) {}
0364
0365 typedef typename XprType::Scalar Scalar;
0366 typedef typename XprType::CoeffReturnType CoeffReturnType;
0367 typedef typename XprType::PacketScalar PacketScalar;
0368
0369 enum {
0370 MaxRowsAtCompileTime = XprType::MaxRowsAtCompileTime,
0371 MaxColsAtCompileTime = XprType::MaxColsAtCompileTime,
0372
0373 Flags = Base::Flags & ~DirectAccessBit,
0374 IsRowMajor = XprType::IsRowMajor,
0375 SizeAtCompileTime = XprType::SizeAtCompileTime,
0376 InnerSizeAtCompileTime = XprType::InnerSizeAtCompileTime
0377 };
0378
0379 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0380 CoeffReturnType coeffByOuterInner(Index outer, Index inner) const
0381 { return Base::coeff(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer); }
0382
0383 template<int LoadMode, typename PacketType>
0384 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0385 PacketType packetByOuterInner(Index outer, Index inner) const
0386 { return Base::template packet<LoadMode,PacketType>(IsRowMajor ? outer : inner, IsRowMajor ? inner : outer); }
0387
0388 };
0389
0390 }
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406 template<typename Derived>
0407 template<typename Func>
0408 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0409 DenseBase<Derived>::redux(const Func& func) const
0410 {
0411 eigen_assert(this->rows()>0 && this->cols()>0 && "you are using an empty matrix");
0412
0413 typedef typename internal::redux_evaluator<Derived> ThisEvaluator;
0414 ThisEvaluator thisEval(derived());
0415
0416
0417
0418 return internal::redux_impl<Func, ThisEvaluator>::run(thisEval, func, derived());
0419 }
0420
0421
0422
0423
0424
0425
0426
0427
0428 template<typename Derived>
0429 template<int NaNPropagation>
0430 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0431 DenseBase<Derived>::minCoeff() const
0432 {
0433 return derived().redux(Eigen::internal::scalar_min_op<Scalar,Scalar, NaNPropagation>());
0434 }
0435
0436
0437
0438
0439
0440
0441
0442
0443 template<typename Derived>
0444 template<int NaNPropagation>
0445 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0446 DenseBase<Derived>::maxCoeff() const
0447 {
0448 return derived().redux(Eigen::internal::scalar_max_op<Scalar,Scalar, NaNPropagation>());
0449 }
0450
0451
0452
0453
0454
0455
0456
0457 template<typename Derived>
0458 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0459 DenseBase<Derived>::sum() const
0460 {
0461 if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
0462 return Scalar(0);
0463 return derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>());
0464 }
0465
0466
0467
0468
0469
0470 template<typename Derived>
0471 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0472 DenseBase<Derived>::mean() const
0473 {
0474 #ifdef __INTEL_COMPILER
0475 #pragma warning push
0476 #pragma warning ( disable : 2259 )
0477 #endif
0478 return Scalar(derived().redux(Eigen::internal::scalar_sum_op<Scalar,Scalar>())) / Scalar(this->size());
0479 #ifdef __INTEL_COMPILER
0480 #pragma warning pop
0481 #endif
0482 }
0483
0484
0485
0486
0487
0488
0489
0490
0491 template<typename Derived>
0492 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0493 DenseBase<Derived>::prod() const
0494 {
0495 if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0))
0496 return Scalar(1);
0497 return derived().redux(Eigen::internal::scalar_product_op<Scalar>());
0498 }
0499
0500
0501
0502
0503
0504
0505
0506 template<typename Derived>
0507 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar
0508 MatrixBase<Derived>::trace() const
0509 {
0510 return derived().diagonal().sum();
0511 }
0512
0513 }
0514
0515 #endif