diff --git a/sycl/include/sycl/reduction.hpp b/sycl/include/sycl/reduction.hpp
index 47dc36c3e0bc6..2144f547cbf61 100644
--- a/sycl/include/sycl/reduction.hpp
+++ b/sycl/include/sycl/reduction.hpp
@@ -147,6 +147,9 @@ __SYCL_EXPORT size_t reduComputeWGSize(size_t NWorkItems, size_t MaxWGSize,
 __SYCL_EXPORT size_t reduGetPreferredWGSize(std::shared_ptr<queue_impl> &Queue,
                                             size_t LocalMemBytesPerWorkItem);
 
+template <typename T, class BinaryOperation, bool IsOptional>
+class ReducerElement;
+
 /// Helper class for accessing reducer-defined types in CRTP
 /// May prove to be useful for other things later
 template <typename Reducer> struct ReducerTraits;
@@ -160,6 +163,7 @@ struct ReducerTraits<reducer<T, BinaryOperation, Dims, Extent,
   static constexpr int dims = Dims;
   static constexpr size_t extent = Extent;
   static constexpr bool has_identity = IdentityContainerT::has_identity;
+  using element_type = ReducerElement<T, BinaryOperation, !has_identity>;
 };
 
 /// Helper class for accessing internal reducer member functions.
@@ -275,10 +279,16 @@ template <class Reducer> class combiner {
   void atomic_combine_impl(T *ReduVarPtr, AtomicFunctor Functor) const {
     auto reducer = static_cast<const Reducer *>(this);
     for (size_t E = 0; E < Extent; ++E) {
+      const auto &ReducerElem = getReducerAccess(*reducer).getElement(E);
+
+      // If the reducer element doesn't have a value we can skip the combine.
+      if (!ReducerElem)
+        continue;
+
       auto AtomicRef = sycl::atomic_ref<T, memory_order::relaxed,
                                         getMemoryScope<Space>(), Space>(
           address_space_cast<Space, access::decorated::no>(ReduVarPtr)[E]);
-      Functor(std::move(AtomicRef), getReducerAccess(*reducer).getElement(E));
+      Functor(std::move(AtomicRef), *ReducerElem);
     }
   }
 
@@ -361,13 +371,6 @@ template <class Reducer> class combiner {
   }
 };
 
-template <typename T, class BinaryOperation, int Dims> class reducer_common {
-public:
-  using value_type = T;
-  using binary_operation = BinaryOperation;
-  static constexpr int dimensions = Dims;
-};
-
 /// Templated class for common functionality of all reduction implementation
 /// classes.
 template <typename T, class BinaryOperation, bool ExplicitIdentity,
@@ -422,6 +425,85 @@ class ReductionIdentityContainer<
   static constexpr bool has_identity = false;
 };
 
+// Class representing a single reducer element. In cases where the identity is
+// unknown this is allowed to not contain a value.
+template <typename T, class BinaryOperation, bool IsOptional>
+class ReducerElement {
+  using value_type = std::conditional_t<IsOptional, std::optional<T>, T>;
+
+  template <bool ExplicitIdentity>
+  constexpr value_type GetInitValue(
+      const ReductionIdentityContainer<T, BinaryOperation, ExplicitIdentity>
+          &IdentityContainer) {
+    constexpr bool ContainerHasIdentity =
+        ReductionIdentityContainer<T, BinaryOperation,
+                                   ExplicitIdentity>::has_identity;
+    static_assert(IsOptional || ContainerHasIdentity);
+    if constexpr (!ContainerHasIdentity)
+      return std::nullopt;
+    else
+      return IdentityContainer.getIdentity();
+  }
+
+public:
+  ReducerElement() = default;
+  ReducerElement(T Value) : MValue{Value} {}
+
+  template <bool ExplicitIdentity>
+  ReducerElement(
+      const ReductionIdentityContainer<T, BinaryOperation, ExplicitIdentity>
+          &IdentityContainer)
+      : MValue(GetInitValue(IdentityContainer)) {}
+
+  ReducerElement &combine(BinaryOperation BinOp, const T &OtherValue) {
+    if constexpr (IsOptional)
+      MValue = MValue ? BinOp(*MValue, OtherValue) : OtherValue;
+    else
+      MValue = BinOp(MValue, OtherValue);
+    return *this;
+  }
+
+  ReducerElement &combine(BinaryOperation BinOp, const ReducerElement &Other) {
+    if constexpr (IsOptional) {
+      if (Other.MValue)
+        return combine(BinOp, *Other.MValue);
+      // If the other value doesn't have a value it is a no-op.
+      return *this;
+    } else {
+      return combine(BinOp, Other.MValue);
+    }
+  }
+
+  constexpr T &operator*() noexcept {
+    if constexpr (IsOptional)
+      return *MValue;
+    else
+      return MValue;
+  }
+  constexpr const T &operator*() const noexcept {
+    if constexpr (IsOptional)
+      return *MValue;
+    else
+      return MValue;
+  }
+
+  constexpr explicit operator bool() const {
+    if constexpr (IsOptional)
+      return MValue.has_value();
+    return true;
+  }
+
+private:
+  value_type MValue;
+};
+
+template <typename T, class BinaryOperation, int Dims> class reducer_common {
+public:
+  using value_type = T;
+  using binary_operation = BinaryOperation;
+  static constexpr int dimensions = Dims;
+};
+
 // Token class to help with the in-place construction of reducers.
 template <class BinaryOperation, typename IdentityContainerT>
 struct ReducerToken {
@@ -450,21 +532,13 @@ class reducer<
                       !detail::IsKnownIdentityOp<T, BinaryOperation>::value>>>,
       public detail::reducer_common<T, BinaryOperation, Dims> {
   static constexpr bool has_identity = IdentityContainerT::has_identity;
-  using internal_value_type =
-      std::conditional_t<has_identity, T, std::optional<T>>;
-
-  constexpr internal_value_type
-  GetInitialValue(const IdentityContainerT &IdentityContainer) {
-    if constexpr (has_identity)
-      return IdentityContainer.getIdentity();
-    else
-      return std::nullopt;
-  }
+  using element_type =
+      detail::ReducerElement<T, BinaryOperation, !has_identity>;
 
 public:
   reducer(const IdentityContainerT &IdentityContainer, BinaryOperation BOp)
-      : MValue(GetInitialValue(IdentityContainer)),
-        MIdentity(IdentityContainer), MBinaryOp(BOp) {}
+      : MValue(IdentityContainer), MIdentity(IdentityContainer),
+        MBinaryOp(BOp) {}
   reducer(
       const detail::ReducerToken<BinaryOperation, IdentityContainerT> &Token)
       : reducer(Token.IdentityContainer, Token.BOp) {}
@@ -475,10 +549,7 @@ class reducer<
   reducer &operator=(reducer &&) = delete;
 
   reducer &combine(const T &Partial) {
-    if constexpr (has_identity)
-      MValue = MBinaryOp(MValue, Partial);
-    else
-      MValue = MValue ? MBinaryOp(*MValue, Partial) : Partial;
+    MValue.combine(MBinaryOp, Partial);
     return *this;
   }
 
@@ -491,20 +562,10 @@ class reducer<
 private:
   template <typename ReducerT> friend class detail::ReducerAccess;
 
-  T &getElement(size_t) {
-    if constexpr (has_identity)
-      return MValue;
-    else
-      return *MValue;
-  }
-  const T &getElement(size_t) const {
-    if constexpr (has_identity)
-      return MValue;
-    else
-      return *MValue;
-  }
+  element_type &getElement(size_t) { return MValue; }
+  const element_type &getElement(size_t) const { return MValue; }
 
-  internal_value_type MValue;
+  detail::ReducerElement<T, BinaryOperation, !has_identity> MValue;
   const IdentityContainerT MIdentity;
   BinaryOperation MBinaryOp;
 };
@@ -526,6 +587,10 @@ class reducer<
                       Dims == 0 && Extent == 1 && View == false &&
                       detail::IsKnownIdentityOp<T, BinaryOperation>::value>>>,
       public detail::reducer_common<T, BinaryOperation, Dims> {
+  static constexpr bool has_identity = IdentityContainerT::has_identity;
+  using element_type =
+      detail::ReducerElement<T, BinaryOperation, !has_identity>;
+
 public:
   reducer() : MValue(getIdentity()) {}
   reducer(const IdentityContainerT & /* Identity */, BinaryOperation)
@@ -541,7 +606,7 @@ class reducer<
 
   reducer &combine(const T &Partial) {
     BinaryOperation BOp;
-    MValue = BOp(MValue, Partial);
+    MValue.combine(BOp, Partial);
     return *this;
   }
 
@@ -554,9 +619,9 @@ class reducer<
     return detail::known_identity_impl<BinaryOperation, T>::value;
   }
 
-  T &getElement(size_t) { return MValue; }
-  const T &getElement(size_t) const { return MValue; }
-  T MValue;
+  element_type &getElement(size_t) { return MValue; }
+  const element_type &getElement(size_t) const { return MValue; }
+  detail::ReducerElement<T, BinaryOperation, !has_identity> MValue;
 };
 
 /// Component of 'reducer' class for array reductions, representing a single
@@ -570,11 +635,11 @@ class reducer<T, BinaryOperation, Dims, Extent, IdentityContainerT, View,
                   std::enable_if_t<Dims == 0 && View == true>>>,
       public detail::reducer_common<T, BinaryOperation, Dims> {
   static constexpr bool has_identity = IdentityContainerT::has_identity;
-  using internal_value_type =
-      std::conditional_t<has_identity, T, std::optional<T>>;
+  using element_type =
+      detail::ReducerElement<T, BinaryOperation, !has_identity>;
 
 public:
-  reducer(internal_value_type &Ref, BinaryOperation BOp)
+  reducer(element_type &Ref, BinaryOperation BOp)
       : MElement(Ref), MBinaryOp(BOp) {}
   reducer(
       const detail::ReducerToken<BinaryOperation, IdentityContainerT> &Token)
@@ -586,17 +651,17 @@ class reducer<T, BinaryOperation, Dims, Extent, IdentityContainerT, View,
   reducer &operator=(reducer &&) = delete;
 
   reducer &combine(const T &Partial) {
-    if constexpr (has_identity)
-      MElement = MBinaryOp(MElement, Partial);
-    else
-      MElement = MElement ? MBinaryOp(*MElement, Partial) : Partial;
+    MElement.combine(MBinaryOp, Partial);
     return *this;
   }
 
 private:
   template <typename ReducerT> friend class detail::ReducerAccess;
 
-  internal_value_type &MElement;
+  element_type &getElement(size_t) { return MElement; }
+  const element_type &getElement(size_t) const { return MElement; }
+
+  element_type &MElement;
   BinaryOperation MBinaryOp;
 };
 
@@ -615,21 +680,13 @@ class reducer<
                       !detail::IsKnownIdentityOp<T, BinaryOperation>::value>>>,
       public detail::reducer_common<T, BinaryOperation, Dims> {
   static constexpr bool has_identity = IdentityContainerT::has_identity;
-  using internal_value_type =
-      std::conditional_t<has_identity, T, std::optional<T>>;
-
-  constexpr internal_value_type
-  GetInitialValue(const IdentityContainerT &IdentityContainer) {
-    if constexpr (has_identity)
-      return IdentityContainer.getIdentity();
-    else
-      return std::nullopt;
-  }
+  using element_type =
+      detail::ReducerElement<T, BinaryOperation, !has_identity>;
 
 public:
   reducer(const IdentityContainerT &IdentityContainer, BinaryOperation BOp)
-      : MValue(GetInitialValue(IdentityContainer)),
-        MIdentity(IdentityContainer), MBinaryOp(BOp) {}
+      : MValue(IdentityContainer), MIdentity(IdentityContainer),
+        MBinaryOp(BOp) {}
   reducer(
       const detail::ReducerToken<BinaryOperation, IdentityContainerT> &Token)
       : reducer(Token.IdentityContainer, Token.BOp) {}
@@ -653,20 +710,10 @@ class reducer<
 private:
   template <typename ReducerT> friend class detail::ReducerAccess;
 
-  T &getElement(size_t E) {
-    if constexpr (has_identity)
-      return MValue[E];
-    else
-      return *(MValue[E]);
-  }
-  const T &getElement(size_t E) const {
-    if constexpr (has_identity)
-      return MValue[E];
-    else
-      return *(MValue[E]);
-  }
+  element_type &getElement(size_t E) { return MValue[E]; }
+  const element_type &getElement(size_t E) const { return MValue[E]; }
 
-  marray<internal_value_type, Extent> MValue;
+  marray<element_type, Extent> MValue;
   const IdentityContainerT MIdentity;
   BinaryOperation MBinaryOp;
 };
@@ -685,6 +732,10 @@ class reducer<
                       Dims == 1 && View == false &&
                       detail::IsKnownIdentityOp<T, BinaryOperation>::value>>>,
       public detail::reducer_common<T, BinaryOperation, Dims> {
+  static constexpr bool has_identity = IdentityContainerT::has_identity;
+  using element_type =
+      detail::ReducerElement<T, BinaryOperation, !has_identity>;
+
 public:
   reducer() : MValue(getIdentity()) {}
   reducer(const IdentityContainerT & /* Identity */, BinaryOperation)
@@ -714,10 +765,10 @@ class reducer<
     return detail::known_identity_impl<BinaryOperation, T>::value;
   }
 
-  T &getElement(size_t E) { return MValue[E]; }
-  const T &getElement(size_t E) const { return MValue[E]; }
+  element_type &getElement(size_t E) { return MValue[E]; }
+  const element_type &getElement(size_t E) const { return MValue[E]; }
 
-  marray<T, Extent> MValue;
+  marray<element_type, Extent> MValue;
 };
 
 namespace detail {
@@ -739,10 +790,8 @@ template <class T> struct get_red_t<T *> {
   using type = T;
 };
 
-template <class T, int AccessorDims, access::mode Mode,
-          access::placeholder IsPH, typename PropList>
-struct get_red_t<
-    accessor<T, AccessorDims, Mode, access::target::device, IsPH, PropList>> {
+template <class T, int Dims, typename AllocatorT>
+struct get_red_t<buffer<T, Dims, AllocatorT>> {
   using type = T;
 };
 
@@ -797,6 +846,8 @@ class reduction_impl_algo {
       detail::ReducerToken<BinaryOperation, identity_container_type>;
   using reducer_type =
       reducer<T, BinaryOperation, Dims, Extent, identity_container_type>;
+  using reducer_element_type =
+      typename ReducerTraits<reducer_type>::element_type;
   using result_type = T;
   using binary_operation = BinaryOperation;
 
@@ -846,9 +897,10 @@ class reduction_impl_algo {
     // If there is only one WG we can avoid creation of temporary buffer with
     // partial sums and write directly into user's reduction variable.
     if constexpr (IsOneWG) {
-      return MRedOut;
+      return getUserRedVarAccess(CGH);
     } else {
-      MOutBufPtr = std::make_shared<buffer<T, 1>>(range<1>(Size));
+      MOutBufPtr =
+          std::make_shared<buffer<reducer_element_type, 1>>(range<1>(Size));
       CGH.addReduction(MOutBufPtr);
       return accessor{*MOutBufPtr, CGH};
     }
@@ -867,17 +919,24 @@ class reduction_impl_algo {
   /// Otherwise, a new buffer is created and accessor to that buffer is
   /// returned.
   auto getWriteAccForPartialReds(size_t Size, handler &CGH) {
-    if constexpr (!is_usm) {
+    static_assert(!has_identity || sizeof(reducer_element_type) == sizeof(T),
+                  "Unexpected size of reducer element type.");
+
+    // We can only use the output memory directly if it is not USM and we have
+    // and identity, i.e. it has a thin element wrapper.
+    if constexpr (!is_usm && has_identity) {
       if (Size == 1) {
-        CGH.associateWithHandler(&MRedOut, access::target::device);
-        return MRedOut;
+        auto ReinterpretRedOut =
+            MRedOut.template reinterpret<reducer_element_type>();
+        return accessor{ReinterpretRedOut, CGH};
       }
     }
 
     // Create a new output buffer and return an accessor to it.
     //
     // Array reductions are performed element-wise to avoid stack growth.
-    MOutBufPtr = std::make_shared<buffer<T, 1>>(range<1>(Size));
+    MOutBufPtr =
+        std::make_shared<buffer<reducer_element_type, 1>>(range<1>(Size));
     CGH.addReduction(MOutBufPtr);
     return accessor{*MOutBufPtr, CGH};
   }
@@ -934,8 +993,8 @@ class reduction_impl_algo {
             }
           });
         } else {
-          associateWithHandler(CopyHandler, &Out, access::target::device);
-          CopyHandler.copy(Mem, Out);
+          accessor OutAcc{Out, CGH};
+          CopyHandler.copy(Mem, OutAcc);
         }
       });
     };
@@ -947,7 +1006,7 @@ class reduction_impl_algo {
       if (initializeToIdentity())
         DoIt(MRedOut);
       else
-        Func(MRedOut);
+        Func(accessor{MRedOut, CGH});
     }
   }
 
@@ -959,7 +1018,7 @@ class reduction_impl_algo {
     std::ignore = CGH;
     assert(!initializeToIdentity() &&
            "Initialize to identity not allowed for identity-less reductions.");
-    Func(MRedOut);
+    Func(accessor{MRedOut, CGH});
   }
 
   const identity_container_type &getIdentityContainer() {
@@ -997,7 +1056,13 @@ class reduction_impl_algo {
   BinaryOperation getBinaryOperation() const { return MBinaryOp; }
   bool initializeToIdentity() const { return InitializeToIdentity; }
 
-  RedOutVar &getUserRedVar() { return MRedOut; }
+  auto getUserRedVarAccess(handler &CGH) {
+    std::ignore = CGH;
+    if constexpr (is_usm)
+      return MRedOut;
+    else
+      return accessor{MRedOut, CGH};
+  }
 
 private:
   // Object holding the identity if available.
@@ -1005,7 +1070,7 @@ class reduction_impl_algo {
 
   // Array reduction is performed element-wise to avoid stack growth, hence
   // 1-dimensional always.
-  std::shared_ptr<buffer<T, 1>> MOutBufPtr;
+  std::shared_ptr<buffer<reducer_element_type, 1>> MOutBufPtr;
 
   BinaryOperation MBinaryOp;
   bool InitializeToIdentity;
@@ -1103,15 +1168,16 @@ void reduSaveFinalResultToUserMem(handler &CGH, Reduction &Redu) {
                 "USM reduction, not a buffer-based one.");
   size_t NElements = Reduction::num_elements;
   auto InAcc = Redu.getReadAccToPreviousPartialReds(CGH);
-  auto UserVarPtr = Redu.getUserRedVar();
+  auto UserVarPtr = Redu.getUserRedVarAccess(CGH);
   bool IsUpdateOfUserVar = !Redu.initializeToIdentity();
   auto BOp = Redu.getBinaryOperation();
   CGH.single_task<KernelName>([=] {
     for (int i = 0; i < NElements; ++i) {
+      auto Elem = InAcc[i];
       if (IsUpdateOfUserVar)
-        UserVarPtr[i] = BOp(UserVarPtr[i], InAcc.get_pointer()[i]);
+        UserVarPtr[i] = BOp(UserVarPtr[i], *Elem);
       else
-        UserVarPtr[i] = InAcc.get_pointer()[i];
+        UserVarPtr[i] = *Elem;
     }
   });
 }
@@ -1121,6 +1187,10 @@ template <typename KernelName, strategy S, class... Ts> struct MainKrn;
 template <typename KernelName, strategy S, class... Ts> struct AuxKrn;
 } // namespace reduction
 
+// Tag structs to help creating unique kernels for multi-reduction cases.
+struct KernelOneWGTag {};
+struct KernelMultipleWGTag {};
+
 /// A helper to pass undefined (sycl::detail::auto_name) names unmodified. We
 /// must do that to avoid name collisions.
 template <template <typename, reduction::strategy, typename...> class MainOrAux,
@@ -1140,6 +1210,10 @@ struct NDRangeReduction<reduction::strategy::local_atomic_and_atomic_cross_wg> {
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    static_assert(Reduction::has_identity,
+                  "Identityless reductions are not supported by the "
+                  "local_atomic_and_atomic_cross_wg strategy.");
+
     std::ignore = Queue;
     using Name = __sycl_reduction_kernel<
         reduction::MainKrn, KernelName,
@@ -1170,7 +1244,7 @@ struct NDRangeReduction<reduction::strategy::local_atomic_and_atomic_cross_wg> {
         workGroupBarrier();
         if (LID == 0) {
           for (size_t E = 0; E < NElements; ++E) {
-            getReducerAccess(Reducer).getElement(E) = GroupSum[E];
+            *getReducerAccess(Reducer).getElement(E) = GroupSum[E];
           }
           Reducer.template atomic_combine(&Out[0]);
         }
@@ -1196,9 +1270,7 @@ struct NDRangeReduction<
     size_t WGSize = NDRange.get_local_range().size();
     size_t NWorkGroups = NDRange.get_group_range().size();
 
-    auto &Out = Redu.getUserRedVar();
-    if constexpr (!Reduction::is_usm)
-      associateWithHandler(CGH, &Out, access::target::device);
+    auto Out = Redu.getUserRedVarAccess(CGH);
 
     auto &PartialSumsBuf = Redu.getTempBuffer(NWorkGroups * NElements, CGH);
     accessor PartialSums(PartialSumsBuf, CGH, sycl::read_write, sycl::no_init);
@@ -1224,7 +1296,7 @@ struct NDRangeReduction<
         // reduce_over_group is only defined for each T, not for span<T, ...>
         size_t LID = NDId.get_local_id(0);
         for (int E = 0; E < NElements; ++E) {
-          auto &RedElem = getReducerAccess(Reducer).getElement(E);
+          auto &RedElem = *getReducerAccess(Reducer).getElement(E);
           RedElem = reduce_over_group(Group, RedElem, BOp);
           if (LID == 0) {
             if (NWorkGroups == 1) {
@@ -1235,7 +1307,7 @@ struct NDRangeReduction<
               Out[E] = RedElem;
             } else {
               PartialSums[NDId.get_group_linear_id() * NElements + E] =
-                  getReducerAccess(Reducer).getElement(E);
+                  *getReducerAccess(Reducer).getElement(E);
             }
           }
         }
@@ -1384,18 +1456,33 @@ template <> struct NDRangeReduction<reduction::strategy::range_basic> {
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    using reducer_type = typename Reduction::reducer_type;
+    using element_type = typename ReducerTraits<reducer_type>::element_type;
+
+    // If reduction has an identity and is not USM, the reducer element is just
+    // a thin wrapper around the result type so the partial sum will use the
+    // output memory iff NWorkGroups == 1. Otherwise, we need to make sure the
+    // right output buffer is written in case NWorkGroups == 1.
+    constexpr bool UsePartialSumForOutput =
+        !Reduction::is_usm && Reduction::has_identity;
+
     std::ignore = Queue;
     size_t NElements = Reduction::num_elements;
     size_t WGSize = NDRange.get_local_range().size();
     size_t NWorkGroups = NDRange.get_group_range().size();
 
-    bool IsUpdateOfUserVar = !Reduction::is_usm && !Redu.initializeToIdentity();
+    bool IsUpdateOfUserVar = !Redu.initializeToIdentity();
     auto PartialSums =
         Redu.getWriteAccForPartialReds(NWorkGroups * NElements, CGH);
-    auto Out = (NWorkGroups == 1)
+    auto Out = [&]() {
+      if constexpr (UsePartialSumForOutput)
+        return (NWorkGroups == 1)
                    ? PartialSums
                    : Redu.getWriteAccForPartialReds(NElements, CGH);
-    local_accessor<typename Reduction::result_type, 1> LocalReds{WGSize, CGH};
+      else
+        return Redu.getUserRedVarAccess(CGH);
+    }();
+    local_accessor<element_type, 1> LocalReds{WGSize, CGH};
     auto NWorkGroupsFinished =
         Redu.getReadWriteAccessorToInitializedGroupsCounter(CGH);
     local_accessor<int, 1> DoReducePartialSumsInLastWG{1, CGH};
@@ -1408,25 +1495,35 @@ template <> struct NDRangeReduction<reduction::strategy::range_basic> {
 
     CGH.parallel_for<Name>(NDRange, Properties, [=](nd_item<1> NDId) {
       // Call user's functions. Reducer.MValue gets initialized there.
-      typename Reduction::reducer_type Reducer =
-          typename Reduction::reducer_type(IdentityContainer, BOp);
+      reducer_type Reducer = reducer_type(IdentityContainer, BOp);
       KernelFunc(NDId, Reducer);
 
+      auto ElementCombiner = [&](element_type &LHS, const element_type &RHS) {
+        return LHS.combine(BOp, RHS);
+      };
+
       // If there are multiple values, reduce each separately
       // This prevents local memory from scaling with elements
       size_t LID = NDId.get_local_linear_id();
       for (int E = 0; E < NElements; ++E) {
 
         doTreeReduction<WorkSizeGuarantees::Equal>(
-            WGSize, NDId, LocalReds, BOp,
+            WGSize, NDId, LocalReds, ElementCombiner,
             [&](size_t) { return getReducerAccess(Reducer).getElement(E); });
 
         if (LID == 0) {
           auto V = LocalReds[0];
-          if (NWorkGroups == 1 && IsUpdateOfUserVar)
-            V = BOp(V, Out[E]);
-          // if NWorkGroups == 1, then PartialsSum and Out point to same memory.
-          PartialSums[NDId.get_group_linear_id() * NElements + E] = V;
+
+          bool IsOneWG = NWorkGroups == 1;
+          if (IsOneWG && IsUpdateOfUserVar)
+            V.combine(BOp, Out[E]);
+
+          // if NWorkGroups == 1 && UsePartialSumForOutput, then PartialsSum
+          // and Out point to same memory.
+          if (UsePartialSumForOutput || !IsOneWG)
+            PartialSums[NDId.get_group_linear_id() * NElements + E] = V;
+          else if (V)
+            Out[E] = *V;
         }
       }
 
@@ -1446,22 +1543,17 @@ template <> struct NDRangeReduction<reduction::strategy::range_basic> {
         // TODO: Opportunity to parallelize across elements
         for (int E = 0; E < NElements; ++E) {
           doTreeReduction<WorkSizeGuarantees::None>(
-              NWorkGroups, NDId, LocalReds, BOp,
+              NWorkGroups, NDId, LocalReds, ElementCombiner,
               [&](size_t I) { return PartialSums[I * NElements + E]; });
           if (LID == 0) {
             auto V = LocalReds[0];
             if (IsUpdateOfUserVar)
-              V = BOp(V, Out[E]);
-            Out[E] = V;
+              V.combine(BOp, Out[E]);
+            Out[E] = *V;
           }
         }
       }
     });
-
-    if constexpr (Reduction::is_usm)
-      reduction::withAuxHandler(CGH, [&](handler &CopyHandler) {
-        reduSaveFinalResultToUserMem<KernelName>(CopyHandler, Redu);
-      });
   }
 };
 
@@ -1472,6 +1564,10 @@ struct NDRangeReduction<reduction::strategy::group_reduce_and_atomic_cross_wg> {
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    static_assert(Reduction::has_identity,
+                  "Identityless reductions are not supported by the "
+                  "group_reduce_and_atomic_cross_wg strategy.");
+
     std::ignore = Queue;
     using Name = __sycl_reduction_kernel<
         reduction::MainKrn, KernelName,
@@ -1486,8 +1582,8 @@ struct NDRangeReduction<reduction::strategy::group_reduce_and_atomic_cross_wg> {
 
         typename Reduction::binary_operation BOp;
         for (int E = 0; E < NElements; ++E) {
-          getReducerAccess(Reducer).getElement(E) = reduce_over_group(
-              NDIt.get_group(), getReducerAccess(Reducer).getElement(E), BOp);
+          auto &ReducerElem = getReducerAccess(Reducer).getElement(E);
+          *ReducerElem = reduce_over_group(NDIt.get_group(), *ReducerElem, BOp);
         }
         if (NDIt.get_local_linear_id() == 0)
           Reducer.atomic_combine(&Out[0]);
@@ -1504,6 +1600,9 @@ struct NDRangeReduction<
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    using reducer_type = typename Reduction::reducer_type;
+    using element_type = typename ReducerTraits<reducer_type>::element_type;
+
     std::ignore = Queue;
     using Name = __sycl_reduction_kernel<
         reduction::MainKrn, KernelName,
@@ -1514,23 +1613,27 @@ struct NDRangeReduction<
 
       // Use local memory to reduce elements in work-groups into zero-th
       // element.
-      local_accessor<typename Reduction::result_type, 1> LocalReds{WGSize, CGH};
+      local_accessor<element_type, 1> LocalReds{WGSize, CGH};
 
       CGH.parallel_for<Name>(NDRange, Properties, [=](nd_item<Dims> NDIt) {
         // Call user's functions. Reducer.MValue gets initialized there.
-        typename Reduction::reducer_type Reducer;
+        reducer_type Reducer;
         KernelFunc(NDIt, Reducer);
 
         size_t WGSize = NDIt.get_local_range().size();
         size_t LID = NDIt.get_local_linear_id();
 
+        typename Reduction::binary_operation BOp;
+        auto ElementCombiner = [&](element_type &LHS, const element_type &RHS) {
+          return LHS.combine(BOp, RHS);
+        };
+
         // If there are multiple values, reduce each separately
         // This prevents local memory from scaling with elements
         for (int E = 0; E < NElements; ++E) {
 
-          typename Reduction::binary_operation BOp;
           doTreeReduction<WorkSizeGuarantees::Equal>(
-              WGSize, NDIt, LocalReds, BOp,
+              WGSize, NDIt, LocalReds, ElementCombiner,
               [&](size_t) { return getReducerAccess(Reducer).getElement(E); });
 
           if (LID == 0)
@@ -1558,6 +1661,10 @@ struct NDRangeReduction<
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    static_assert(Reduction::has_identity,
+                  "Identityless reductions are not supported by the "
+                  "group_reduce_and_multiple_kernels strategy.");
+
     // Before running the kernels, check that device has enough local memory
     // to hold local arrays that may be required for the reduction algorithm.
     // TODO: If the work-group-size is limited by the local memory, then
@@ -1600,11 +1707,11 @@ struct NDRangeReduction<
       typename Reduction::binary_operation BOp;
       for (int E = 0; E < NElements; ++E) {
         typename Reduction::result_type PSum;
-        PSum = getReducerAccess(Reducer).getElement(E);
+        PSum = *getReducerAccess(Reducer).getElement(E);
         PSum = reduce_over_group(NDIt.get_group(), PSum, BOp);
         if (NDIt.get_local_linear_id() == 0) {
           if (IsUpdateOfUserVar)
-            PSum = BOp(Out[E], PSum);
+            PSum = BOp(*Out[E], PSum);
           Out[WGID * NElements + E] = PSum;
         }
       }
@@ -1663,13 +1770,13 @@ struct NDRangeReduction<
           for (int E = 0; E < NElements; ++E) {
             typename Reduction::result_type PSum =
                 (HasUniformWG || (GID < NWorkItems))
-                    ? In[GID * NElements + E]
+                    ? *In[GID * NElements + E]
                     : ReducerAccess<typename Reduction::reducer_type>::
                           getIdentityStatic();
             PSum = reduce_over_group(NDIt.get_group(), PSum, BOp);
             if (NDIt.get_local_linear_id() == 0) {
               if (IsUpdateOfUserVar)
-                PSum = BOp(Out[E], PSum);
+                PSum = BOp(*Out[E], PSum);
               Out[WGID * NElements + E] = PSum;
             }
           }
@@ -1692,11 +1799,13 @@ template <> struct NDRangeReduction<reduction::strategy::basic> {
   static void run(handler &CGH, std::shared_ptr<detail::queue_impl> &Queue,
                   nd_range<Dims> NDRange, PropertiesT &Properties,
                   Reduction &Redu, KernelType &KernelFunc) {
+    using element_type = typename Reduction::reducer_element_type;
+
     constexpr bool HFR = Reduction::has_fast_reduce;
-    size_t OneElemSize = HFR ? 0 : sizeof(typename Reduction::result_type);
+    size_t OneElemSize = HFR ? 0 : sizeof(element_type);
     // TODO: currently the maximal work group size is determined for the given
-    // queue/device, while it may be safer to use queries to the kernel compiled
-    // for the device.
+    // queue/device, while it may be safer to use queries to the kernel
+    // compiled for the device.
     size_t MaxWGSize = reduGetMaxWGSize(Queue, OneElemSize);
     if (NDRange.get_local_range().size() > MaxWGSize)
       throw sycl::runtime_error("The implementation handling parallel_for with"
@@ -1705,55 +1814,87 @@ template <> struct NDRangeReduction<reduction::strategy::basic> {
                                     std::to_string(MaxWGSize),
                                 PI_ERROR_INVALID_WORK_GROUP_SIZE);
 
-    size_t NElements = Reduction::num_elements;
-    size_t WGSize = NDRange.get_local_range().size();
     size_t NWorkGroups = NDRange.get_group_range().size();
-    auto Out = Redu.getWriteAccForPartialReds(NWorkGroups * NElements, CGH);
 
-    bool IsUpdateOfUserVar =
-        !Reduction::is_usm && !Redu.initializeToIdentity() && NWorkGroups == 1;
+    bool IsUpdateOfUserVar = !Redu.initializeToIdentity();
+    std::ignore = IsUpdateOfUserVar;
 
-    // Use local memory to reduce elements in work-groups into 0-th element.
-    local_accessor<typename Reduction::result_type, 1> LocalReds{WGSize, CGH};
+    // The type of the Out "accessor" differs between scenarios when there is
+    // just one WorkGroup and when there are multiple. Use this lambda to write
+    // the code just once.
+    auto First = [&](auto KernelTag) {
+      // We can deduce IsOneWG from the tag type.
+      constexpr bool IsOneWG =
+          std::is_same_v<std::remove_reference_t<decltype(KernelTag)>,
+                         KernelOneWGTag>;
 
-    auto BOp = Redu.getBinaryOperation();
-    auto IdentityContainer = Redu.getIdentityContainer();
+      constexpr size_t NElements = Reduction::num_elements;
 
-    using Name = __sycl_reduction_kernel<reduction::MainKrn, KernelName,
-                                         reduction::strategy::basic>;
+      size_t WGSize = NDRange.get_local_range().size();
 
-    CGH.parallel_for<Name>(NDRange, Properties, [=](nd_item<Dims> NDIt) {
-      // Call user's functions. Reducer.MValue gets initialized there.
-      typename Reduction::reducer_type Reducer =
-          typename Reduction::reducer_type(IdentityContainer, BOp);
-      KernelFunc(NDIt, Reducer);
+      auto Out = [&]() {
+        if constexpr (IsOneWG)
+          return Redu.getUserRedVarAccess(CGH);
+        else
+          return Redu.getWriteAccForPartialReds(NWorkGroups * NElements, CGH);
+      }();
 
-      size_t WGSize = NDIt.get_local_range().size();
-      size_t LID = NDIt.get_local_linear_id();
+      // Use local memory to reduce elements in work-groups into 0-th element.
+      local_accessor<element_type, 1> LocalReds{WGSize, CGH};
 
-      // If there are multiple values, reduce each separately
-      // This prevents local memory from scaling with elements
-      for (int E = 0; E < NElements; ++E) {
+      auto BOp = Redu.getBinaryOperation();
+      auto IdentityContainer = Redu.getIdentityContainer();
 
-        doTreeReduction<WorkSizeGuarantees::Equal>(
-            WGSize, NDIt, LocalReds, BOp,
-            [&](size_t) { return getReducerAccess(Reducer).getElement(E); });
+      using Name = __sycl_reduction_kernel<reduction::MainKrn, KernelName,
+                                           reduction::strategy::basic,
+                                           decltype(KernelTag)>;
 
-        // Compute the partial sum/reduction for the work-group.
-        if (LID == 0) {
-          size_t GrID = NDIt.get_group_linear_id();
-          typename Reduction::result_type PSum = LocalReds[0];
-          if (IsUpdateOfUserVar)
-            PSum = BOp(Out[0], PSum);
-          Out[GrID * NElements + E] = PSum;
-        }
+      CGH.parallel_for<Name>(NDRange, Properties, [=](nd_item<Dims> NDIt) {
+        // Call user's functions. Reducer.MValue gets initialized there.
+        typename Reduction::reducer_type Reducer =
+            typename Reduction::reducer_type(IdentityContainer, BOp);
+        KernelFunc(NDIt, Reducer);
+
+        size_t WGSize = NDIt.get_local_range().size();
+        size_t LID = NDIt.get_local_linear_id();
 
-        // Ensure item 0 is finished with LocalReds before next iteration
-        if (E != NElements - 1) {
-          NDIt.barrier();
+        auto ElementCombiner = [&](element_type &LHS, const element_type &RHS) {
+          return LHS.combine(BOp, RHS);
+        };
+
+        // If there are multiple values, reduce each separately
+        // This prevents local memory from scaling with elements
+        for (int E = 0; E < NElements; ++E) {
+
+          doTreeReduction<WorkSizeGuarantees::Equal>(
+              WGSize, NDIt, LocalReds, ElementCombiner,
+              [&](size_t) { return getReducerAccess(Reducer).getElement(E); });
+
+          // Compute the partial sum/reduction for the work-group.
+          if (LID == 0) {
+            element_type PSum = LocalReds[0];
+            if constexpr (IsOneWG) {
+              if (IsUpdateOfUserVar)
+                PSum.combine(BOp, Out[E]);
+              Out[E] = *PSum;
+            } else {
+              size_t GrID = NDIt.get_group_linear_id();
+              Out[GrID * NElements + E] = PSum;
+            }
+          }
+
+          // Ensure item 0 is finished with LocalReds before next iteration
+          if (E != NElements - 1) {
+            NDIt.barrier();
+          }
         }
-      }
-    });
+      });
+    };
+
+    if (NWorkGroups == 1)
+      First(KernelOneWGTag{});
+    else
+      First(KernelMultipleWGTag{});
 
     reduction::finalizeHandler(CGH);
 
@@ -1772,76 +1913,97 @@ template <> struct NDRangeReduction<reduction::strategy::basic> {
                                 PI_ERROR_INVALID_WORK_GROUP_SIZE);
     size_t NWorkItems = NDRange.get_group_range().size();
     while (NWorkItems > 1) {
-      reduction::withAuxHandler(CGH, [&](handler &AuxHandler) {
-        size_t NElements = Reduction::num_elements;
-        size_t NWorkGroups;
-        size_t WGSize = reduComputeWGSize(NWorkItems, MaxWGSize, NWorkGroups);
-
-        // The last work-group may be not fully loaded with work, or the work
-        // group size may be not power of two. Those two cases considered
-        // inefficient as they require additional code and checks in the kernel.
-        bool HasUniformWG = NWorkGroups * WGSize == NWorkItems;
-
-        // Get read accessor to the buffer that was used as output
-        // in the previous kernel.
-        auto In = Redu.getReadAccToPreviousPartialReds(AuxHandler);
-        auto Out =
-            Redu.getWriteAccForPartialReds(NWorkGroups * NElements, AuxHandler);
-
-        bool IsUpdateOfUserVar = !Reduction::is_usm &&
-                                 !Redu.initializeToIdentity() &&
-                                 NWorkGroups == 1;
-
-        bool UniformPow2WG = HasUniformWG && (WGSize & (WGSize - 1)) == 0;
-        // Use local memory to reduce elements in work-groups into 0-th element.
-        local_accessor<typename Reduction::result_type, 1> LocalReds{
-            WGSize, AuxHandler};
-
-        auto BOp = Redu.getBinaryOperation();
-        using Name = __sycl_reduction_kernel<reduction::AuxKrn, KernelName,
-                                             reduction::strategy::basic>;
-
-        range<1> GlobalRange = {UniformPow2WG ? NWorkItems
-                                              : NWorkGroups * WGSize};
-        nd_range<1> Range{GlobalRange, range<1>(WGSize)};
-        AuxHandler.parallel_for<Name>(Range, [=](nd_item<1> NDIt) {
-          size_t WGSize = NDIt.get_local_range().size();
-          size_t LID = NDIt.get_local_linear_id();
-          size_t GID = NDIt.get_global_linear_id();
-          size_t GrID = NDIt.get_group_linear_id();
+      size_t NWorkGroups;
+      size_t WGSize = reduComputeWGSize(NWorkItems, MaxWGSize, NWorkGroups);
+
+      auto Rest = [&](auto KernelTag) {
+        reduction::withAuxHandler(CGH, [&](handler &AuxHandler) {
+          // We can deduce IsOneWG from the tag type.
+          constexpr bool IsOneWG =
+              std::is_same_v<std::remove_reference_t<decltype(KernelTag)>,
+                             KernelOneWGTag>;
+
+          constexpr size_t NElements = Reduction::num_elements;
+
+          // The last work-group may be not fully loaded with work, or the work
+          // group size may be not power of two. Those two cases considered
+          // inefficient as they require additional code and checks in the
+          // kernel.
+          bool HasUniformWG = NWorkGroups * WGSize == NWorkItems;
+
+          // Get read accessor to the buffer that was used as output
+          // in the previous kernel.
+          auto In = Redu.getReadAccToPreviousPartialReds(AuxHandler);
+
+          auto Out = [&]() {
+            if constexpr (IsOneWG)
+              return Redu.getUserRedVarAccess(AuxHandler);
+            else
+              return Redu.getWriteAccForPartialReds(NWorkGroups * NElements,
+                                                    AuxHandler);
+          }();
+
+          bool UniformPow2WG = HasUniformWG && (WGSize & (WGSize - 1)) == 0;
+          // Use local memory to reduce elements in work-groups into 0-th
+          // element.
+          local_accessor<element_type, 1> LocalReds{WGSize, AuxHandler};
+
+          auto BOp = Redu.getBinaryOperation();
+          using Name = __sycl_reduction_kernel<reduction::AuxKrn, KernelName,
+                                               reduction::strategy::basic,
+                                               decltype(KernelTag)>;
+
+          range<1> GlobalRange = {UniformPow2WG ? NWorkItems
+                                                : NWorkGroups * WGSize};
+          nd_range<1> Range{GlobalRange, range<1>(WGSize)};
+          AuxHandler.parallel_for<Name>(Range, [=](nd_item<1> NDIt) {
+            size_t WGSize = NDIt.get_local_range().size();
+            size_t LID = NDIt.get_local_linear_id();
+            size_t GID = NDIt.get_global_linear_id();
+            size_t GrID = NDIt.get_group_linear_id();
 
-          for (int E = 0; E < NElements; ++E) {
             // The last work-group may not have enough work for all its items.
             size_t RemainingWorkSize =
                 sycl::min(WGSize, NWorkItems - GrID * WGSize);
 
-            doTreeReduction<WorkSizeGuarantees::LessOrEqual>(
-                RemainingWorkSize, NDIt, LocalReds, BOp,
-                [&](size_t) { return In[GID * NElements + E]; });
-
-            // Compute the partial sum/reduction for the work-group.
-            if (LID == 0) {
-              typename Reduction::result_type PSum = LocalReds[0];
-              if (IsUpdateOfUserVar)
-                PSum = BOp(Out[0], PSum);
-              Out[GrID * NElements + E] = PSum;
-            }
-
-            // Ensure item 0 is finished with LocalReds before next iteration
-            if (E != NElements - 1) {
-              NDIt.barrier();
+            auto ElementCombiner = [&](element_type &LHS,
+                                       const element_type &RHS) {
+              return LHS.combine(BOp, RHS);
+            };
+
+            for (int E = 0; E < NElements; ++E) {
+
+              doTreeReduction<WorkSizeGuarantees::LessOrEqual>(
+                  RemainingWorkSize, NDIt, LocalReds, ElementCombiner,
+                  [&](size_t) { return In[GID * NElements + E]; });
+
+              // Compute the partial sum/reduction for the work-group.
+              if (LID == 0) {
+                element_type PSum = LocalReds[0];
+                if constexpr (IsOneWG) {
+                  if (IsUpdateOfUserVar)
+                    PSum.combine(BOp, Out[E]);
+                  Out[E] = *PSum;
+                } else {
+                  Out[GrID * NElements + E] = PSum;
+                }
+              }
+
+              // Ensure item 0 is finished with LocalReds before next iteration
+              if (E != NElements - 1) {
+                NDIt.barrier();
+              }
             }
-          }
+          });
+          NWorkItems = NWorkGroups;
         });
-        NWorkItems = NWorkGroups;
-      });
-    } // end while (NWorkItems > 1)
+      };
 
-    if constexpr (Reduction::is_usm) {
-      reduction::withAuxHandler(CGH, [&](handler &CopyHandler) {
-        reduSaveFinalResultToUserMem<KernelName>(CopyHandler, Redu);
-      });
-    }
+      if (NWorkGroups == 1)
+        Rest(KernelOneWGTag{});
+      else
+        Rest(KernelMultipleWGTag{});
+    } // end while (NWorkItems > 1)
   }
 };
 
@@ -1873,25 +2035,28 @@ auto getLastCombine(OutAccT OutAcc, LocalAccT LocalAcc, BOPT BOP,
   }
 }
 
-template <typename... Reductions, typename... OutAccT, typename... LocalAccT,
-          typename... BOPsT, typename... Ts, size_t... Is>
+template <bool IsOneWG, typename... Reductions, typename... OutAccT,
+          typename... LocalAccT, typename... BOPsT, typename... Ts,
+          size_t... Is>
 void writeReduSumsToOutAccs(
-    bool IsOneWG, size_t OutAccIndex, ReduTupleT<OutAccT...> OutAccs,
+    size_t OutAccIndex, ReduTupleT<OutAccT...> OutAccs,
     ReduTupleT<LocalAccT...> LocalAccs, ReduTupleT<BOPsT...> BOPs,
     ReduTupleT<Ts...> IdentityVals,
     std::array<bool, sizeof...(Reductions)> IsInitializeToIdentity,
     std::index_sequence<Is...>) {
-  // Add the initial value of user's variable to the final result.
-  if (IsOneWG) {
+  if constexpr (IsOneWG) {
+    // Add the initial value of user's variable to the final result.
+    // After this we know there will be a value in the 0th element.
     ((std::get<Is>(LocalAccs)[0] = getLastCombine(
           std::get<Is>(OutAccs), std::get<Is>(LocalAccs), std::get<Is>(BOPs),
           std::get<Is>(IdentityVals), IsInitializeToIdentity[Is])),
      ...);
+    ((std::get<Is>(OutAccs)[OutAccIndex] = *std::get<Is>(LocalAccs)[0]), ...);
+  } else {
+    // The partial sums for the work-group are stored in 0-th elements of local
+    // accessors. Simply write those sums to output accessors.
+    ((std::get<Is>(OutAccs)[OutAccIndex] = std::get<Is>(LocalAccs)[0]), ...);
   }
-
-  // The partial sums for the work-group are stored in 0-th elements of local
-  // accessors. Simply write those sums to output accessors.
-  ((std::get<Is>(OutAccs)[OutAccIndex] = std::get<Is>(LocalAccs)[0]), ...);
 }
 
 // Concatenate an empty sequence.
@@ -1977,13 +2142,35 @@ struct IsArrayReduction {
   };
 };
 
+template <typename ElementType, typename BOPT>
+constexpr auto makeAdjustedBOP(BOPT &BOP) {
+  return [&](ElementType &LHS, const ElementType &RHS) {
+    return LHS.combine(BOP, RHS);
+  };
+}
+
+template <typename... Reductions, typename... BOPsT, size_t... Is>
+constexpr auto makeAdjustedBOPs(ReduTupleT<BOPsT...> &BOPsTuple,
+                                std::index_sequence<Is...>) {
+  return makeReduTupleT(
+      makeAdjustedBOP<typename std::tuple_element_t<
+          Is, std::tuple<Reductions...>>::reducer_element_type>(
+          std::get<Is>(BOPsTuple))...);
+}
+
+template <typename... Reductions, typename... BOPsT>
+constexpr auto makeAdjustedBOPs(ReduTupleT<BOPsT...> &BOPsTuple) {
+  return makeAdjustedBOPs<Reductions...>(
+      BOPsTuple, std::make_index_sequence<sizeof...(Reductions)>{});
+}
+
 /// All scalar reductions are processed together; there is one loop of log2(N)
 /// steps, and each reduction uses its own storage.
-template <typename... Reductions, int Dims, typename... LocalAccT,
+template <bool IsOneWG, typename... Reductions, int Dims, typename... LocalAccT,
           typename... OutAccT, typename... ReducerT, typename... Ts,
           typename... BOPsT, size_t... Is>
 void reduCGFuncImplScalar(
-    bool IsOneWG, nd_item<Dims> NDIt, ReduTupleT<LocalAccT...> LocalAccsTuple,
+    nd_item<Dims> NDIt, ReduTupleT<LocalAccT...> LocalAccsTuple,
     ReduTupleT<OutAccT...> OutAccsTuple, std::tuple<ReducerT...> &ReducersTuple,
     ReduTupleT<Ts...> IdentitiesTuple, ReduTupleT<BOPsT...> BOPsTuple,
     std::array<bool, sizeof...(Reductions)> InitToIdentityProps,
@@ -1995,45 +2182,62 @@ void reduCGFuncImplScalar(
         getReducerAccess(std::get<Is>(ReducersTuple)).getElement(0)),
    ...);
 
-  doTreeReductionOnTuple(WGSize, LID, LocalAccsTuple, BOPsTuple, ReduIndices);
+  // We apply tree-reduction on reducer elements so we adjust the operations
+  // to combine these.
+  auto AdjustedBOPsTuple = makeAdjustedBOPs<Reductions...>(BOPsTuple);
+
+  doTreeReductionOnTuple(WGSize, LID, LocalAccsTuple, AdjustedBOPsTuple,
+                         ReduIndices);
 
   // Compute the partial sum/reduction for the work-group.
   if (LID == 0) {
     size_t GrID = NDIt.get_group_linear_id();
-    writeReduSumsToOutAccs<Reductions...>(
-        IsOneWG, GrID, OutAccsTuple, LocalAccsTuple, BOPsTuple, IdentitiesTuple,
+    writeReduSumsToOutAccs<IsOneWG, Reductions...>(
+        GrID, OutAccsTuple, LocalAccsTuple, AdjustedBOPsTuple, IdentitiesTuple,
         InitToIdentityProps, ReduIndices);
   }
 }
 
 /// Each array reduction is processed separately.
-template <typename Reduction, int Dims, typename LocalAccT, typename OutAccT,
-          typename ReducerT, typename BOPT>
-void reduCGFuncImplArrayHelper(bool IsOneWG, nd_item<Dims> NDIt,
-                               LocalAccT LocalReds, OutAccT Out,
-                               ReducerT &Reducer, BOPT BOp,
+template <bool IsOneWG, typename Reduction, int Dims, typename LocalAccT,
+          typename OutAccT, typename ReducerT, typename BOPT>
+void reduCGFuncImplArrayHelper(nd_item<Dims> NDIt, LocalAccT LocalReds,
+                               OutAccT Out, ReducerT &Reducer, BOPT BOp,
                                bool IsInitializeToIdentity) {
+  using element_type = typename Reduction::reducer_element_type;
+
   size_t WGSize = NDIt.get_local_range().size();
   size_t LID = NDIt.get_local_linear_id();
 
+  auto ElementCombiner = [&](element_type &LHS, const element_type &RHS) {
+    return LHS.combine(BOp, RHS);
+  };
+
   // If there are multiple values, reduce each separately
   // This prevents local memory from scaling with elements
   auto NElements = Reduction::num_elements;
   for (size_t E = 0; E < NElements; ++E) {
     doTreeReduction<WorkSizeGuarantees::Equal>(
-        WGSize, NDIt, LocalReds, BOp,
+        WGSize, NDIt, LocalReds, ElementCombiner,
         [&](size_t) { return getReducerAccess(Reducer).getElement(E); });
 
     // Add the initial value of user's variable to the final result.
     if (LID == 0) {
       size_t GrID = NDIt.get_group_linear_id();
-      if constexpr (Reduction::has_identity) {
-        Out[GrID * NElements + E] =
-            IsOneWG ? BOp(LocalReds[0],
-                          IsInitializeToIdentity ? Reducer.identity() : Out[E])
-                    : LocalReds[0];
+      size_t OutIdx = GrID * NElements + E;
+      if constexpr (IsOneWG) {
+        // If there is only a single work-group, the output will be an actual
+        // value rather than a potentially optional value.
+        if constexpr (Reduction::has_identity) {
+          Out[OutIdx] = *ElementCombiner(LocalReds[0], IsInitializeToIdentity
+                                                           ? Reducer.identity()
+                                                           : Out[E]);
+        } else {
+          Out[OutIdx] = *LocalReds[0];
+        }
       } else {
-        Out[GrID * NElements + E] = LocalReds[0];
+        // Otherwise we propagate a potentially optional value.
+        Out[OutIdx] = LocalReds[0];
       }
     }
 
@@ -2044,18 +2248,18 @@ void reduCGFuncImplArrayHelper(bool IsOneWG, nd_item<Dims> NDIt,
   }
 }
 
-template <typename... Reductions, int Dims, typename... LocalAccT,
+template <bool IsOneWG, typename... Reductions, int Dims, typename... LocalAccT,
           typename... OutAccT, typename... ReducerT, typename... BOPsT,
           size_t... Is>
 void reduCGFuncImplArray(
-    bool IsOneWG, nd_item<Dims> NDIt, ReduTupleT<LocalAccT...> LocalAccsTuple,
+    nd_item<Dims> NDIt, ReduTupleT<LocalAccT...> LocalAccsTuple,
     ReduTupleT<OutAccT...> OutAccsTuple, std::tuple<ReducerT...> &ReducersTuple,
     ReduTupleT<BOPsT...> BOPsTuple,
     std::array<bool, sizeof...(Reductions)> InitToIdentityProps,
     std::index_sequence<Is...>) {
   using ReductionPack = std::tuple<Reductions...>;
-  (reduCGFuncImplArrayHelper<std::tuple_element_t<Is, ReductionPack>>(
-       IsOneWG, NDIt, std::get<Is>(LocalAccsTuple), std::get<Is>(OutAccsTuple),
+  (reduCGFuncImplArrayHelper<IsOneWG, std::tuple_element_t<Is, ReductionPack>>(
+       NDIt, std::get<Is>(LocalAccsTuple), std::get<Is>(OutAccsTuple),
        std::get<Is>(ReducersTuple), std::get<Is>(BOPsTuple),
        InitToIdentityProps[Is]),
    ...);
@@ -2087,15 +2291,13 @@ void reduCGFuncMulti(handler &CGH, KernelType KernelFunc,
   auto ArrayIs = filterSequence<Reductions...>(ArrayPredicate, ReduIndices);
 
   auto LocalAccsTuple = makeReduTupleT(
-      local_accessor<typename Reductions::result_type, 1>{WGSize, CGH}...);
-
-  size_t NWorkGroups = Range.get_group_range().size();
-  bool IsOneWG = NWorkGroups == 1;
+      local_accessor<typename Reductions::reducer_element_type, 1>{WGSize,
+                                                                   CGH}...);
 
   // The type of the Out "accessor" differs between scenarios when there is just
   // one WorkGroup and when there are multiple. Use this lambda to write the
   // code just once.
-  auto Rest = [&](auto OutAccsTuple) {
+  auto Rest = [&](auto KernelTag, auto OutAccsTuple) {
     auto IdentitiesTuple =
         makeReduTupleT(std::get<Is>(ReduTuple).getIdentityContainer()...);
     auto BOPsTuple =
@@ -2105,9 +2307,14 @@ void reduCGFuncMulti(handler &CGH, KernelType KernelFunc,
 
     using Name = __sycl_reduction_kernel<reduction::MainKrn, KernelName,
                                          reduction::strategy::multi,
-                                         decltype(OutAccsTuple)>;
+                                         decltype(KernelTag)>;
 
     CGH.parallel_for<Name>(Range, Properties, [=](nd_item<Dims> NDIt) {
+      // We can deduce IsOneWG from the tag type.
+      constexpr bool IsOneWG =
+          std::is_same_v<std::remove_reference_t<decltype(KernelTag)>,
+                         KernelOneWGTag>;
+
       // Pass all reductions to user's lambda in the same order as supplied
       // Each reducer initializes its own storage
       auto ReduIndices = std::index_sequence_for<Reductions...>();
@@ -2123,48 +2330,51 @@ void reduCGFuncMulti(handler &CGH, KernelType KernelFunc,
       // reduCGFuncImplScalar<Reductions...>(NDIt, LocalAccsTuple, OutAccsTuple,
       // ReducersTuple, IdentitiesTuple, BOPsTuple, InitToIdentityProps,
       // ReduIndices);
-      reduCGFuncImplScalar<Reductions...>(
-          IsOneWG, NDIt, LocalAccsTuple, OutAccsTuple, ReducersTuple,
-          IdentitiesTuple, BOPsTuple, InitToIdentityProps, ScalarIs);
+      reduCGFuncImplScalar<IsOneWG, Reductions...>(
+          NDIt, LocalAccsTuple, OutAccsTuple, ReducersTuple, IdentitiesTuple,
+          BOPsTuple, InitToIdentityProps, ScalarIs);
 
       // Combine and write-back the results of any array reductions
       // These are handled separately to minimize temporary storage and account
       // for the fact that each array reduction may have a different number of
       // elements to reduce (i.e. a different extent).
-      reduCGFuncImplArray<Reductions...>(IsOneWG, NDIt, LocalAccsTuple,
-                                         OutAccsTuple, ReducersTuple, BOPsTuple,
-                                         InitToIdentityProps, ArrayIs);
+      reduCGFuncImplArray<IsOneWG, Reductions...>(
+          NDIt, LocalAccsTuple, OutAccsTuple, ReducersTuple, BOPsTuple,
+          InitToIdentityProps, ArrayIs);
     });
   };
 
-  if (IsOneWG)
-    Rest(createReduOutAccs<true>(NWorkGroups, CGH, ReduTuple, ReduIndices));
+  size_t NWorkGroups = Range.get_group_range().size();
+  if (NWorkGroups == 1)
+    Rest(KernelOneWGTag{},
+         createReduOutAccs<true>(NWorkGroups, CGH, ReduTuple, ReduIndices));
   else
-    Rest(createReduOutAccs<false>(NWorkGroups, CGH, ReduTuple, ReduIndices));
+    Rest(KernelMultipleWGTag{},
+         createReduOutAccs<false>(NWorkGroups, CGH, ReduTuple, ReduIndices));
 }
 
+// TODO: Is this still needed?
 template <typename... Reductions, size_t... Is>
 void associateReduAccsWithHandler(handler &CGH,
                                   std::tuple<Reductions...> &ReduTuple,
                                   std::index_sequence<Is...>) {
   auto ProcessOne = [&CGH](auto Redu) {
-    if constexpr (!decltype(Redu)::is_usm) {
-      associateWithHandler(CGH, &Redu.getUserRedVar(), access::target::device);
-    }
+    if constexpr (!decltype(Redu)::is_usm)
+      Redu.getUserRedVarAccess(CGH);
   };
   (ProcessOne(std::get<Is>(ReduTuple)), ...);
 }
 
 /// All scalar reductions are processed together; there is one loop of log2(N)
 /// steps, and each reduction uses its own storage.
-template <typename... Reductions, int Dims, typename... LocalAccT,
+template <bool IsOneWG, typename... Reductions, int Dims, typename... LocalAccT,
           typename... InAccT, typename... OutAccT, typename... Ts,
           typename... BOPsT, size_t... Is>
 void reduAuxCGFuncImplScalar(
-    bool IsOneWG, nd_item<Dims> NDIt, size_t LID, size_t GID,
-    size_t RemainingWorkSize, ReduTupleT<LocalAccT...> LocalAccsTuple,
-    ReduTupleT<InAccT...> InAccsTuple, ReduTupleT<OutAccT...> OutAccsTuple,
-    ReduTupleT<Ts...> IdentitiesTuple, ReduTupleT<BOPsT...> BOPsTuple,
+    nd_item<Dims> NDIt, size_t LID, size_t GID, size_t RemainingWorkSize,
+    ReduTupleT<LocalAccT...> LocalAccsTuple, ReduTupleT<InAccT...> InAccsTuple,
+    ReduTupleT<OutAccT...> OutAccsTuple, ReduTupleT<Ts...> IdentitiesTuple,
+    ReduTupleT<BOPsT...> BOPsTuple,
     std::array<bool, sizeof...(Reductions)> InitToIdentityProps,
     std::index_sequence<Is...> ReduIndices) {
   // The end work-group may have less work than the rest, so we only need to
@@ -2172,45 +2382,59 @@ void reduAuxCGFuncImplScalar(
   if (LID < RemainingWorkSize)
     ((std::get<Is>(LocalAccsTuple)[LID] = std::get<Is>(InAccsTuple)[GID]), ...);
 
-  doTreeReductionOnTuple(RemainingWorkSize, LID, LocalAccsTuple, BOPsTuple,
-                         ReduIndices);
+  // We apply tree-reduction on reducer elements so we adjust the operations
+  // to combine these.
+  auto AdjustedBOPsTuple = makeAdjustedBOPs<Reductions...>(BOPsTuple);
+
+  doTreeReductionOnTuple(RemainingWorkSize, LID, LocalAccsTuple,
+                         AdjustedBOPsTuple, ReduIndices);
 
   // Compute the partial sum/reduction for the work-group.
   if (LID == 0) {
     size_t GrID = NDIt.get_group_linear_id();
-    writeReduSumsToOutAccs<Reductions...>(
-        IsOneWG, GrID, OutAccsTuple, LocalAccsTuple, BOPsTuple, IdentitiesTuple,
+    writeReduSumsToOutAccs<IsOneWG, Reductions...>(
+        GrID, OutAccsTuple, LocalAccsTuple, AdjustedBOPsTuple, IdentitiesTuple,
         InitToIdentityProps, ReduIndices);
   }
 }
 
-template <typename Reduction, int Dims, typename LocalAccT, typename InAccT,
-          typename OutAccT, typename T, typename BOPT>
-void reduAuxCGFuncImplArrayHelper(bool IsOneWG, nd_item<Dims> NDIt, size_t LID,
-                                  size_t GID, size_t RemainingWorkSize,
-                                  LocalAccT LocalReds, InAccT In, OutAccT Out,
-                                  T IdentityContainer, BOPT BOp,
-                                  bool IsInitializeToIdentity) {
+template <bool IsOneWG, typename Reduction, int Dims, typename LocalAccT,
+          typename InAccT, typename OutAccT, typename T, typename BOPT>
+void reduAuxCGFuncImplArrayHelper(nd_item<Dims> NDIt, size_t LID, size_t GID,
+                                  size_t RemainingWorkSize, LocalAccT LocalReds,
+                                  InAccT In, OutAccT Out, T IdentityContainer,
+                                  BOPT BOp, bool IsInitializeToIdentity) {
+  using element_type = typename Reduction::reducer_element_type;
+  auto ElementCombiner = [&](element_type &LHS, const element_type &RHS) {
+    return LHS.combine(BOp, RHS);
+  };
 
   // If there are multiple values, reduce each separately
   // This prevents local memory from scaling with elements
   auto NElements = Reduction::num_elements;
   for (size_t E = 0; E < NElements; ++E) {
     doTreeReduction<WorkSizeGuarantees::LessOrEqual>(
-        RemainingWorkSize, NDIt, LocalReds, BOp,
+        RemainingWorkSize, NDIt, LocalReds, ElementCombiner,
         [&](size_t) { return In[GID * NElements + E]; });
 
     // Add the initial value of user's variable to the final result.
     if (LID == 0) {
       size_t GrID = NDIt.get_group_linear_id();
-      if constexpr (Reduction::has_identity) {
-        Out[GrID * NElements + E] =
-            IsOneWG ? BOp(LocalReds[0], IsInitializeToIdentity
-                                            ? IdentityContainer.getIdentity()
-                                            : Out[E])
-                    : LocalReds[0];
+      size_t OutIdx = GrID * NElements + E;
+      if constexpr (IsOneWG) {
+        // If there is only a single work-group, the output will be an actual
+        // value rather than a potentially optional value.
+        if constexpr (Reduction::has_identity) {
+          Out[OutIdx] = *ElementCombiner(LocalReds[0],
+                                         IsInitializeToIdentity
+                                             ? IdentityContainer.getIdentity()
+                                             : Out[E]);
+        } else {
+          Out[OutIdx] = *LocalReds[0];
+        }
       } else {
-        Out[GrID * NElements + E] = LocalReds[0];
+        // Otherwise we propagate a potentially optional value.
+        Out[OutIdx] = LocalReds[0];
       }
     }
 
@@ -2221,19 +2445,20 @@ void reduAuxCGFuncImplArrayHelper(bool IsOneWG, nd_item<Dims> NDIt, size_t LID,
   }
 }
 
-template <typename... Reductions, int Dims, typename... LocalAccT,
+template <bool IsOneWG, typename... Reductions, int Dims, typename... LocalAccT,
           typename... InAccT, typename... OutAccT, typename... Ts,
           typename... BOPsT, size_t... Is>
 void reduAuxCGFuncImplArray(
-    bool IsOneWG, nd_item<Dims> NDIt, size_t LID, size_t GID,
-    size_t RemainingWorkSize, ReduTupleT<LocalAccT...> LocalAccsTuple,
-    ReduTupleT<InAccT...> InAccsTuple, ReduTupleT<OutAccT...> OutAccsTuple,
-    ReduTupleT<Ts...> IdentitiesTuple, ReduTupleT<BOPsT...> BOPsTuple,
+    nd_item<Dims> NDIt, size_t LID, size_t GID, size_t RemainingWorkSize,
+    ReduTupleT<LocalAccT...> LocalAccsTuple, ReduTupleT<InAccT...> InAccsTuple,
+    ReduTupleT<OutAccT...> OutAccsTuple, ReduTupleT<Ts...> IdentitiesTuple,
+    ReduTupleT<BOPsT...> BOPsTuple,
     std::array<bool, sizeof...(Reductions)> InitToIdentityProps,
     std::index_sequence<Is...>) {
   using ReductionPack = std::tuple<Reductions...>;
-  (reduAuxCGFuncImplArrayHelper<std::tuple_element_t<Is, ReductionPack>>(
-       IsOneWG, NDIt, LID, GID, RemainingWorkSize, std::get<Is>(LocalAccsTuple),
+  (reduAuxCGFuncImplArrayHelper<IsOneWG,
+                                std::tuple_element_t<Is, ReductionPack>>(
+       NDIt, LID, GID, RemainingWorkSize, std::get<Is>(LocalAccsTuple),
        std::get<Is>(InAccsTuple), std::get<Is>(OutAccsTuple),
        std::get<Is>(IdentitiesTuple), std::get<Is>(BOPsTuple),
        InitToIdentityProps[Is]),
@@ -2252,7 +2477,6 @@ size_t reduAuxCGFunc(handler &CGH, size_t NWorkItems, size_t MaxWGSize,
   size_t WGSize = reduComputeWGSize(NWorkItems, MaxWGSize, NWorkGroups);
 
   bool Pow2WG = (WGSize & (WGSize - 1)) == 0;
-  bool IsOneWG = NWorkGroups == 1;
   bool HasUniformWG = Pow2WG && (NWorkGroups * WGSize == NWorkItems);
 
   // Like reduCGFuncImpl, we also have to split out scalar and array reductions
@@ -2263,7 +2487,8 @@ size_t reduAuxCGFunc(handler &CGH, size_t NWorkItems, size_t MaxWGSize,
   auto ArrayIs = filterSequence<Reductions...>(ArrayPredicate, ReduIndices);
 
   auto LocalAccsTuple = makeReduTupleT(
-      local_accessor<typename Reductions::result_type, 1>{WGSize, CGH}...);
+      local_accessor<typename Reductions::reducer_element_type, 1>{WGSize,
+                                                                   CGH}...);
   auto InAccsTuple = makeReduTupleT(
       std::get<Is>(ReduTuple).getReadAccToPreviousPartialReds(CGH)...);
 
@@ -2287,6 +2512,11 @@ size_t reduAuxCGFunc(handler &CGH, size_t NWorkItems, size_t MaxWGSize,
     range<1> GlobalRange = {HasUniformWG ? NWorkItems : NWorkGroups * WGSize};
     nd_range<1> Range{GlobalRange, range<1>(WGSize)};
     CGH.parallel_for<Name>(Range, [=](nd_item<1> NDIt) {
+      // We can deduce IsOneWG from the predicate type.
+      constexpr bool IsOneWG =
+          std::is_same_v<std::remove_reference_t<decltype(Predicate)>,
+                         IsNonUsmReductionPredicate>;
+
       size_t WGSize = NDIt.get_local_range().size();
       size_t RemainingWorkSize =
           sycl::min(WGSize, NWorkItems - WGSize * NDIt.get_group_linear_id());
@@ -2294,14 +2524,14 @@ size_t reduAuxCGFunc(handler &CGH, size_t NWorkItems, size_t MaxWGSize,
       size_t GID = NDIt.get_global_linear_id();
 
       // Handle scalar and array reductions
-      reduAuxCGFuncImplScalar<Reductions...>(
-          IsOneWG, NDIt, LID, GID, RemainingWorkSize, LocalAccsTuple,
-          InAccsTuple, OutAccsTuple, IdentitiesTuple, BOPsTuple,
-          InitToIdentityProps, ScalarIs);
-      reduAuxCGFuncImplArray<Reductions...>(
-          IsOneWG, NDIt, LID, GID, RemainingWorkSize, LocalAccsTuple,
-          InAccsTuple, OutAccsTuple, IdentitiesTuple, BOPsTuple,
-          InitToIdentityProps, ArrayIs);
+      reduAuxCGFuncImplScalar<IsOneWG, Reductions...>(
+          NDIt, LID, GID, RemainingWorkSize, LocalAccsTuple, InAccsTuple,
+          OutAccsTuple, IdentitiesTuple, BOPsTuple, InitToIdentityProps,
+          ScalarIs);
+      reduAuxCGFuncImplArray<IsOneWG, Reductions...>(
+          NDIt, LID, GID, RemainingWorkSize, LocalAccsTuple, InAccsTuple,
+          OutAccsTuple, IdentitiesTuple, BOPsTuple, InitToIdentityProps,
+          ArrayIs);
     });
   };
   if (NWorkGroups == 1)
@@ -2564,7 +2794,7 @@ auto reduction(buffer<T, 1, AllocatorT> Var, handler &CGH,
   bool InitializeToIdentity =
       PropList.has_property<property::reduction::initialize_to_identity>();
   return detail::make_reduction<BinaryOperation, 0, 1, false>(
-      accessor{Var, CGH}, Combiner, InitializeToIdentity);
+      Var, Combiner, InitializeToIdentity);
 }
 
 /// Constructs a reduction object using the reduction variable referenced by
@@ -2590,7 +2820,7 @@ auto reduction(buffer<T, 1, AllocatorT> Var, handler &CGH, const T &Identity,
   bool InitializeToIdentity =
       PropList.has_property<property::reduction::initialize_to_identity>();
   return detail::make_reduction<BinaryOperation, 0, 1, true>(
-      accessor{Var, CGH}, Identity, Combiner, InitializeToIdentity);
+      Var, Identity, Combiner, InitializeToIdentity);
 }
 
 /// Constructs a reduction object using the reduction variable referenced by
diff --git a/sycl/test-e2e/Reduction/reduction_big_data.cpp b/sycl/test-e2e/Reduction/reduction_big_data.cpp
index 7146975af534d..ff07266b7e34f 100644
--- a/sycl/test-e2e/Reduction/reduction_big_data.cpp
+++ b/sycl/test-e2e/Reduction/reduction_big_data.cpp
@@ -49,8 +49,9 @@ int test(queue &Q, T Identity) {
 
   // Initialize.
   BinaryOperation BOp;
-  T CorrectOut;
-  initInputData(InBuf, CorrectOut, BOp, NWorkItems);
+  std::optional<T> CorrectOutOpt;
+  initInputData(InBuf, CorrectOutOpt, BOp, NWorkItems);
+  T CorrectOut = *CorrectOutOpt;
 
   // Compute.
   Q.submit([&](handler &CGH) {
diff --git a/sycl/test-e2e/Reduction/reduction_ctor.cpp b/sycl/test-e2e/Reduction/reduction_ctor.cpp
index 6ee4fa48b3525..d7eb0f9c3ba8d 100644
--- a/sycl/test-e2e/Reduction/reduction_ctor.cpp
+++ b/sycl/test-e2e/Reduction/reduction_ctor.cpp
@@ -25,7 +25,7 @@ void test_reducer(Reduction &Redu, T A, T B) {
 
   typename Reduction::binary_operation BOp;
   T ExpectedValue = BOp(A, B);
-  assert(ExpectedValue == detail::ReducerAccess{Reducer}.getElement(0) &&
+  assert(ExpectedValue == *detail::ReducerAccess{Reducer}.getElement(0) &&
          "Wrong result of binary operation.");
   assert(
       toBool(Reducer.identity() == Redu.getIdentityContainer().getIdentity()) &&
@@ -40,7 +40,7 @@ void test_reducer(Reduction &Redu, T Identity, BinaryOperation BOp, T A, T B) {
 
   T ExpectedValue = BOp(A, B);
   assert(
-      toBool(ExpectedValue == detail::ReducerAccess{Reducer}.getElement(0)) &&
+      toBool(ExpectedValue == *detail::ReducerAccess{Reducer}.getElement(0)) &&
       "Wrong result of binary operation.");
   assert(
       toBool(Reducer.identity() == Redu.getIdentityContainer().getIdentity()) &&
diff --git a/sycl/test-e2e/Reduction/reduction_nd_N_vars.cpp b/sycl/test-e2e/Reduction/reduction_nd_N_vars.cpp
index 5b72d6ce3464c..0b8e04409f065 100644
--- a/sycl/test-e2e/Reduction/reduction_nd_N_vars.cpp
+++ b/sycl/test-e2e/Reduction/reduction_nd_N_vars.cpp
@@ -48,7 +48,9 @@ struct Red {
   }
 
   void init() {
-    initInputData(InBuf, CorrectOut, BOp, NWorkItems);
+    std::optional<T> CorrectOutOpt;
+    initInputData(InBuf, CorrectOutOpt, BOp, NWorkItems);
+    CorrectOut = *CorrectOutOpt;
     if (!PropList.template has_property<
             property::reduction::initialize_to_identity>())
       CorrectOut = BOp(CorrectOut, InitVal);
diff --git a/sycl/test-e2e/Reduction/reduction_nd_reducer_skip.cpp b/sycl/test-e2e/Reduction/reduction_nd_reducer_skip.cpp
new file mode 100644
index 0000000000000..cf9a61d2ea0af
--- /dev/null
+++ b/sycl/test-e2e/Reduction/reduction_nd_reducer_skip.cpp
@@ -0,0 +1,66 @@
+// RUN: %clangxx -fsycl -fsycl-targets=%sycl_triple %s -o %t.out
+// RUN: %CPU_RUN_PLACEHOLDER %t.out
+// RUN: %GPU_RUN_PLACEHOLDER %t.out
+// RUN: %ACC_RUN_PLACEHOLDER %t.out
+//
+// Group algorithms are not supported on Nvidia.
+// XFAIL: hip_nvidia
+
+// This test performs basic checks of parallel_for(nd_range, reduction, func)
+// with reductions initialized with a one element buffer. Additionally, some
+// reducers will not be written to.
+
+#include "reduction_utils.hpp"
+
+using namespace sycl;
+
+int NumErrors = 0;
+
+template <typename T> class SkipEvenName;
+template <typename T> class SkipOddName;
+template <typename T> class SkipAllName;
+
+template <typename Name, typename T, class BinaryOperation>
+void tests(queue &Q, T Identity, T Init, BinaryOperation BOp, size_t WGSize,
+           size_t NWItems) {
+  nd_range<1> NDRange(range<1>{NWItems}, range<1>{WGSize});
+  NumErrors += test<SkipEvenName<Name>, T>(Q, Identity, Init, BOp, NDRange,
+                                           property_list{}, SkipEvenOp{});
+  NumErrors += test<SkipOddName<Name>, T>(Q, Identity, Init, BOp, NDRange,
+                                          property_list{}, SkipOddOp{});
+  NumErrors += test<SkipAllName<Name>, T>(Q, Identity, Init, BOp, NDRange,
+                                          property_list{}, SkipAllOp{});
+}
+
+int main() {
+  queue Q;
+  printDeviceInfo(Q);
+
+  // Check some non power-of-two work-group sizes.
+  tests<class A1, int>(Q, 0, 99, std::plus<int>{}, 1, 7);
+  tests<class A2, int>(Q, 0, 99, std::plus<int>{}, 49, 49 * 5);
+
+  // Try some power-of-two work-group sizes.
+  tests<class B1, int>(Q, 0, 99, std::plus<>{}, 1, 32);
+  tests<class B2, int>(Q, 1, 99, std::multiplies<>{}, 4, 32);
+  tests<class B3, int>(Q, 0, 99, std::bit_or<>{}, 8, 128);
+  tests<class B4, int>(Q, 0, 99, std::bit_xor<>{}, 16, 256);
+  tests<class B5, int>(Q, ~0, 99, std::bit_and<>{}, 32, 256);
+  tests<class B6, int>(Q, (std::numeric_limits<int>::max)(), -99,
+                       ext::oneapi::minimum<>{}, 64, 256);
+  tests<class B7, int>(Q, (std::numeric_limits<int>::min)(), 99,
+                       ext::oneapi::maximum<>{}, 128, 256);
+  tests<class B8, int>(Q, 0, 99, std::plus<>{}, 256, 256);
+
+  // Check with various types.
+  tests<class C1, float>(Q, 1, 99, std::multiplies<>{}, 8, 24);
+  tests<class C2, short>(Q, 0x7fff, -99, ext::oneapi::minimum<>{}, 8, 256);
+  tests<class C3, unsigned char>(Q, 0, 99, ext::oneapi::maximum<>{}, 8, 256);
+
+  // Check with CUSTOM type.
+  using CV = CustomVec<long long>;
+  tests<class D1, CV>(Q, CV(0), CV(99), CustomVecPlus<long long>{}, 8, 256);
+
+  printFinalStatus(NumErrors);
+  return NumErrors;
+}
diff --git a/sycl/test-e2e/Reduction/reduction_range_1d_reducer_skip.cpp b/sycl/test-e2e/Reduction/reduction_range_1d_reducer_skip.cpp
new file mode 100644
index 0000000000000..12258afbd4670
--- /dev/null
+++ b/sycl/test-e2e/Reduction/reduction_range_1d_reducer_skip.cpp
@@ -0,0 +1,65 @@
+// RUN: %clangxx -fsycl -fsycl-targets=%sycl_triple %s -o %t.out -Xsycl-target-backend=nvptx64-nvidia-cuda --cuda-gpu-arch=sm_60
+// RUN: %CPU_RUN_PLACEHOLDER %t.out
+// RUN: %GPU_RUN_PLACEHOLDER %t.out
+// RUN: %ACC_RUN_PLACEHOLDER %t.out
+
+// This test performs basic checks of parallel_for(range<1>, reduction, func)
+// with reductions initialized with a one element buffer. Additionally, some
+// reducers will not be written to.
+
+#include "reduction_utils.hpp"
+
+using namespace sycl;
+
+int NumErrors = 0;
+
+template <typename T> class SkipEvenName;
+template <typename T> class SkipOddName;
+template <typename T> class SkipAllName;
+
+template <typename Name, typename T, typename... ArgTys>
+void tests(ArgTys &&...Args) {
+  NumErrors += test<SkipEvenName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                           property_list{}, SkipEvenOp{});
+  NumErrors += test<SkipOddName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipOddOp{});
+  NumErrors += test<SkipAllName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipAllOp{});
+}
+
+int main() {
+  queue Q;
+  printDeviceInfo(Q);
+  size_t MaxWGSize =
+      Q.get_device().get_info<info::device::max_work_group_size>();
+
+  constexpr access::mode RW = access::mode::read_write;
+  // Fast-reduce and Fast-atomics. Try various range types/sizes.
+  tests<class A1, int>(Q, 0, 99, std::plus<int>{}, range<1>(1));
+  tests<class A2, int>(Q, 0, 99, std::plus<>{}, range<1>(2));
+  tests<class A3, int>(Q, 0, 99, std::plus<>{}, range<1>(7));
+  tests<class A4, int>(Q, 0, 99, std::plus<>{}, range<1>(64));
+  tests<class A5, int>(Q, 0, 99, std::plus<>{}, range<1>(MaxWGSize * 2));
+  tests<class A6, int>(Q, 0, 99, std::plus<>{}, range<1>(MaxWGSize * 2 + 5));
+
+  // Check with CUSTOM type.
+  tests<class B1, CustomVec<long long>>(Q, 0, 99, CustomVecPlus<long long>{},
+                                        range<1>(256));
+  tests<class B2, CustomVec<long long>>(Q, 0, 99, CustomVecPlus<long long>{},
+                                        range<1>(MaxWGSize * 3));
+  tests<class B3, CustomVec<long long>>(Q, 99, CustomVecPlus<long long>{},
+                                        range<1>(72));
+
+  // Check with identityless operations.
+  tests<class C1, int>(Q, 99, PlusWithoutIdentity<int>{}, range<1>(1));
+  tests<class C2, int>(Q, 99, PlusWithoutIdentity<int>{}, range<1>(2));
+  tests<class C3, int>(Q, 99, PlusWithoutIdentity<int>{}, range<1>(7));
+  tests<class C4, int>(Q, 99, PlusWithoutIdentity<int>{}, range<1>(64));
+  tests<class C5, int>(Q, 99, PlusWithoutIdentity<int>{},
+                       range<1>(MaxWGSize * 2));
+  tests<class C6, int>(Q, 99, PlusWithoutIdentity<int>{},
+                       range<1>(MaxWGSize * 2 + 5));
+
+  printFinalStatus(NumErrors);
+  return NumErrors;
+}
diff --git a/sycl/test-e2e/Reduction/reduction_range_2d_dw_reducer_skip.cpp b/sycl/test-e2e/Reduction/reduction_range_2d_dw_reducer_skip.cpp
new file mode 100644
index 0000000000000..1598a90bcb63f
--- /dev/null
+++ b/sycl/test-e2e/Reduction/reduction_range_2d_dw_reducer_skip.cpp
@@ -0,0 +1,69 @@
+// RUN: %clangxx -fsycl -fsycl-targets=%sycl_triple %s -o %t.out -Xsycl-target-backend=nvptx64-nvidia-cuda --cuda-gpu-arch=sm_60
+// RUN: %CPU_RUN_PLACEHOLDER %t.out
+// RUN: %GPU_RUN_PLACEHOLDER %t.out
+
+// TODO: accelerator may not suport atomics required by the current
+// implementation. Enable testing when implementation is fixed.
+// RUNx: %ACC_RUN_PLACEHOLDER %t.out
+
+// This test performs basic checks of parallel_for(range<2>, reduction, func)
+// with reductions initialized with a one element buffer. Additionally, some
+// reducers will not be written to.
+
+#include "reduction_utils.hpp"
+
+using namespace sycl;
+
+int NumErrors = 0;
+
+template <typename T> class SkipEvenName;
+template <typename T> class SkipOddName;
+template <typename T> class SkipAllName;
+
+template <typename Name, typename T, typename... ArgTys>
+void tests(ArgTys &&...Args) {
+  NumErrors += test<SkipEvenName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                           property_list{}, SkipEvenOp{});
+  NumErrors += test<SkipOddName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipOddOp{});
+  NumErrors += test<SkipAllName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipAllOp{});
+}
+
+int main() {
+  queue Q;
+  printDeviceInfo(Q);
+  size_t MaxWGSize =
+      Q.get_device().get_info<info::device::max_work_group_size>();
+
+  tests<class A1, int>(Q, 0, 99, std::plus<>{}, range<2>{1, 1});
+  tests<class A2, int>(Q, 0, 99, std::plus<>{}, range<2>{2, 2});
+  tests<class A3, int>(Q, 0, 99, std::plus<>{}, range<2>{2, 3});
+  tests<class A4, int>(Q, 0, 99, std::plus<>{}, range<2>{MaxWGSize, 1});
+  tests<class A5, int64_t>(Q, 0, 99, std::plus<>{}, range<2>{1, MaxWGSize});
+  tests<class A6, int64_t>(Q, 0, 99, std::plus<>{}, range<2>{2, MaxWGSize * 2});
+  tests<class A7, int64_t>(Q, 0, 99, std::plus<>{}, range<2>{MaxWGSize * 3, 7});
+  tests<class A8, int64_t>(Q, 0, 99, std::plus<>{}, range<2>{3, MaxWGSize * 3});
+
+  tests<class B1, CustomVec<long long>>(Q, 0, 99, CustomVecPlus<long long>{},
+                                        range<2>{33, MaxWGSize});
+  tests<class B2, CustomVec<long long>>(Q, 99, CustomVecPlus<long long>{},
+                                        range<2>{33, MaxWGSize});
+
+  tests<class C1, int>(Q, 99, PlusWithoutIdentity<int>{}, range<2>{1, 1});
+  tests<class C2, int>(Q, 99, PlusWithoutIdentity<int>{}, range<2>{2, 2});
+  tests<class C3, int>(Q, 99, PlusWithoutIdentity<int>{}, range<2>{2, 3});
+  tests<class C4, int>(Q, 99, PlusWithoutIdentity<int>{},
+                       range<2>{MaxWGSize, 1});
+  tests<class C5, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<2>{1, MaxWGSize});
+  tests<class C6, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<2>{2, MaxWGSize * 2});
+  tests<class C7, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<2>{MaxWGSize * 3, 7});
+  tests<class C8, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<2>{3, MaxWGSize * 3});
+
+  printFinalStatus(NumErrors);
+  return NumErrors;
+}
diff --git a/sycl/test-e2e/Reduction/reduction_range_3d_rw.cpp b/sycl/test-e2e/Reduction/reduction_range_3d_rw.cpp
index 522669e310b93..f29c828bea9f7 100644
--- a/sycl/test-e2e/Reduction/reduction_range_3d_rw.cpp
+++ b/sycl/test-e2e/Reduction/reduction_range_3d_rw.cpp
@@ -67,14 +67,12 @@ int main() {
   tests<class D2, int>(Q, 99, PlusWithoutIdentity<int>{}, range<3>{2, 2, 2});
   tests<class D3, int>(Q, 99, PlusWithoutIdentity<int>{}, range<3>{2, 3, 4});
 
-  /* Temporarily disabled
   tests<class D4, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
                            range<3>{1, 1, MaxWGSize + 1});
   tests<class D5, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
                            range<3>{1, MaxWGSize + 1, 1});
   tests<class D6, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
                            range<3>{MaxWGSize + 1, 1, 1});
-  */
 
   tests<class D7, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
                            range<3>{2, 5, MaxWGSize * 2});
@@ -83,6 +81,27 @@ int main() {
   tests<class D9, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
                            range<3>{MaxWGSize * 3, 8, 4});
 
+  tests<class E1, int>(Q, 99, MultipliesWithoutIdentity<int>{},
+                       range<3>{1, 1, 1});
+  tests<class E2, int>(Q, 99, MultipliesWithoutIdentity<int>{},
+                       range<3>{2, 2, 2});
+  tests<class E3, int>(Q, 99, MultipliesWithoutIdentity<int>{},
+                       range<3>{2, 3, 4});
+
+  tests<class E4, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{1, 1, MaxWGSize + 1});
+  tests<class E5, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{1, MaxWGSize + 1, 1});
+  tests<class E6, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{MaxWGSize + 1, 1, 1});
+
+  tests<class E7, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{2, 5, MaxWGSize * 2});
+  tests<class E8, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{3, MaxWGSize * 3, 2});
+  tests<class E9, int64_t>(Q, 99, MultipliesWithoutIdentity<int64_t>{},
+                           range<3>{MaxWGSize * 3, 8, 4});
+
   printFinalStatus(NumErrors);
   return NumErrors;
 }
diff --git a/sycl/test-e2e/Reduction/reduction_range_3d_rw_reducer_skip.cpp b/sycl/test-e2e/Reduction/reduction_range_3d_rw_reducer_skip.cpp
new file mode 100644
index 0000000000000..089720b17e917
--- /dev/null
+++ b/sycl/test-e2e/Reduction/reduction_range_3d_rw_reducer_skip.cpp
@@ -0,0 +1,82 @@
+// RUN: %clangxx -fsycl -fsycl-targets=%sycl_triple %s -o %t.out -Xsycl-target-backend=nvptx64-nvidia-cuda --cuda-gpu-arch=sm_60
+// RUN: %CPU_RUN_PLACEHOLDER %t.out
+// RUN: %GPU_RUN_PLACEHOLDER %t.out
+
+// TODO: accelerator may not suport atomics required by the current
+// implementation. Enable testing when implementation is fixed.
+// RUNx: %ACC_RUN_PLACEHOLDER %t.out
+
+// This test performs basic checks of parallel_for(range<3>, reduction, func)
+// with reductions initialized with a one element buffer. Additionally, some
+// reducers will not be written to.
+
+#include "reduction_utils.hpp"
+
+using namespace sycl;
+
+int NumErrors = 0;
+
+template <typename T> class SkipEvenName;
+template <typename T> class SkipOddName;
+template <typename T> class SkipAllName;
+
+template <typename Name, typename T, typename... ArgTys>
+void tests(ArgTys &&...Args) {
+  NumErrors += test<SkipEvenName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                           property_list{}, SkipEvenOp{});
+  NumErrors += test<SkipOddName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipOddOp{});
+  NumErrors += test<SkipAllName<Name>, T>(std::forward<ArgTys>(Args)...,
+                                          property_list{}, SkipAllOp{});
+}
+
+int main() {
+  queue Q;
+  printDeviceInfo(Q);
+  size_t MaxWGSize =
+      Q.get_device().get_info<info::device::max_work_group_size>();
+
+  tests<class A1, int>(Q, 0, 99, std::plus<>{}, range<3>{1, 1, 1});
+  tests<class A2, int>(Q, 0, 99, std::plus<>{}, range<3>{2, 2, 2});
+  tests<class A3, int>(Q, 0, 99, std::plus<>{}, range<3>{2, 3, 4});
+
+  tests<class A4, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{1, 1, MaxWGSize + 1});
+  tests<class A5, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{1, MaxWGSize + 1, 1});
+  tests<class A6, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{MaxWGSize + 1, 1, 1});
+
+  tests<class A7, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{2, 5, MaxWGSize * 2});
+  tests<class A8, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{3, MaxWGSize * 3, 2});
+  tests<class A9, int64_t>(Q, 0, 99, std::plus<>{},
+                           range<3>{MaxWGSize * 3, 8, 4});
+
+  tests<class B1, CustomVec<long long>>(Q, 0, 99, CustomVecPlus<long long>{},
+                                        range<3>{2, 33, MaxWGSize});
+  tests<class B2, CustomVec<long long>>(Q, 99, CustomVecPlus<long long>{},
+                                        range<3>{2, 33, MaxWGSize});
+
+  tests<class C1, int>(Q, 99, PlusWithoutIdentity<int>{}, range<3>{1, 1, 1});
+  tests<class C2, int>(Q, 99, PlusWithoutIdentity<int>{}, range<3>{2, 2, 2});
+  tests<class C3, int>(Q, 99, PlusWithoutIdentity<int>{}, range<3>{2, 3, 4});
+
+  tests<class C4, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{1, 1, MaxWGSize + 1});
+  tests<class C5, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{1, MaxWGSize + 1, 1});
+  tests<class C6, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{MaxWGSize + 1, 1, 1});
+
+  tests<class C7, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{2, 5, MaxWGSize * 2});
+  tests<class C8, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{3, MaxWGSize * 3, 2});
+  tests<class C9, int64_t>(Q, 99, PlusWithoutIdentity<int64_t>{},
+                           range<3>{MaxWGSize * 3, 8, 4});
+
+  printFinalStatus(NumErrors);
+  return NumErrors;
+}
diff --git a/sycl/test-e2e/Reduction/reduction_range_N_vars.cpp b/sycl/test-e2e/Reduction/reduction_range_N_vars.cpp
index 290227208113b..ed7e8d8b83688 100644
--- a/sycl/test-e2e/Reduction/reduction_range_N_vars.cpp
+++ b/sycl/test-e2e/Reduction/reduction_range_N_vars.cpp
@@ -48,7 +48,9 @@ struct Red {
   }
 
   void init() {
-    initInputData(InBuf, CorrectOut, BOp, NWorkItems);
+    std::optional<T> CorrectOutOpt;
+    initInputData(InBuf, CorrectOutOpt, BOp, NWorkItems);
+    CorrectOut = *CorrectOutOpt;
     if (!PropList.template has_property<
             property::reduction::initialize_to_identity>())
       CorrectOut = BOp(CorrectOut, InitVal);
diff --git a/sycl/test-e2e/Reduction/reduction_usm.cpp b/sycl/test-e2e/Reduction/reduction_usm.cpp
index 2a53361b13637..2e6f25d273626 100644
--- a/sycl/test-e2e/Reduction/reduction_usm.cpp
+++ b/sycl/test-e2e/Reduction/reduction_usm.cpp
@@ -38,12 +38,12 @@ int test(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
   }
 
   // Initialize.
-  T CorrectOut;
+  std::optional<T> CorrectOutOpt;
   BinaryOperation BOp;
 
   buffer<T, 1> InBuf(NWItems);
-  initInputData(InBuf, CorrectOut, BOp, NWItems);
-  CorrectOut = BOp(CorrectOut, Init);
+  initInputData(InBuf, CorrectOutOpt, BOp, NWItems);
+  T CorrectOut = BOp(*CorrectOutOpt, Init);
 
   // Compute.
   Q.submit([&](handler &CGH) {
diff --git a/sycl/test-e2e/Reduction/reduction_usm_dw.cpp b/sycl/test-e2e/Reduction/reduction_usm_dw.cpp
index 1460fac3907b5..f94c564e23e68 100644
--- a/sycl/test-e2e/Reduction/reduction_usm_dw.cpp
+++ b/sycl/test-e2e/Reduction/reduction_usm_dw.cpp
@@ -39,11 +39,12 @@ int test(queue &Q, T Identity, T Init, size_t WGSize, size_t NWItems,
   }
 
   // Initialize.
-  T CorrectOut;
+  std::optional<T> CorrectOutOpt;
   BinaryOperation BOp;
 
   buffer<T, 1> InBuf(NWItems);
-  initInputData(InBuf, CorrectOut, BOp, NWItems);
+  initInputData(InBuf, CorrectOutOpt, BOp, NWItems);
+  T CorrectOut = *CorrectOutOpt;
 
   // Compute.
   Q.submit([&](handler &CGH) {
diff --git a/sycl/test-e2e/Reduction/reduction_utils.hpp b/sycl/test-e2e/Reduction/reduction_utils.hpp
index d5875fbbba475..74c0f0913cce1 100644
--- a/sycl/test-e2e/Reduction/reduction_utils.hpp
+++ b/sycl/test-e2e/Reduction/reduction_utils.hpp
@@ -4,11 +4,29 @@
 
 using namespace sycl;
 
+struct AllIdOp {
+  constexpr bool operator()(size_t Idx) const { return true; }
+};
+
+struct SkipAllOp {
+  constexpr bool operator()(size_t Idx) const { return false; }
+};
+
+struct SkipEvenOp {
+  constexpr bool operator()(size_t Idx) const { return Idx % 2; }
+};
+
+struct SkipOddOp {
+  constexpr bool operator()(size_t Idx) const { return (Idx + 1) % 2; }
+};
+
 /// Initializes the buffer<1> \p 'InBuf' buffer with pseudo-random values,
-/// computes the write the reduction value \p 'ExpectedOut'.
-template <typename T, class BinaryOperation>
-void initInputData(buffer<T, 1> &InBuf, T &ExpectedOut, BinaryOperation BOp,
-                   range<1> Range) {
+/// computes the write the reduction value \p 'ExpectedOut'. Linearized IDs are
+/// filtered in \p 'ExpectedOut' using \p 'IdFilterFunc'.
+template <typename T, class BinaryOperation, typename IdFilterFuncT = AllIdOp>
+void initInputData(buffer<T, 1> &InBuf, std::optional<T> &ExpectedOut,
+                   BinaryOperation BOp, range<1> Range,
+                   IdFilterFuncT IdFilterFunc = {}) {
   size_t N = Range.size();
   assert(N != 0);
   auto In = InBuf.template get_access<access::mode::write>();
@@ -27,15 +45,18 @@ void initInputData(buffer<T, 1> &InBuf, T &ExpectedOut, BinaryOperation BOp,
       In[I] = I;
     else
       In[I] = ((I + 1) % 5) + 1.1;
-    ExpectedOut = I == 0 ? In[I] : BOp(ExpectedOut, In[I]);
+    if (IdFilterFunc(I))
+      ExpectedOut = ExpectedOut ? BOp(*ExpectedOut, In[I]) : In[I];
   }
 };
 
 /// Initializes the buffer<2> \p 'InBuf' buffer with pseudo-random values,
-/// computes the write the reduction value \p 'ExpectedOut'.
-template <typename T, class BinaryOperation>
-void initInputData(buffer<T, 2> &InBuf, T &ExpectedOut, BinaryOperation BOp,
-                   range<2> Range) {
+/// computes the write the reduction value \p 'ExpectedOut'. Linearized IDs are
+/// filtered in \p 'ExpectedOut' using \p 'IdFilterFunc'.
+template <typename T, class BinaryOperation, typename IdFilterFuncT = AllIdOp>
+void initInputData(buffer<T, 2> &InBuf, std::optional<T> &ExpectedOut,
+                   BinaryOperation BOp, range<2> Range,
+                   IdFilterFuncT IdFilterFunc = {}) {
   assert(Range.size() != 0);
   auto In = InBuf.template get_access<access::mode::write>();
   for (int J = 0; J < Range[0]; ++J) {
@@ -54,16 +75,19 @@ void initInputData(buffer<T, 2> &InBuf, T &ExpectedOut, BinaryOperation BOp,
         In[J][I] = I + J;
       else
         In[J][I] = ((I + 1 + J) % 5) + 1.1;
-      ExpectedOut = (I == 0 && J == 0) ? In[J][I] : BOp(ExpectedOut, In[J][I]);
+      if (IdFilterFunc(I + J * Range[1]))
+        ExpectedOut = ExpectedOut ? BOp(*ExpectedOut, In[J][I]) : In[J][I];
     }
   }
 };
 
 /// Initializes the buffer<3> \p 'InBuf' buffer with pseudo-random values,
-/// computes the write the reduction value \p 'ExpectedOut'.
-template <typename T, class BinaryOperation>
-void initInputData(buffer<T, 3> &InBuf, T &ExpectedOut, BinaryOperation BOp,
-                   range<3> Range) {
+/// computes the write the reduction value \p 'ExpectedOut'. Linearized IDs are
+/// filtered in \p 'ExpectedOut' using \p 'IdFilterFunc'.
+template <typename T, class BinaryOperation, typename IdFilterFuncT = AllIdOp>
+void initInputData(buffer<T, 3> &InBuf, std::optional<T> &ExpectedOut,
+                   BinaryOperation BOp, range<3> Range,
+                   IdFilterFuncT IdFilterFunc = {}) {
   assert(Range.size() != 0);
   auto In = InBuf.template get_access<access::mode::write>();
   for (int K = 0; K < Range[0]; ++K) {
@@ -83,9 +107,9 @@ void initInputData(buffer<T, 3> &InBuf, T &ExpectedOut, BinaryOperation BOp,
           In[K][J][I] = I + J + K;
         else
           In[K][J][I] = ((I + 1 + J + K * 3) % 5) + 1.1;
-        ExpectedOut = (I == 0 && J == 0 && K == 0)
-                          ? In[K][J][I]
-                          : BOp(ExpectedOut, In[K][J][I]);
+        if (IdFilterFunc(I + J * Range[2] + K * Range[1] * Range[2]))
+          ExpectedOut =
+              ExpectedOut ? BOp(*ExpectedOut, In[K][J][I]) : In[K][J][I];
       }
     }
   }
@@ -141,6 +165,10 @@ template <class T> struct PlusWithoutIdentity {
   T operator()(const T &A, const T &B) const { return A + B; }
 };
 
+template <class T> struct MultipliesWithoutIdentity {
+  T operator()(const T &A, const T &B) const { return A * B; }
+};
+
 template <typename T> T getMinimumFPValue() {
   return std::numeric_limits<T>::has_infinity
              ? static_cast<T>(-std::numeric_limits<T>::infinity())
@@ -283,10 +311,10 @@ auto init_to_identity() {
 
 template <typename Name, typename T, bool HasIdentity, class BinaryOperation,
           template <int> typename RangeTy, int Dims,
-          typename PropListTy = property_list>
+          typename PropListTy = property_list, typename IdFilterFuncT = AllIdOp>
 int testInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
               BinaryOperation BOp, const RangeTy<Dims> &Range,
-              PropListTy PropList = {}) {
+              PropListTy PropList = {}, IdFilterFuncT IdFilterFunc = {}) {
   constexpr bool IsRange = std::is_same_v<range<Dims>, RangeTy<Dims>>;
   constexpr bool IsNDRange = std::is_same_v<nd_range<Dims>, RangeTy<Dims>>;
   static_assert(IsRange || IsNDRange);
@@ -326,11 +354,11 @@ int testInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
   buffer<T, 1> OutBuf(1);
 
   // Initialize.
-  T CorrectOut;
-  initInputData(InBuf, CorrectOut, BOp, GlobalRange);
+  std::optional<T> CorrectOut;
+  initInputData(InBuf, CorrectOut, BOp, GlobalRange, IdFilterFunc);
   if (!PropList.template has_property<
           property::reduction::initialize_to_identity>()) {
-    CorrectOut = BOp(CorrectOut, Init);
+    CorrectOut = CorrectOut ? BOp(*CorrectOut, Init) : Init;
   }
 
   // The value assigned here must be discarded (if IsReadWrite is true).
@@ -352,42 +380,48 @@ int testInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
     auto In = InBuf.template get_access<access::mode::read>(CGH);
     auto Redu = CreateReduction();
     if constexpr (IsRange)
-      CGH.parallel_for<Name>(
-          Range, Redu, [=](id<Dims> Id, auto &Sum) { Sum.combine(In[Id]); });
+      CGH.parallel_for<Name>(Range, Redu, [=](item<Dims> Id, auto &Sum) {
+        if (IdFilterFunc(Id.get_linear_id()))
+          Sum.combine(In[Id]);
+      });
     else
       CGH.parallel_for<Name>(Range, Redu, [=](nd_item<Dims> NDIt, auto &Sum) {
-        Sum.combine(In[NDIt.get_global_id()]);
+        if (IdFilterFunc(NDIt.get_global_linear_id()))
+          Sum.combine(In[NDIt.get_global_linear_id()]);
       });
   });
 
   // Check correctness.
   auto Out = OutBuf.template get_access<access::mode::read>();
   T ComputedOut = *(Out.get_pointer());
-  return checkResults(Q, BOp, Range, ComputedOut, CorrectOut);
+  return checkResults(Q, BOp, Range, ComputedOut, *CorrectOut);
 }
 
 template <typename Name, typename T, class BinaryOperation,
           template <int> typename RangeTy, int Dims,
-          typename PropListTy = property_list>
+          typename PropListTy = property_list, typename IdFilterFuncT = AllIdOp>
 int test(queue &Q, T Identity, T Init, BinaryOperation BOp,
-         const RangeTy<Dims> &Range, PropListTy PropList = {}) {
+         const RangeTy<Dims> &Range, PropListTy PropList = {},
+         IdFilterFuncT IdFilterFunc = {}) {
   return testInner<Name>(Q, OptionalIdentity(Identity), Init, BOp, Range,
                          PropList);
 }
 
 template <typename Name, typename T, class BinaryOperation,
           template <int> typename RangeTy, int Dims,
-          typename PropListTy = property_list>
+          typename PropListTy = property_list, typename IdFilterFuncT = AllIdOp>
 int test(queue &Q, T Init, BinaryOperation BOp, const RangeTy<Dims> &Range,
-         PropListTy PropList = {}) {
+         PropListTy PropList = {}, IdFilterFuncT IdFilterFunc = {}) {
   return testInner<Name>(Q, OptionalIdentity<T>(), Init, BOp, Range, PropList);
 }
 
 template <typename Name, typename T, bool HasIdentity, class BinaryOperation,
-          int Dims, typename PropListTy = property_list>
+          int Dims, typename PropListTy = property_list,
+          typename IdFilterFuncT = AllIdOp>
 int testUSMInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
                  BinaryOperation BOp, const range<Dims> &Range,
-                 usm::alloc AllocType, PropListTy PropList = {}) {
+                 usm::alloc AllocType, PropListTy PropList = {},
+                 IdFilterFuncT IdFilterFunc = {}) {
   printTestLabel<T, BinaryOperation>(Range);
 
   auto Dev = Q.get_device();
@@ -434,12 +468,12 @@ int testUSMInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
   }
 
   // Initialize.
-  T CorrectOut;
+  std::optional<T> CorrectOut;
   buffer<T, Dims> InBuf(Range);
-  initInputData(InBuf, CorrectOut, BOp, Range);
+  initInputData(InBuf, CorrectOut, BOp, Range, IdFilterFunc);
   if (!PropList.template has_property<
           property::reduction::initialize_to_identity>()) {
-    CorrectOut = BOp(CorrectOut, Init);
+    CorrectOut = CorrectOut ? BOp(*CorrectOut, Init) : Init;
   }
 
   // Compute.
@@ -457,7 +491,10 @@ int testUSMInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
      auto In = InBuf.template get_access<access::mode::read>(CGH);
      auto Redu = CreateReduction();
      CGH.parallel_for<TName<Name, class Test>>(
-         Range, Redu, [=](id<Dims> Id, auto &Sum) { Sum.combine(In[Id]); });
+         Range, Redu, [=](item<Dims> Id, auto &Sum) {
+           if (IdFilterFunc(Id.get_linear_id()))
+             Sum.combine(In[Id]);
+         });
    }).wait();
 
   // Check correctness.
@@ -476,24 +513,25 @@ int testUSMInner(queue &Q, OptionalIdentity<T, HasIdentity> Identity, T Init,
 
   std::string AllocStr =
       "AllocMode=" + std::to_string(static_cast<int>(AllocType));
-  int Error = checkResults(Q, BOp, Range, ComputedOut, CorrectOut, AllocStr);
+  int Error = checkResults(Q, BOp, Range, ComputedOut, *CorrectOut, AllocStr);
   free(ReduVarPtr, Q.get_context());
   return Error;
 }
 
 template <typename Name, typename T, class BinaryOperation, int Dims,
-          typename PropListTy = property_list>
+          typename PropListTy = property_list, typename IdFilterFuncT = AllIdOp>
 int testUSM(queue &Q, T Identity, T Init, BinaryOperation BOp,
             const range<Dims> &Range, usm::alloc AllocType,
-            PropListTy PropList = {}) {
+            property_list PropList = {}, IdFilterFuncT IdFilterFunc = {}) {
   return testUSMInner<Name>(Q, OptionalIdentity(Identity), Init, BOp, Range,
                             AllocType, PropList);
 }
 
 template <typename Name, typename T, class BinaryOperation, int Dims,
-          typename PropListTy = property_list>
+          typename PropListTy = property_list, typename IdFilterFuncT = AllIdOp>
 int testUSM(queue &Q, T Init, BinaryOperation BOp, const range<Dims> &Range,
-            usm::alloc AllocType, PropListTy PropList = {}) {
+            usm::alloc AllocType, property_list PropList = {},
+            IdFilterFuncT IdFilterFunc = {}) {
   return testUSMInner<Name>(Q, OptionalIdentity<T>(), Init, BOp, Range,
                             AllocType, PropList);
-}
+}
\ No newline at end of file