1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com>
5 // Copyright (C) 2016 Mehdi Goli, Codeplay Software Ltd <eigen@codeplay.com>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_CXX11_TENSOR_TENSOR_REDUCTION_H
12 #define EIGEN_CXX11_TENSOR_TENSOR_REDUCTION_H
13 
14 namespace Eigen {
15 
16 /** \class TensorReduction
17   * \ingroup CXX11_Tensor_Module
18   *
19   * \brief Tensor reduction class.
20   *
21   */
22 
23 namespace internal {
24   template<typename Op, typename Dims, typename XprType,template <class> class MakePointer_ >
25   struct traits<TensorReductionOp<Op, Dims, XprType, MakePointer_> >
26  : traits<XprType>
27 {
28   typedef traits<XprType> XprTraits;
29   typedef typename XprTraits::Scalar Scalar;
30   typedef typename XprTraits::StorageKind StorageKind;
31   typedef typename XprTraits::Index Index;
32   typedef typename XprType::Nested Nested;
33   static const int NumDimensions = XprTraits::NumDimensions - array_size<Dims>::value;
34   static const int Layout = XprTraits::Layout;
35 
36   template <class T> struct MakePointer {
37     // Intermediate typedef to workaround MSVC issue.
38     typedef MakePointer_<T> MakePointerT;
39     typedef typename MakePointerT::Type Type;
40   };
41 };
42 
43 template<typename Op, typename Dims, typename XprType, template <class> class MakePointer_>
44 struct eval<TensorReductionOp<Op, Dims, XprType, MakePointer_>, Eigen::Dense>
45 {
46   typedef const TensorReductionOp<Op, Dims, XprType, MakePointer_>& type;
47 };
48 
49 template<typename Op, typename Dims, typename XprType, template <class> class MakePointer_>
50 struct nested<TensorReductionOp<Op, Dims, XprType, MakePointer_>, 1, typename eval<TensorReductionOp<Op, Dims, XprType, MakePointer_> >::type>
51 {
52   typedef TensorReductionOp<Op, Dims, XprType, MakePointer_> type;
53 };
54 
55 
56 template <typename OutputDims> struct DimInitializer {
57   template <typename InputDims, typename ReducedDims> EIGEN_DEVICE_FUNC
58   static void run(const InputDims& input_dims,
59                   const array<bool, internal::array_size<InputDims>::value>& reduced,
60                   OutputDims* output_dims, ReducedDims* reduced_dims) {
61     const int NumInputDims = internal::array_size<InputDims>::value;
62     int outputIndex = 0;
63     int reduceIndex = 0;
64     for (int i = 0; i < NumInputDims; ++i) {
65       if (reduced[i]) {
66         (*reduced_dims)[reduceIndex] = input_dims[i];
67         ++reduceIndex;
68       } else {
69         (*output_dims)[outputIndex] = input_dims[i];
70         ++outputIndex;
71       }
72     }
73   }
74 };
75 
76 template <> struct DimInitializer<Sizes<> > {
77   template <typename InputDims, typename Index, size_t Rank> EIGEN_DEVICE_FUNC
78   static void run(const InputDims& input_dims, const array<bool, Rank>&,
79                   Sizes<>*, array<Index, Rank>* reduced_dims) {
80     const int NumInputDims = internal::array_size<InputDims>::value;
81     for (int i = 0; i < NumInputDims; ++i) {
82       (*reduced_dims)[i] = input_dims[i];
83     }
84   }
85 };
86 
87 
88 template <typename ReducedDims, int NumTensorDims, int Layout>
89 struct are_inner_most_dims {
90   static const bool value = false;
91 };
92 template <typename ReducedDims, int NumTensorDims, int Layout>
93 struct preserve_inner_most_dims {
94   static const bool value = false;
95 };
96 
97 #if EIGEN_HAS_CONSTEXPR && EIGEN_HAS_VARIADIC_TEMPLATES
98 template <typename ReducedDims, int NumTensorDims>
99 struct are_inner_most_dims<ReducedDims, NumTensorDims, ColMajor>{
100   static const bool tmp1 = indices_statically_known_to_increase<ReducedDims>();
101   static const bool tmp2 = index_statically_eq<ReducedDims>(0, 0);
102   static const bool tmp3 = index_statically_eq<ReducedDims>(array_size<ReducedDims>::value-1, array_size<ReducedDims>::value-1);
103   static const bool value = tmp1 & tmp2 & tmp3;
104 };
105 template <typename ReducedDims, int NumTensorDims>
106 struct are_inner_most_dims<ReducedDims, NumTensorDims, RowMajor>{
107   static const bool tmp1 = indices_statically_known_to_increase<ReducedDims>();
108   static const bool tmp2 = index_statically_eq<ReducedDims>(0, NumTensorDims - array_size<ReducedDims>::value);
109   static const bool tmp3 = index_statically_eq<ReducedDims>(array_size<ReducedDims>::value - 1, NumTensorDims - 1);
110   static const bool value = tmp1 & tmp2 & tmp3;
111 
112 };
113 template <typename ReducedDims, int NumTensorDims>
114 struct preserve_inner_most_dims<ReducedDims, NumTensorDims, ColMajor>{
115   static const bool tmp1 = indices_statically_known_to_increase<ReducedDims>();
116   static const bool tmp2 = index_statically_gt<ReducedDims>(0, 0);
117   static const bool value = tmp1 & tmp2;
118 
119 };
120 template <typename ReducedDims, int NumTensorDims>
121 struct preserve_inner_most_dims<ReducedDims, NumTensorDims, RowMajor>{
122   static const bool tmp1 = indices_statically_known_to_increase<ReducedDims>();
123   static const bool tmp2 = index_statically_lt<ReducedDims>(array_size<ReducedDims>::value - 1, NumTensorDims - 1);
124   static const bool value = tmp1 & tmp2;
125 };
126 #endif
127 
128 
129 template <int DimIndex, typename Self, typename Op>
130 struct GenericDimReducer {
131   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self& self, typename Self::Index firstIndex, Op& reducer, typename Self::CoeffReturnType* accum) {
132     EIGEN_STATIC_ASSERT((DimIndex > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
133     for (int j = 0; j < self.m_reducedDims[DimIndex]; ++j) {
134       const typename Self::Index input = firstIndex + j * self.m_reducedStrides[DimIndex];
135       GenericDimReducer<DimIndex-1, Self, Op>::reduce(self, input, reducer, accum);
136     }
137   }
138 };
139 template <typename Self, typename Op>
140 struct GenericDimReducer<0, Self, Op> {
141   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self& self, typename Self::Index firstIndex, Op& reducer, typename Self::CoeffReturnType* accum) {
142     for (int j = 0; j < self.m_reducedDims[0]; ++j) {
143       const typename Self::Index input = firstIndex + j * self.m_reducedStrides[0];
144       reducer.reduce(self.m_impl.coeff(input), accum);
145     }
146   }
147 };
148 template <typename Self, typename Op>
149 struct GenericDimReducer<-1, Self, Op> {
150   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self& self, typename Self::Index index, Op& reducer, typename Self::CoeffReturnType* accum) {
151     reducer.reduce(self.m_impl.coeff(index), accum);
152   }
153 };
154 
155 template <typename Self, typename Op, bool Vectorizable = (Self::InputPacketAccess & Op::PacketAccess)>
156 struct InnerMostDimReducer {
157   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Self::CoeffReturnType reduce(const Self& self, typename Self::Index firstIndex, typename Self::Index numValuesToReduce, Op& reducer) {
158     typename Self::CoeffReturnType accum = reducer.initialize();
159     for (typename Self::Index j = 0; j < numValuesToReduce; ++j) {
160       reducer.reduce(self.m_impl.coeff(firstIndex + j), &accum);
161     }
162     return reducer.finalize(accum);
163   }
164 };
165 
166 template <typename Self, typename Op>
167 struct InnerMostDimReducer<Self, Op, true> {
168   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Self::CoeffReturnType reduce(const Self& self, typename Self::Index firstIndex, typename Self::Index numValuesToReduce, Op& reducer) {
169     const int packetSize = internal::unpacket_traits<typename Self::PacketReturnType>::size;
170     const typename Self::Index VectorizedSize = (numValuesToReduce / packetSize) * packetSize;
171     typename Self::PacketReturnType p = reducer.template initializePacket<typename Self::PacketReturnType>();
172     for (typename Self::Index j = 0; j < VectorizedSize; j += packetSize) {
173       reducer.reducePacket(self.m_impl.template packet<Unaligned>(firstIndex + j), &p);
174     }
175     typename Self::CoeffReturnType accum = reducer.initialize();
176     for (typename Self::Index j = VectorizedSize; j < numValuesToReduce; ++j) {
177       reducer.reduce(self.m_impl.coeff(firstIndex + j), &accum);
178     }
179     return reducer.finalizeBoth(accum, p);
180   }
181 };
182 
183 template <int DimIndex, typename Self, typename Op, bool vectorizable = (Self::InputPacketAccess & Op::PacketAccess)>
184 struct InnerMostDimPreserver {
185   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self&, typename Self::Index, Op&, typename Self::PacketReturnType*) {
186     eigen_assert(false && "should never be called");
187   }
188 };
189 
190 template <int DimIndex, typename Self, typename Op>
191 struct InnerMostDimPreserver<DimIndex, Self, Op, true> {
192   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self& self, typename Self::Index firstIndex, Op& reducer, typename Self::PacketReturnType* accum) {
193     EIGEN_STATIC_ASSERT((DimIndex > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
194     for (typename Self::Index j = 0; j < self.m_reducedDims[DimIndex]; ++j) {
195       const typename Self::Index input = firstIndex + j * self.m_reducedStrides[DimIndex];
196       InnerMostDimPreserver<DimIndex-1, Self, Op>::reduce(self, input, reducer, accum);
197     }
198   }
199 };
200 
201 template <typename Self, typename Op>
202 struct InnerMostDimPreserver<0, Self, Op, true> {
203   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self& self, typename Self::Index firstIndex, Op& reducer, typename Self::PacketReturnType* accum) {
204     for (typename Self::Index j = 0; j < self.m_reducedDims[0]; ++j) {
205       const typename Self::Index input = firstIndex + j * self.m_reducedStrides[0];
206       reducer.reducePacket(self.m_impl.template packet<Unaligned>(input), accum);
207     }
208   }
209 };
210 template <typename Self, typename Op>
211 struct InnerMostDimPreserver<-1, Self, Op, true> {
212   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void reduce(const Self&, typename Self::Index, Op&, typename Self::PacketReturnType*) {
213     eigen_assert(false && "should never be called");
214   }
215 };
216 
217 // Default full reducer
218 template <typename Self, typename Op, typename Device, bool Vectorizable = (Self::InputPacketAccess & Op::PacketAccess)>
219 struct FullReducer {
220   static const bool HasOptimizedImplementation = false;
221 
222   static EIGEN_DEVICE_FUNC void run(const Self& self, Op& reducer, const Device&, typename Self::CoeffReturnType* output) {
223     const typename Self::Index num_coeffs = array_prod(self.m_impl.dimensions());
224     *output = InnerMostDimReducer<Self, Op, Vectorizable>::reduce(self, 0, num_coeffs, reducer);
225   }
226 };
227 
228 
229 #ifdef EIGEN_USE_THREADS
230 // Multithreaded full reducers
231 template <typename Self, typename Op,
232           bool Vectorizable = (Self::InputPacketAccess & Op::PacketAccess)>
233 struct FullReducerShard {
234   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(const Self& self, typename Self::Index firstIndex,
235                   typename Self::Index numValuesToReduce, Op& reducer,
236                   typename Self::CoeffReturnType* output) {
237     *output = InnerMostDimReducer<Self, Op, Vectorizable>::reduce(
238         self, firstIndex, numValuesToReduce, reducer);
239   }
240 };
241 
242 // Multithreaded full reducer
243 template <typename Self, typename Op, bool Vectorizable>
244 struct FullReducer<Self, Op, ThreadPoolDevice, Vectorizable> {
245   static const bool HasOptimizedImplementation = !Op::IsStateful;
246   static const int PacketSize =
247       unpacket_traits<typename Self::PacketReturnType>::size;
248 
249   // launch one reducer per thread and accumulate the result.
250   static void run(const Self& self, Op& reducer, const ThreadPoolDevice& device,
251                   typename Self::CoeffReturnType* output) {
252     typedef typename Self::Index Index;
253     const Index num_coeffs = array_prod(self.m_impl.dimensions());
254     if (num_coeffs == 0) {
255       *output = reducer.finalize(reducer.initialize());
256       return;
257     }
258     const TensorOpCost cost =
259         self.m_impl.costPerCoeff(Vectorizable) +
260         TensorOpCost(0, 0, internal::functor_traits<Op>::Cost, Vectorizable,
261                      PacketSize);
262     const int num_threads = TensorCostModel<ThreadPoolDevice>::numThreads(
263         num_coeffs, cost, device.numThreads());
264     if (num_threads == 1) {
265       *output =
266           InnerMostDimReducer<Self, Op, Vectorizable>::reduce(self, 0, num_coeffs, reducer);
267       return;
268     }
269     const Index blocksize =
270         std::floor<Index>(static_cast<float>(num_coeffs) / num_threads);
271     const Index numblocks = blocksize > 0 ? num_coeffs / blocksize : 0;
272     eigen_assert(num_coeffs >= numblocks * blocksize);
273 
274     Barrier barrier(internal::convert_index<unsigned int>(numblocks));
275     MaxSizeVector<typename Self::CoeffReturnType> shards(numblocks, reducer.initialize());
276     for (Index i = 0; i < numblocks; ++i) {
277       device.enqueue_with_barrier(&barrier, &FullReducerShard<Self, Op, Vectorizable>::run,
278                                   self, i * blocksize, blocksize, reducer,
279                                   &shards[i]);
280     }
281     typename Self::CoeffReturnType finalShard;
282     if (numblocks * blocksize < num_coeffs) {
283       finalShard = InnerMostDimReducer<Self, Op, Vectorizable>::reduce(
284           self, numblocks * blocksize, num_coeffs - numblocks * blocksize,
285           reducer);
286     } else {
287       finalShard = reducer.initialize();
288     }
289     barrier.Wait();
290 
291     for (Index i = 0; i < numblocks; ++i) {
292       reducer.reduce(shards[i], &finalShard);
293     }
294     *output = reducer.finalize(finalShard);
295   }
296 };
297 
298 #endif
299 
300 
301 // Default inner reducer
302 template <typename Self, typename Op, typename Device>
303 struct InnerReducer {
304   static const bool HasOptimizedImplementation = false;
305 
306   EIGEN_DEVICE_FUNC static bool run(const Self&, Op&, const Device&, typename Self::CoeffReturnType*, typename Self::Index, typename Self::Index) {
307     eigen_assert(false && "Not implemented");
308     return true;
309   }
310 };
311 
312 // Default outer reducer
313 template <typename Self, typename Op, typename Device>
314 struct OuterReducer {
315   static const bool HasOptimizedImplementation = false;
316 
317   EIGEN_DEVICE_FUNC static bool run(const Self&, Op&, const Device&, typename Self::CoeffReturnType*, typename Self::Index, typename Self::Index) {
318     eigen_assert(false && "Not implemented");
319     return true;
320   }
321 };
322 
323 
324 #if defined(EIGEN_USE_GPU) && defined(__CUDACC__)
325 template <int B, int N, typename S, typename R, typename I>
326 __global__ void FullReductionKernel(R, const S, I, typename S::CoeffReturnType*, unsigned int*);
327 
328 
329 #ifdef EIGEN_HAS_CUDA_FP16
330 template <typename S, typename R, typename I>
331 __global__ void ReductionInitFullReduxKernelHalfFloat(R, const S, I, half2*);
332 template <int B, int N, typename S, typename R, typename I>
333 __global__ void FullReductionKernelHalfFloat(R, const S, I, half*, half2*);
334 template <int NPT, typename S, typename R, typename I>
335 __global__ void InnerReductionKernelHalfFloat(R, const S, I, I, half*);
336 
337 #endif
338 
339 template <int NPT, typename S, typename R, typename I>
340 __global__ void InnerReductionKernel(R, const S, I, I, typename S::CoeffReturnType*);
341 
342 template <int NPT, typename S, typename R, typename I>
343 __global__ void OuterReductionKernel(R, const S, I, I, typename S::CoeffReturnType*);
344 #endif
345 
346 }  // end namespace internal
347 
348 
349 template <typename Op, typename Dims, typename XprType,  template <class> class MakePointer_>
350 class TensorReductionOp : public TensorBase<TensorReductionOp<Op, Dims, XprType, MakePointer_>, ReadOnlyAccessors> {
351   public:
352     typedef typename Eigen::internal::traits<TensorReductionOp>::Scalar Scalar;
353     typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
354     typedef typename internal::remove_const<typename XprType::CoeffReturnType>::type CoeffReturnType;
355     typedef typename Eigen::internal::nested<TensorReductionOp>::type Nested;
356     typedef typename Eigen::internal::traits<TensorReductionOp>::StorageKind StorageKind;
357     typedef typename Eigen::internal::traits<TensorReductionOp>::Index Index;
358 
359     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
360     TensorReductionOp(const XprType& expr, const Dims& dims) : m_expr(expr), m_dims(dims)
361     { }
362     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
363     TensorReductionOp(const XprType& expr, const Dims& dims, const Op& reducer) : m_expr(expr), m_dims(dims), m_reducer(reducer)
364     { }
365 
366     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
367     const XprType& expression() const { return m_expr; }
368     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
369     const Dims& dims() const { return m_dims; }
370     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
371     const Op& reducer() const { return m_reducer; }
372 
373   protected:
374     typename XprType::Nested m_expr;
375     const Dims m_dims;
376     const Op m_reducer;
377 };
378 
379 
380 // Eval as rvalue
381 template<typename Op, typename Dims, typename ArgType, template <class> class MakePointer_, typename Device>
382 struct TensorEvaluator<const TensorReductionOp<Op, Dims, ArgType, MakePointer_>, Device>
383 {
384   typedef TensorReductionOp<Op, Dims, ArgType, MakePointer_> XprType;
385   typedef typename XprType::Index Index;
386   typedef ArgType ChildType;
387   typedef typename TensorEvaluator<ArgType, Device>::Dimensions InputDimensions;
388   static const int NumInputDims = internal::array_size<InputDimensions>::value;
389   static const int NumReducedDims = internal::array_size<Dims>::value;
390   static const int NumOutputDims = NumInputDims - NumReducedDims;
391   typedef typename internal::conditional<NumOutputDims==0, Sizes<>, DSizes<Index, NumOutputDims> >::type Dimensions;
392   typedef typename XprType::Scalar Scalar;
393   typedef TensorEvaluator<const TensorReductionOp<Op, Dims, ArgType, MakePointer_>, Device> Self;
394   static const bool InputPacketAccess = TensorEvaluator<ArgType, Device>::PacketAccess;
395   typedef typename internal::remove_const<typename XprType::CoeffReturnType>::type CoeffReturnType;
396   typedef typename PacketType<CoeffReturnType, Device>::type PacketReturnType;
397   static const int PacketSize = internal::unpacket_traits<PacketReturnType>::size;
398 
399   enum {
400     IsAligned = false,
401     PacketAccess = Self::InputPacketAccess && Op::PacketAccess,
402     Layout = TensorEvaluator<ArgType, Device>::Layout,
403     CoordAccess = false,  // to be implemented
404     RawAccess = false
405   };
406 
407   static const bool ReducingInnerMostDims = internal::are_inner_most_dims<Dims, NumInputDims, Layout>::value;
408   static const bool PreservingInnerMostDims = internal::preserve_inner_most_dims<Dims, NumInputDims, Layout>::value;
409   static const bool RunningFullReduction = (NumOutputDims==0);
410 
411   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorEvaluator(const XprType& op, const Device& device)
412       : m_impl(op.expression(), device), m_reducer(op.reducer()), m_result(NULL), m_device(device), m_xpr_dims(op.dims())
413   {
414     EIGEN_STATIC_ASSERT((NumInputDims >= NumReducedDims), YOU_MADE_A_PROGRAMMING_MISTAKE);
415     EIGEN_STATIC_ASSERT((!ReducingInnerMostDims | !PreservingInnerMostDims | (NumReducedDims == NumInputDims)),
416                         YOU_MADE_A_PROGRAMMING_MISTAKE);
417 
418     // Build the bitmap indicating if an input dimension is reduced or not.
419     for (int i = 0; i < NumInputDims; ++i) {
420       m_reduced[i] = false;
421     }
422     for (int i = 0; i < NumReducedDims; ++i) {
423       eigen_assert(op.dims()[i] >= 0);
424       eigen_assert(op.dims()[i] < NumInputDims);
425       m_reduced[op.dims()[i]] = true;
426     }
427 
428     const typename TensorEvaluator<ArgType, Device>::Dimensions& input_dims = m_impl.dimensions();
429     internal::DimInitializer<Dimensions>::run(input_dims, m_reduced, &m_dimensions, &m_reducedDims);
430 
431     // Precompute output strides.
432     if (NumOutputDims > 0) {
433       if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
434         m_outputStrides[0] = 1;
435         for (int i = 1; i < NumOutputDims; ++i) {
436           m_outputStrides[i] = m_outputStrides[i - 1] * m_dimensions[i - 1];
437         }
438       } else {
439         m_outputStrides.back() = 1;
440         for (int i = NumOutputDims - 2; i >= 0; --i) {
441           m_outputStrides[i] = m_outputStrides[i + 1] * m_dimensions[i + 1];
442         }
443       }
444     }
445 
446     // Precompute input strides.
447     if (NumInputDims > 0) {
448       array<Index, NumInputDims> input_strides;
449       if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
450         input_strides[0] = 1;
451         for (int i = 1; i < NumInputDims; ++i) {
452           input_strides[i] = input_strides[i-1] * input_dims[i-1];
453         }
454       } else {
455         input_strides.back() = 1;
456         for (int i = NumInputDims - 2; i >= 0; --i) {
457           input_strides[i] = input_strides[i + 1] * input_dims[i + 1];
458         }
459       }
460 
461       int outputIndex = 0;
462       int reduceIndex = 0;
463       for (int i = 0; i < NumInputDims; ++i) {
464         if (m_reduced[i]) {
465           m_reducedStrides[reduceIndex] = input_strides[i];
466           ++reduceIndex;
467         } else {
468           m_preservedStrides[outputIndex] = input_strides[i];
469           ++outputIndex;
470         }
471       }
472     }
473 
474     // Special case for full reductions
475     if (NumOutputDims == 0) {
476       m_preservedStrides[0] = internal::array_prod(input_dims);
477     }
478   }
479 
480   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { return m_dimensions; }
481 
482   EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool evalSubExprsIfNeeded(typename MakePointer_<CoeffReturnType>::Type data) {
483     m_impl.evalSubExprsIfNeeded(NULL);
484 
485     // Use the FullReducer if possible.
486     if ((RunningFullReduction && RunningOnSycl) ||(RunningFullReduction &&
487         internal::FullReducer<Self, Op, Device>::HasOptimizedImplementation &&
488         ((RunningOnGPU && (m_device.majorDeviceVersion() >= 3)) ||
489          !RunningOnGPU))) {
490       bool need_assign = false;
491       if (!data) {
492         m_result = static_cast<CoeffReturnType*>(m_device.allocate(sizeof(CoeffReturnType)));
493         data = m_result;
494         need_assign = true;
495       }
496       Op reducer(m_reducer);
497       internal::FullReducer<Self, Op, Device>::run(*this, reducer, m_device, data);
498       return need_assign;
499     }
500     else if(RunningOnSycl){
501       const Index num_values_to_reduce = internal::array_prod(m_reducedDims);
502       const Index num_coeffs_to_preserve = internal::array_prod(m_dimensions);
503       if (!data) {
504         data = static_cast<CoeffReturnType*>(m_device.allocate(sizeof(CoeffReturnType) * num_coeffs_to_preserve));
505         m_result = data;
506       }
507       Op reducer(m_reducer);
508       internal::InnerReducer<Self, Op, Device>::run(*this, reducer, m_device, data, num_values_to_reduce, num_coeffs_to_preserve);
509       return (m_result != NULL);
510     }
511 
512     // Attempt to use an optimized reduction.
513     else if (RunningOnGPU && (m_device.majorDeviceVersion() >= 3)) {
514       bool reducing_inner_dims = true;
515       for (int i = 0; i < NumReducedDims; ++i) {
516         if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
517           reducing_inner_dims &= m_reduced[i];
518         } else {
519           reducing_inner_dims &= m_reduced[NumInputDims - 1 - i];
520         }
521       }
522       if (internal::InnerReducer<Self, Op, Device>::HasOptimizedImplementation &&
523           (reducing_inner_dims || ReducingInnerMostDims)) {
524         const Index num_values_to_reduce = internal::array_prod(m_reducedDims);
525         const Index num_coeffs_to_preserve = internal::array_prod(m_dimensions);
526         if (!data) {
527           if (num_coeffs_to_preserve < 1024 && num_values_to_reduce > num_coeffs_to_preserve && num_values_to_reduce > 128) {
528             data = static_cast<CoeffReturnType*>(m_device.allocate(sizeof(CoeffReturnType) * num_coeffs_to_preserve));
529             m_result = data;
530           }
531           else {
532             return true;
533           }
534         }
535         Op reducer(m_reducer);
536         if (internal::InnerReducer<Self, Op, Device>::run(*this, reducer, m_device, data, num_values_to_reduce, num_coeffs_to_preserve)) {
537           if (m_result) {
538             m_device.deallocate(m_result);
539             m_result = NULL;
540           }
541           return true;
542         } else {
543           return (m_result != NULL);
544         }
545       }
546 
547       bool preserving_inner_dims = true;
548       for (int i = 0; i < NumReducedDims; ++i) {
549         if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
550           preserving_inner_dims &= m_reduced[NumInputDims - 1 - i];
551         } else {
552           preserving_inner_dims &= m_reduced[i];
553         }
554       }
555       if (internal::OuterReducer<Self, Op, Device>::HasOptimizedImplementation &&
556           preserving_inner_dims) {
557         const Index num_values_to_reduce = internal::array_prod(m_reducedDims);
558         const Index num_coeffs_to_preserve = internal::array_prod(m_dimensions);
559         if (!data) {
560           if (num_coeffs_to_preserve < 1024 && num_values_to_reduce > num_coeffs_to_preserve && num_values_to_reduce > 32) {
561             data = static_cast<CoeffReturnType*>(m_device.allocate(sizeof(CoeffReturnType) * num_coeffs_to_preserve));
562             m_result = data;
563           }
564           else {
565             return true;
566           }
567         }
568         Op reducer(m_reducer);
569         if (internal::OuterReducer<Self, Op, Device>::run(*this, reducer, m_device, data, num_values_to_reduce, num_coeffs_to_preserve)) {
570           if (m_result) {
571             m_device.deallocate(m_result);
572             m_result = NULL;
573           }
574           return true;
575         } else {
576           return (m_result != NULL);
577         }
578       }
579     }
580     return true;
581   }
582 
583   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void cleanup() {
584     m_impl.cleanup();
585     if (m_result) {
586       m_device.deallocate(m_result);
587       m_result = NULL;
588     }
589   }
590 
591   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE CoeffReturnType coeff(Index index) const
592   {
593     if ((RunningOnSycl || RunningFullReduction || RunningOnGPU) && m_result) {
594       return *(m_result + index);
595     }
596     Op reducer(m_reducer);
597     if (ReducingInnerMostDims || RunningFullReduction) {
598       const Index num_values_to_reduce =
599         (static_cast<int>(Layout) == static_cast<int>(ColMajor)) ? m_preservedStrides[0] : m_preservedStrides[NumPreservedStrides - 1];
600       return internal::InnerMostDimReducer<Self, Op>::reduce(*this, firstInput(index),
601                                                              num_values_to_reduce, reducer);
602     } else {
603       typename Self::CoeffReturnType accum = reducer.initialize();
604       internal::GenericDimReducer<NumReducedDims-1, Self, Op>::reduce(*this, firstInput(index), reducer, &accum);
605       return reducer.finalize(accum);
606     }
607   }
608 
609   // TODO(bsteiner): provide a more efficient implementation.
610   template<int LoadMode>
611   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PacketReturnType packet(Index index) const
612   {
613     EIGEN_STATIC_ASSERT((PacketSize > 1), YOU_MADE_A_PROGRAMMING_MISTAKE)
614     eigen_assert(index + PacketSize - 1 < Index(internal::array_prod(dimensions())));
615 
616     if (RunningOnGPU && m_result) {
617       return internal::pload<PacketReturnType>(m_result + index);
618     }
619 
620     EIGEN_ALIGN_MAX typename internal::remove_const<CoeffReturnType>::type values[PacketSize];
621     if (ReducingInnerMostDims) {
622       const Index num_values_to_reduce =
623         (static_cast<int>(Layout) == static_cast<int>(ColMajor)) ? m_preservedStrides[0] : m_preservedStrides[NumPreservedStrides - 1];
624       const Index firstIndex = firstInput(index);
625       for (Index i = 0; i < PacketSize; ++i) {
626         Op reducer(m_reducer);
627         values[i] = internal::InnerMostDimReducer<Self, Op>::reduce(*this, firstIndex + i * num_values_to_reduce,
628                                                                     num_values_to_reduce, reducer);
629       }
630     } else if (PreservingInnerMostDims) {
631       const Index firstIndex = firstInput(index);
632       const int innermost_dim = (static_cast<int>(Layout) == static_cast<int>(ColMajor)) ? 0 : NumOutputDims - 1;
633       // TBD: extend this the the n innermost dimensions that we preserve.
634       if (((firstIndex % m_dimensions[innermost_dim]) + PacketSize - 1) < m_dimensions[innermost_dim]) {
635         Op reducer(m_reducer);
636         typename Self::PacketReturnType accum = reducer.template initializePacket<typename Self::PacketReturnType>();
637         internal::InnerMostDimPreserver<NumReducedDims-1, Self, Op>::reduce(*this, firstIndex, reducer, &accum);
638         return reducer.finalizePacket(accum);
639       } else {
640         for (int i = 0; i < PacketSize; ++i) {
641           values[i] = coeff(index + i);
642         }
643       }
644     } else {
645       for (int i = 0; i < PacketSize; ++i) {
646         values[i] = coeff(index + i);
647       }
648     }
649     PacketReturnType rslt = internal::pload<PacketReturnType>(values);
650     return rslt;
651   }
652 
653   // Must be called after evalSubExprsIfNeeded().
654   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorOpCost costPerCoeff(bool vectorized) const {
655     if (RunningFullReduction && m_result) {
656       return TensorOpCost(sizeof(CoeffReturnType), 0, 0, vectorized, PacketSize);
657     } else {
658       const Index num_values_to_reduce = internal::array_prod(m_reducedDims);
659       const double compute_cost = num_values_to_reduce * internal::functor_traits<Op>::Cost;
660       return m_impl.costPerCoeff(vectorized) * num_values_to_reduce +
661           TensorOpCost(0, 0, compute_cost, vectorized, PacketSize);
662     }
663   }
664 
665   EIGEN_DEVICE_FUNC typename MakePointer_<Scalar>::Type data() const { return m_result; }
666   /// required by sycl in order to extract the accessor
667   const TensorEvaluator<ArgType, Device>& impl() const { return m_impl; }
668   /// added for sycl in order to construct the buffer from the sycl device
669   const Device& device() const{return m_device;}
670   /// added for sycl in order to re-construct the reduction eval on the device for the sub-kernel
671   const Dims& xprDims() const {return m_xpr_dims;}
672 
673 
674   private:
675   template <int, typename, typename> friend struct internal::GenericDimReducer;
676   template <typename, typename, bool> friend struct internal::InnerMostDimReducer;
677   template <int, typename, typename, bool> friend struct internal::InnerMostDimPreserver;
678   template <typename S, typename O, typename D, bool V> friend struct internal::FullReducer;
679 #ifdef EIGEN_USE_THREADS
680   template <typename S, typename O, bool V> friend struct internal::FullReducerShard;
681 #endif
682 #if defined(EIGEN_USE_GPU) && defined(__CUDACC__)
683   template <int B, int N, typename S, typename R, typename I> friend void internal::FullReductionKernel(R, const S, I, typename S::CoeffReturnType*, unsigned int*);
684 #ifdef EIGEN_HAS_CUDA_FP16
685   template <typename S, typename R, typename I> friend void internal::ReductionInitFullReduxKernelHalfFloat(R, const S, I, half2*);
686   template <int B, int N, typename S, typename R, typename I> friend void internal::FullReductionKernelHalfFloat(R, const S, I, half*, half2*);
687   template <int NPT, typename S, typename R, typename I> friend void internal::InnerReductionKernelHalfFloat(R, const S, I, I, half*);
688 #endif
689   template <int NPT, typename S, typename R, typename I> friend void internal::InnerReductionKernel(R, const S, I, I, typename S::CoeffReturnType*);
690 
691   template <int NPT, typename S, typename R, typename I> friend void internal::OuterReductionKernel(R, const S, I, I, typename S::CoeffReturnType*);
692 #endif
693 
694   template <typename S, typename O, typename D> friend struct internal::InnerReducer;
695 
696   // Returns the Index in the input tensor of the first value that needs to be
697   // used to compute the reduction at output index "index".
698   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index firstInput(Index index) const {
699     if (ReducingInnerMostDims) {
700       if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
701         return index * m_preservedStrides[0];
702       } else {
703         return index * m_preservedStrides[NumPreservedStrides - 1];
704       }
705     }
706     // TBD: optimize the case where we preserve the innermost dimensions.
707     Index startInput = 0;
708     if (static_cast<int>(Layout) == static_cast<int>(ColMajor)) {
709       for (int i = NumOutputDims - 1; i > 0; --i) {
710         // This is index_i in the output tensor.
711         const Index idx = index / m_outputStrides[i];
712         startInput += idx * m_preservedStrides[i];
713         index -= idx * m_outputStrides[i];
714       }
715       if (PreservingInnerMostDims) {
716         eigen_assert(m_preservedStrides[0] == 1);
717         startInput += index;
718       } else {
719         startInput += index * m_preservedStrides[0];
720       }
721     } else {
722       for (int i = 0; i < NumOutputDims - 1; ++i) {
723         // This is index_i in the output tensor.
724         const Index idx = index / m_outputStrides[i];
725         startInput += idx * m_preservedStrides[i];
726         index -= idx * m_outputStrides[i];
727       }
728       if (PreservingInnerMostDims) {
729         eigen_assert(m_preservedStrides[NumPreservedStrides - 1] == 1);
730         startInput += index;
731       } else {
732         startInput += index * m_preservedStrides[NumPreservedStrides - 1];
733       }
734     }
735     return startInput;
736   }
737 
738   // Bitmap indicating if an input dimension is reduced or not.
739   array<bool, NumInputDims> m_reduced;
740   // Dimensions of the output of the operation.
741   Dimensions m_dimensions;
742   // Precomputed strides for the output tensor.
743   array<Index, NumOutputDims> m_outputStrides;
744   // Subset of strides of the input tensor for the non-reduced dimensions.
745   // Indexed by output dimensions.
746   static const int NumPreservedStrides = max_n_1<NumOutputDims>::size;
747   array<Index, NumPreservedStrides> m_preservedStrides;
748 
749   // Subset of strides of the input tensor for the reduced dimensions.
750   // Indexed by reduced dimensions.
751   array<Index, NumReducedDims> m_reducedStrides;
752   // Size of the input dimensions that are reduced.
753   // Indexed by reduced dimensions.
754   array<Index, NumReducedDims> m_reducedDims;
755 
756   // Evaluator for the input expression.
757   TensorEvaluator<ArgType, Device> m_impl;
758 
759   // Operation to apply for computing the reduction.
760   Op m_reducer;
761 
762   // For full reductions
763 #if defined(EIGEN_USE_GPU) && defined(__CUDACC__)
764   static const bool RunningOnGPU = internal::is_same<Device, Eigen::GpuDevice>::value;
765   static const bool RunningOnSycl = false;
766 #elif defined(EIGEN_USE_SYCL)
767 static const bool RunningOnSycl = internal::is_same<typename internal::remove_all<Device>::type, Eigen::SyclDevice>::value;
768 static const bool RunningOnGPU = false;
769 #else
770   static const bool RunningOnGPU = false;
771   static const bool RunningOnSycl = false;
772 #endif
773   typename MakePointer_<CoeffReturnType>::Type m_result;
774 
775   const Device& m_device;
776   const Dims& m_xpr_dims;
777 };
778 
779 } // end namespace Eigen
780 
781 #endif // EIGEN_CXX11_TENSOR_TENSOR_REDUCTION_H
782