AMReX_Reduce.H

Go to the documentation of this file.

#ifndef AMREX_REDUCE_H_
#define AMREX_REDUCE_H_
#include <AMReX_Config.H>
#include <AMReX.H>
 
#include <AMReX_Concepts.H>
#include <AMReX_Gpu.H>
#include <AMReX_Arena.H>
#include <AMReX_OpenMP.H>
#include <AMReX_MFIter.H>
#include <AMReX_TypeList.H>
#include <AMReX_ValLocPair.H>
 
#include <algorithm>
#include <concepts>
#include <functional>
#include <limits>
#include <utility>
 
namespace amrex {
 
namespace Reduce {
 
// The declaration of these functions are here to work around doxygen issues.

template <typename T, std::integral N>
T Sum (N n, T const* v, T init_val = 0);

template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Sum (N n, F const& f, T init_val = 0);

template <typename T, std::integral N>
T Min (N n, T const* v, T init_val = std::numeric_limits<T>::max());

template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Min (N n, F const& f, T init_val = std::numeric_limits<T>::max());

template <typename T, std::integral N>
T Max (N n, T const* v, T init_val = std::numeric_limits<T>::lowest());

template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Max (N n, F const& f, T init_val = std::numeric_limits<T>::lowest());

template <typename T, std::integral N>
std::pair<T,T> MinMax (N n, T const* v);

template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
std::pair<T,T> MinMax (N n, F const& f);

template <typename T, std::integral N, typename P>
bool AnyOf (N n, T const* v, P const& pred);

template <typename P, int dim>
bool AnyOf (BoxND<dim> const& box, P const& pred);
 
}

namespace Reduce::detail {
 
// CPU reductions keep nacc partial accumulators so that the inner loop
// vectorizes.
inline constexpr int nacc = 8;
 
// Pragmas for the loop over the nacc accumulators. With GCC, `unroll 4`
// keeps it vectorizable and novector plus a full unroll keeps it scalar on
// purpose. Clang needs neither, and nvcc does not know them.
#if defined(AMREX_CXX_GCC) && !defined(__CUDACC__)
#define AMREX_PRAGMA_NACC_VECTORIZE _Pragma("GCC unroll 4")
#if __GNUC__ >= 14
#define AMREX_PRAGMA_NACC_NOVECTOR _Pragma("GCC novector")
#else
#define AMREX_PRAGMA_NACC_NOVECTOR
#endif
#define AMREX_PRAGMA_NACC_SCALAR AMREX_PRAGMA_NACC_NOVECTOR _Pragma("GCC unroll 8")
#else
#define AMREX_PRAGMA_NACC_VECTORIZE
#define AMREX_PRAGMA_NACC_NOVECTOR
#define AMREX_PRAGMA_NACC_SCALAR
#endif
static_assert(nacc == 8, "The AMREX_PRAGMA_NACC unroll factors assume nacc == 8");
 
// Split a 1D range into chunks of a multiple of nacc, enough for all threads.
template <std::integral N>
std::pair<N,N> chunks_1d (N n)
{
    auto const nthreads = OpenMP::get_max_threads();
    N const chunk = std::max(static_cast<N>(nacc), static_cast<N>((n / (4*nthreads) / nacc) * nacc));
    N const nchunks = n / chunk + ((n % chunk != 0) ? 1 : 0);
    return {chunk, nchunks};
}
 
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <std::size_t Idx, typename T, typename P>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void for_each_parallel (T& d, T const& s, Gpu::Handler const& h)
    {
        P().parallel_update(amrex::get<Idx>(d), amrex::get<Idx>(s), h);
    }
 
    template <std::size_t Idx, typename T, typename P, typename P1, typename... Ps>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void for_each_parallel (T& d, T const& s, Gpu::Handler const& h)
    {
        P().parallel_update(amrex::get<Idx>(d), amrex::get<Idx>(s), h);
        for_each_parallel<Idx+1,T,P1,Ps...>(d, s, h);
    }
#else
    template <std::size_t Idx, typename T, typename P>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void for_each_parallel (T& d, T const& s)
    {
        P().parallel_update(amrex::get<Idx>(d), amrex::get<Idx>(s));
    }
 
    template <std::size_t Idx, typename T, typename P, typename P1, typename... Ps>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void for_each_parallel (T& d, T const& s)
    {
        P().parallel_update(amrex::get<Idx>(d), amrex::get<Idx>(s));
        for_each_parallel<Idx+1,T,P1,Ps...>(d, s);
    }
#endif
#endif
 
    template <std::size_t Idx, typename T, typename P>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void for_each_local (T& d, T const& s)
    {
        P().local_update(amrex::get<Idx>(d), amrex::get<Idx>(s));
    }
 
    template <std::size_t Idx, typename T, typename P, typename P1, typename... Ps>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void for_each_local (T& d, T const& s)
    {
        P().local_update(amrex::get<Idx>(d), amrex::get<Idx>(s));
        for_each_local<Idx+1,T,P1,Ps...>(d, s);
    }
 
    template <std::size_t Idx, typename T, typename P>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    constexpr void for_each_init (T& t)
    {
        P().init(amrex::get<Idx>(t));
    }
 
    template <std::size_t Idx, typename T, typename P, typename P1, typename... Ps>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    constexpr void for_each_init (T& t)
    {
        P().init(amrex::get<Idx>(t));
        for_each_init<Idx+1,T,P1,Ps...>(t);
    }
}

struct ReduceOpSum
{
 
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <typename T>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s, Gpu::Handler const& h) const noexcept {
        T r = Gpu::blockReduceSum(s,h);
        if (h.threadIdx() == 0) { d += r; }
    }
#else
    template <typename T, int MT=AMREX_GPU_MAX_THREADS>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s) const noexcept {
        T r = Gpu::blockReduceSum<MT>(s);
        if (threadIdx.x == 0) { d += r; }
    }
#endif
#endif
 
    template <typename T>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void local_update (T& d, T const& s) const noexcept { d += s; }
 
    template <typename T>
    constexpr void init (T& t) const noexcept { t = 0; }
};

struct ReduceOpMin
{
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <typename T>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s, Gpu::Handler const& h) const noexcept {
        T r = Gpu::blockReduceMin(s,h);
        if (h.threadIdx() == 0) { d = amrex::min(d,r); }
    }
#else
    template <typename T, int MT=AMREX_GPU_MAX_THREADS>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s) const noexcept {
        T r = Gpu::blockReduceMin<MT>(s);
        if (threadIdx.x == 0) { d = amrex::min(d,r); }
    }
#endif
#endif
 
    template <typename T>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void local_update (T& d, T const& s) const noexcept { d = amrex::min(d,s); }
 
    template <typename T>
    requires (std::numeric_limits<T>::is_specialized)
    constexpr void init (T& t) const noexcept { t = std::numeric_limits<T>::max(); }
 
    template <typename T>
    requires (!std::numeric_limits<T>::is_specialized)
    constexpr void init (T& t) const noexcept { t = T::max(); }
};

struct ReduceOpMax
{
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <typename T>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s, Gpu::Handler const& h) const noexcept {
        T r = Gpu::blockReduceMax(s,h);
        if (h.threadIdx() == 0) { d = amrex::max(d,r); }
    }
#else
    template <typename T, int MT=AMREX_GPU_MAX_THREADS>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T const& s) const noexcept {
        T r = Gpu::blockReduceMax<MT>(s);
        if (threadIdx.x == 0) { d = amrex::max(d,r); }
    }
#endif
#endif
 
    template <typename T>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void local_update (T& d, T const& s) const noexcept { d = amrex::max(d,s); }
 
    template <typename T>
    requires (std::numeric_limits<T>::is_specialized)
    constexpr void init (T& t) const noexcept { t = std::numeric_limits<T>::lowest(); }
 
    template <typename T>
    requires (!std::numeric_limits<T>::is_specialized)
    constexpr void init (T& t) const noexcept { t = T::lowest(); }
};

struct ReduceOpLogicalAnd
{
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <std::integral T>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T s, Gpu::Handler const& h) const noexcept {
        T r = Gpu::blockReduceLogicalAnd(s,h);
        if (h.threadIdx() == 0) { d = d && r; }
    }
#else
    template <std::integral T, int MT=AMREX_GPU_MAX_THREADS>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T s) const noexcept {
        T r = Gpu::blockReduceLogicalAnd<MT>(s);
        if (threadIdx.x == 0) { d = d && r; }
    }
#endif
#endif
 
    template <std::integral T>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void local_update (T& d, T s) const noexcept { d = d && s; }
 
    template <std::integral T>
    constexpr void init (T& t) const noexcept { t = true; }
};

struct ReduceOpLogicalOr
{
#ifdef AMREX_USE_GPU
#ifdef AMREX_USE_SYCL
    template <std::integral T>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T s, Gpu::Handler const& h) const noexcept {
        T r = Gpu::blockReduceLogicalOr(s,h);
        if (h.threadIdx() == 0) { d = d || r; }
    }
#else
    template <std::integral T, int MT=AMREX_GPU_MAX_THREADS>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void parallel_update (T& d, T s) const noexcept {
        T r = Gpu::blockReduceLogicalOr<MT>(s);
        if (threadIdx.x == 0) { d = d || r; }
    }
#endif
#endif
 
    template <std::integral T>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    void local_update (T& d, T s) const noexcept { d = d || s; }
 
    template <std::integral T>
    constexpr void init (T& t) const noexcept { t = false; }
};

namespace Reduce::detail {
// Whether the accumulator loop of a CPU reduction over Ps is vectorized.
// A lone min or max stays scalar with GCC 14 or newer on x86 without
// SSE4.1, where the vector form is slower. Older GCC has no novector.
template <typename... Ps>
inline constexpr bool vectorize_nacc_loop =
#if defined(AMREX_CXX_GCC) && !defined(__CUDACC__) && (__GNUC__ >= 14) && \
    (defined(__x86_64__) || defined(__i386__)) && !defined(__SSE4_1__)
    sizeof...(Ps) != 1 ||
    !((std::is_same_v<Ps,ReduceOpMin> || std::is_same_v<Ps,ReduceOpMax>) && ...);
#else
    true;
#endif
 
// The nacc partial accumulators of a CPU reduction. GCC keeps them in
// vector registers only when each tuple member has its own array, clang
// only when each lane holds a whole tuple.
template <typename Tup> struct Lanes;
 
template <typename... Ts>
struct Lanes<GpuTuple<Ts...>>
{
#ifdef AMREX_CXX_GCC
    GpuTuple<GpuArray<Ts,nacc>...> data;
    template <int K> auto& get (int lane) noexcept { return amrex::get<K>(data)[lane]; }
    template <int K> [[nodiscard]] auto const& get (int lane) const noexcept { return amrex::get<K>(data)[lane]; }
#else
    GpuTuple<Ts...> data[nacc];
    template <int K> auto& get (int lane) noexcept { return amrex::get<K>(data[lane]); }
    template <int K> [[nodiscard]] auto const& get (int lane) const noexcept { return amrex::get<K>(data[lane]); }
#endif
};
 
template <int K, typename Tup>
AMREX_FORCE_INLINE
void lanes_set (Lanes<Tup>& acc, Tup const& v)
{
    for (int lane = 0; lane < nacc; ++lane) { acc.template get<K>(lane) = amrex::get<K>(v); }
    if constexpr (K+1 < GpuTupleSize<Tup>::value) { lanes_set<K+1>(acc, v); }
}
 
template <int K, typename Tup, typename P, typename... Ps>
AMREX_FORCE_INLINE
void lanes_update (Lanes<Tup>& acc, int lane, Tup const& s)
{
    P().local_update(acc.template get<K>(lane), amrex::get<K>(s));
    if constexpr (sizeof...(Ps) > 0) { lanes_update<K+1,Tup,Ps...>(acc, lane, s); }
}
 
template <int K, typename Tup, typename P, typename... Ps>
AMREX_FORCE_INLINE
void lanes_fold (Lanes<Tup> const& acc, Tup& r)
{
    for (int lane = 0; lane < nacc; ++lane) { P().local_update(amrex::get<K>(r), acc.template get<K>(lane)); }
    if constexpr (sizeof...(Ps) > 0) { lanes_fold<K+1,Tup,Ps...>(acc, r); }
}
 
// One row of a CPU reduction over the ops Ps: full blocks of nacc cells go
// into the lanes, the rest into rem. value(i) is the tuple for cell i.
template <typename... Ps, typename Tup, std::integral N, typename F>
AMREX_FORCE_INLINE
void reduce_row (N ilo, N ihi, Lanes<Tup>& acc, Tup& rem, F const& value)
{
    N i = ilo;
    AMREX_PRAGMA_NACC_NOVECTOR
    for (; i + (nacc-1) <= ihi; i += nacc) {
        if constexpr (vectorize_nacc_loop<Ps...>) {
            AMREX_PRAGMA_NACC_VECTORIZE
            for (int lane = 0; lane < nacc; ++lane) { lanes_update<0,Tup,Ps...>(acc, lane, value(static_cast<N>(i+lane))); }
        } else {
            AMREX_PRAGMA_NACC_SCALAR
            for (int lane = 0; lane < nacc; ++lane) { lanes_update<0,Tup,Ps...>(acc, lane, value(static_cast<N>(i+lane))); }
        }
    }
    AMREX_PRAGMA_NACC_NOVECTOR
    for (; i <= ihi; ++i) { for_each_local<0,Tup,Ps...>(rem, value(i)); }
}
}
#undef AMREX_PRAGMA_NACC_VECTORIZE
#undef AMREX_PRAGMA_NACC_NOVECTOR
#undef AMREX_PRAGMA_NACC_SCALAR
 
template <typename... Ps> class ReduceOps;
 
#ifdef AMREX_USE_GPU

template <typename... Ts>
class ReduceData
{
public:
    using Type = GpuTuple<Ts...>;
 
    template <typename... Ps>
    explicit ReduceData (ReduceOps<Ps...>& reduce_op)
        : m_max_blocks(Gpu::Device::maxBlocksPerLaunch()),
          m_host_tuple((Type*)(The_Pinned_Arena()->alloc(sizeof(Type)))),
          m_device_tuple((Type*)(The_Arena()->alloc((AMREX_GPU_MAX_STREAMS)
                                                    * m_max_blocks * sizeof(Type)))),
          m_fn_value([&reduce_op,this] () -> Type { return this->value(reduce_op); })
    {
        reduce_op.resetResultReadiness();
        static_assert(std::is_trivially_copyable<Type>(),
                      "ReduceData::Type must be trivially copyable");
        static_assert(std::is_trivially_destructible<Type>(),
                      "ReduceData::Type must be trivially destructible");
 
        new (m_host_tuple) Type();
        m_nblocks.fill(0);
    }
 
    ~ReduceData () {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            !m_used_external_stream || m_value_called,
            "ReduceData used on an external GPU stream must call value() before destruction.");
        The_Pinned_Arena()->free(m_host_tuple);
        The_Arena()->free(m_device_tuple);
    }
 
    ReduceData (ReduceData<Ts...> const&) = delete;
    ReduceData (ReduceData<Ts...> &&) = delete;
    void operator= (ReduceData<Ts...> const&) = delete;
    void operator= (ReduceData<Ts...> &&) = delete;
 
    Type value ()
    {
        Type r = m_fn_value();
        m_value_called = true;
        return r;
    }
 
    template <typename... Ps>
    Type value (ReduceOps<Ps...> & reduce_op)
    {
        Type r = reduce_op.value(*this);
        m_value_called = true;
        return r;
    }
 
    Type* devicePtr () { return m_device_tuple; }
    Type* devicePtr (gpuStream_t const& s) {
        return m_device_tuple+streamIndexChecked(s)*m_max_blocks;
    }
 
    Type* hostPtr () { return m_host_tuple; }
 
    GpuArray<int,AMREX_GPU_MAX_STREAMS>& nBlocks () { return m_nblocks; }
    int& nBlocks (gpuStream_t const& s) { return m_nblocks[streamIndexChecked(s)]; }
 
    int maxBlocks () const { return m_max_blocks; }
 
    int maxStreamIndex () const { return m_max_stream_index; }
    void updateMaxStreamIndex (gpuStream_t const& s) {
        m_max_stream_index = std::max(m_max_stream_index,streamIndexChecked(s));
    }
 
    void markValueCalled () noexcept { m_value_called = true; }
 
private:
    int streamIndexChecked (gpuStream_t const& s)
    {
        int const idx = Gpu::Device::streamIndex(s);
        m_used_external_stream = m_used_external_stream || Gpu::Device::usingExternalStream();
        if (idx == 0) {
            if (m_stream_index_zero_set) {
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_stream_index_zero == s,
                    "ReduceData cannot be reused across different external GPU streams "
                    "or between an external GPU stream and AMReX stream 0.");
            } else {
                m_stream_index_zero = s;
                m_stream_index_zero_set = true;
            }
        }
        return idx;
    }
 
    int m_max_blocks;
    int m_max_stream_index = 0;
    Type* m_host_tuple = nullptr;
    Type* m_device_tuple = nullptr;
    GpuArray<int,AMREX_GPU_MAX_STREAMS> m_nblocks;
    gpuStream_t m_stream_index_zero{};
    bool m_stream_index_zero_set = false;
    bool m_used_external_stream = false;
    bool m_value_called = false;
    std::function<Type()> m_fn_value;
};

namespace Reduce::detail {
 
    // call_f_intvect_box
 
    template <typename F, int dim>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    auto call_f_intvect_box (F const& f, IntVectND<dim> iv, IndexTypeND<dim>) noexcept ->
        decltype(amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv))
    {
        return amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv);
    }
 
    template <typename F, int dim>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    auto call_f_intvect_box (F const& f, IntVectND<dim> iv, IndexTypeND<dim> t) noexcept ->
        decltype(f(BoxND<dim>(iv, iv, t)))
    {
        return f(BoxND<dim>(iv, iv, t));
    }
 
    // call_f_intvect_n
    template <typename F, typename T, int dim>
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    auto call_f_intvect_n (F const& f, IntVectND<dim> iv, T n) noexcept ->
        decltype(amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv, n))
    {
        return amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv, n);
    }
 
    // mf_call_f
 
    struct iterate_box {};
    struct iterate_box_comp {};
 
    template <typename Tag, typename F, typename T, typename... Ps>
    requires (std::same_as<iterate_box,Tag>)
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void mf_call_f (F const& f, int ibox, int i, int j, int k, int, T& r) noexcept
    {
        auto const& pr = f(ibox,i,j,k);
        Reduce::detail::for_each_local<0, T, Ps...>(r, pr);
    }
 
    template <typename Tag, typename F, typename T, typename... Ps>
    requires (std::same_as<iterate_box_comp,Tag>)
    AMREX_GPU_DEVICE AMREX_FORCE_INLINE
    void mf_call_f (F const& f, int ibox, int i, int j, int k, int ncomp, T& r) noexcept
    {
        for (int n = 0; n < ncomp; ++n) {
            auto const& pr = f(ibox,i,j,k,n);
            Reduce::detail::for_each_local<0, T, Ps...>(r, pr);
        }
    }
}

template <typename... Ps>
class ReduceOps
{
public:

 
    // This is public for CUDA
    template <typename Tag, typename MF, typename D, typename F>
    void eval_mf (Tag, MF const& mf, IntVect const& nghost, int ncomp, D& reduce_data, F const& f)
    {
        using ReduceTuple = typename D::Type;
        const int nboxes = mf.local_size();
        if (nboxes > 0) {
            auto const& parforinfo = mf.getParForInfo(nghost);
            auto nblocks_per_box = parforinfo.getNBlocksPerBox(AMREX_GPU_MAX_THREADS);
            AMREX_ASSERT(Long(nblocks_per_box)*Long(nboxes) < Long(std::numeric_limits<int>::max()));
            const int nblocks = nblocks_per_box * nboxes;
            const BoxIndexer* dp_boxes = parforinfo.getBoxes();
 
            auto const& stream = Gpu::gpuStream();
            auto pdst = reduce_data.devicePtr(stream);
            int nblocks_ec = std::min(nblocks, reduce_data.maxBlocks());
            AMREX_ASSERT(Long(nblocks_ec)*2 <= Long(std::numeric_limits<int>::max()));
            int& nblocks_ref = reduce_data.nBlocks(stream);
            auto old_nblocks = static_cast<unsigned int>(nblocks_ref);
            nblocks_ref = amrex::max(nblocks_ref, nblocks_ec);
            reduce_data.updateMaxStreamIndex(stream);
 
#ifdef AMREX_USE_SYCL
            // device reduce needs local(i.e., shared) memory
            constexpr std::size_t shared_mem_bytes = sizeof(ReduceTuple)
                * std::max(Gpu::Device::warp_size, AMREX_GPU_MAX_THREADS/Gpu::Device::warp_size);
            amrex::launch<AMREX_GPU_MAX_THREADS>(nblocks_ec, shared_mem_bytes, stream,
                          [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept
            {
                Dim1 blockIdx {gh.blockIdx()};
                Dim1 threadIdx{gh.threadIdx()};
#else
            amrex::launch_global<AMREX_GPU_MAX_THREADS>
                <<<nblocks_ec, AMREX_GPU_MAX_THREADS, 0, stream>>>
                ([=] AMREX_GPU_DEVICE () noexcept
            {
#endif
                ReduceTuple r;
                Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(r);
                ReduceTuple& dst = pdst[blockIdx.x];
                if (threadIdx.x == 0 && blockIdx.x >= old_nblocks) {
                    dst = r;
                }
                for (int iblock = blockIdx.x; iblock < nblocks; iblock += nblocks_ec) {
                    int ibox = iblock / nblocks_per_box;
                    auto icell = std::uint64_t(iblock-ibox*nblocks_per_box)*AMREX_GPU_MAX_THREADS + threadIdx.x;
 
                    BoxIndexer const& indexer = dp_boxes[ibox];
                    if (icell < indexer.numPts()) {
                        auto [i, j, k] = indexer(icell);
                        Reduce::detail::mf_call_f<Tag, F, ReduceTuple, Ps...>
                            (f, ibox, i, j, k, ncomp, r);
                    }
                }
#ifdef AMREX_USE_SYCL
                Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r, gh);
#else
                Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r);
#endif
            });
        }
    }
 
    // This is public for CUDA
    template <typename Tag, int dim, typename D, typename F>
    void eval_box (Tag, BoxND<dim> const& box, int ncomp, D& reduce_data, F const& f)
    {
        if (box.isEmpty()) { return; }
        using ReduceTuple = typename D::Type;
        auto const& stream = Gpu::gpuStream();
        auto dp = reduce_data.devicePtr(stream);
        int& nblocks = reduce_data.nBlocks(stream);
        const BoxIndexerND<dim> indexer(box);
        IndexTypeND<dim> ixtype = box.ixType();
        constexpr int nitems_per_thread = 4;
        Long nblocks_ec = (box.numPts() + nitems_per_thread*AMREX_GPU_MAX_THREADS-1)
            / (nitems_per_thread*AMREX_GPU_MAX_THREADS);
        nblocks_ec = std::min<Long>(nblocks_ec, reduce_data.maxBlocks());
        reduce_data.updateMaxStreamIndex(stream);
#ifdef AMREX_USE_SYCL
        // device reduce needs local(i.e., shared) memory
        constexpr std::size_t shared_mem_bytes = sizeof(ReduceTuple)
            * std::max(Gpu::Device::warp_size, AMREX_GPU_MAX_THREADS/Gpu::Device::warp_size);
        amrex::launch<AMREX_GPU_MAX_THREADS>(nblocks_ec, shared_mem_bytes, stream,
        [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept
        {
            Dim1 blockIdx {gh.blockIdx()};
            Dim1 threadIdx{gh.threadIdx()};
            Dim1 gridDim  {gh.gridDim()};
#else
        amrex::launch<AMREX_GPU_MAX_THREADS>(nblocks_ec, 0, stream,
        [=] AMREX_GPU_DEVICE () noexcept
        {
#endif
            ReduceTuple r;
            Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(r);
            ReduceTuple& dst = *(dp+blockIdx.x);
            if (threadIdx.x == 0 && static_cast<int>(blockIdx.x) >= nblocks) {
                dst = r;
            }
            for (std::uint64_t icell = std::uint64_t(AMREX_GPU_MAX_THREADS)*blockIdx.x+threadIdx.x,
                 stride = std::uint64_t(AMREX_GPU_MAX_THREADS)*gridDim.x;
                 icell < indexer.numPts();
                 icell += stride)
            {
                auto iv = indexer.intVect(icell);
                amrex::ignore_unused(f,ncomp,ixtype); // work around first-capture
                if constexpr (std::is_same_v<Reduce::detail::iterate_box,Tag>) {
                    auto pr = Reduce::detail::call_f_intvect_box(f, iv, ixtype);
                    Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, pr);
                } else {
                    for (int n = 0; n < ncomp; ++n) {
                        auto pr = Reduce::detail::call_f_intvect_n(f, iv, n);
                        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, pr);
                    }
                }
            }
#ifdef AMREX_USE_SYCL
            Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r, gh);
#else
            Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r);
#endif
        });
        nblocks = std::max(nblocks, static_cast<int>(nblocks_ec));
    }

 
    template <FabArrayType MF, typename D, typename F>
#ifndef AMREX_USE_CUDA
    requires (IsCallable<F, int, int, int, int>::value)
#endif
    void eval (MF const& mf, IntVect const& nghost, D& reduce_data, F&& f)
    {
        using ReduceTuple = typename D::Type;
        const int nboxes = mf.local_size();
        if (nboxes == 0) {
            return;
        } else if (!mf.isFusingCandidate()) {
            for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
                Box const& b = amrex::grow(mfi.validbox(), nghost);
                const int li = mfi.LocalIndex();
                this->eval(b, reduce_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> ReduceTuple
                {
                    return f(li, i, j, k);
                });
            }
        } else {
            eval_mf(Reduce::detail::iterate_box{},
                    mf, nghost, 0, reduce_data, std::forward<F>(f));
        }
    }
 
    template <FabArrayType MF, typename D, typename F>
#ifndef AMREX_USE_CUDA
    requires (IsCallable<F, int, int, int, int, int>::value)
#endif
    void eval (MF const& mf, IntVect const& nghost, int ncomp, D& reduce_data, F&& f)
    {
        using ReduceTuple = typename D::Type;
 
        const int nboxes = mf.local_size();
 
        if (nboxes == 0) {
            return;
        } else if (!mf.isFusingCandidate()) {
            for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
                Box const& b = amrex::grow(mfi.validbox(), nghost);
                const int li = mfi.LocalIndex();
                this->eval(b, ncomp, reduce_data,
                [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept -> ReduceTuple
                {
                    return f(li, i, j, k, n);
                });
            }
        } else if (ncomp == 1) {
            // Without the loop over components, the kernel reads memory faster.
            eval_mf(Reduce::detail::iterate_box{}, mf, nghost, 0, reduce_data,
            [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> ReduceTuple
            {
                return f(box_no, i, j, k, 0);
            });
        } else {
            eval_mf(Reduce::detail::iterate_box_comp{},
                    mf, nghost, ncomp, reduce_data, std::forward<F>(f));
        }
    }
 
    template <typename D, typename F, int dim>
    void eval (BoxND<dim> const& box, D & reduce_data, F const& f)
    {
        eval_box(Reduce::detail::iterate_box{}, box, 0, reduce_data, f);
    }
 
    template <std::integral N, typename D, typename F, int dim>
    void eval (BoxND<dim> const& box, N ncomp, D & reduce_data, F const& f)
    {
        if (ncomp == 1) {
            // Without the loop over components, the kernel reads memory faster.
            eval_box(Reduce::detail::iterate_box{}, box, 0, reduce_data,
            [=] AMREX_GPU_DEVICE (IntVectND<dim> const& iv) noexcept
            {
                return Reduce::detail::call_f_intvect_n(f, iv, 0);
            });
        } else {
            eval_box(Reduce::detail::iterate_box_comp{}, box, ncomp, reduce_data, f);
        }
    }
 
    template <std::integral N, typename D, typename F>
    void eval (N n, D & reduce_data, F const& f)
    {
        if (n <= 0) { return; }
        using ReduceTuple = typename D::Type;
        auto const& stream = Gpu::gpuStream();
        auto dp = reduce_data.devicePtr(stream);
        int& nblocks = reduce_data.nBlocks(stream);
        constexpr int nitems_per_thread = 4;
        Long nblocks_ec_long = (Long(n) + nitems_per_thread*AMREX_GPU_MAX_THREADS-1)
            / (nitems_per_thread*AMREX_GPU_MAX_THREADS);
        int nblocks_ec = static_cast<int>(std::min<Long>(nblocks_ec_long,
                                                         reduce_data.maxBlocks()));
        reduce_data.updateMaxStreamIndex(stream);
#ifdef AMREX_USE_SYCL
        // device reduce needs local(i.e., shared) memory
        constexpr std::size_t shared_mem_bytes = sizeof(ReduceTuple)
            * std::max(Gpu::Device::warp_size, AMREX_GPU_MAX_THREADS/Gpu::Device::warp_size);
        amrex::launch<AMREX_GPU_MAX_THREADS>(nblocks_ec, shared_mem_bytes, stream,
        [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept
        {
            Dim1 blockIdx {gh.blockIdx()};
            Dim1 threadIdx{gh.threadIdx()};
            Dim1 gridDim  {gh.gridDim()};
#else
        amrex::launch<AMREX_GPU_MAX_THREADS>(nblocks_ec, 0, stream,
        [=] AMREX_GPU_DEVICE () noexcept
        {
#endif
            ReduceTuple r;
            Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(r);
            ReduceTuple& dst = *(dp+blockIdx.x);
            if (threadIdx.x == 0 && static_cast<int>(blockIdx.x) >= nblocks) {
                dst = r;
            }
            for (Long i = Long(AMREX_GPU_MAX_THREADS)*blockIdx.x+threadIdx.x,
                 stride = Long(AMREX_GPU_MAX_THREADS)*gridDim.x;
                 i < Long(n);
                 i += stride)
            {
                auto pr = f(N(i));
                Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r,pr);
            }
#ifdef AMREX_USE_SYCL
            Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r, gh);
#else
            Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r);
#endif
        });
        nblocks = amrex::max(nblocks, nblocks_ec);
    }
 
    template <typename D>
    typename D::Type value (D & reduce_data)
    {
        auto hp = reduce_data.hostPtr();
 
        if (m_result_is_ready) {
            reduce_data.markValueCalled();
            return *hp;
        }
 
        using ReduceTuple = typename D::Type;
        auto const& stream = Gpu::gpuStream();
        auto dp = reduce_data.devicePtr();
        auto const& nblocks = reduce_data.nBlocks();
#if defined(AMREX_USE_SYCL)
        if (reduce_data.maxStreamIndex() == 0 && nblocks[0] <= 4096) {
            const int N = nblocks[0];
            if (N == 0) {
                Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(*hp);
            } else {
                Gpu::PinnedVector<ReduceTuple> tmp(N);
                Gpu::dtoh_memcpy_async(tmp.data(), dp, sizeof(ReduceTuple)*N);
                Gpu::streamSynchronize();
                for (int i = 1; i < N; ++i) {
                    Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(tmp[0], tmp[i]);
                }
                *hp = tmp[0];
            }
        } else
#endif
        {
            int maxblocks = reduce_data.maxBlocks();
#ifdef AMREX_USE_SYCL
            // device reduce needs local(i.e., shared) memory
            constexpr std::size_t shared_mem_bytes = sizeof(ReduceTuple)
                * std::max(Gpu::Device::warp_size, AMREX_GPU_MAX_THREADS/Gpu::Device::warp_size);
#ifndef AMREX_NO_SYCL_REDUCE_WORKAROUND
            // xxxxx SYCL todo: reduce bug workaround
            Gpu::DeviceVector<ReduceTuple> dtmp(1);
            auto presult = dtmp.data();
#else
            auto presult = hp;
#endif
            amrex::launch<AMREX_GPU_MAX_THREADS>(1, shared_mem_bytes, stream,
            [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept
            {
                ReduceTuple r;
                Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(r);
                ReduceTuple dst = r;
                for (int istream = 0, nstreams = nblocks.size(); istream < nstreams; ++istream) {
                    auto dp_stream = dp+istream*maxblocks;
                    for (int i = gh.item->get_global_id(0), stride = gh.item->get_global_range(0);
                         i < nblocks[istream]; i += stride) {
                        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, dp_stream[i]);
                    }
                }
                Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r, gh);
                if (gh.threadIdx() == 0) { *presult = dst; }
            });
#ifndef AMREX_NO_SYCL_REDUCE_WORKAROUND
            Gpu::dtoh_memcpy_async(hp, dtmp.data(), sizeof(ReduceTuple));
#endif
#else
            amrex::launch<AMREX_GPU_MAX_THREADS>(1, 0, stream,
            [=] AMREX_GPU_DEVICE () noexcept
            {
                ReduceTuple r;
                Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(r);
                ReduceTuple dst = r;
                for (int istream = 0, nstreams = nblocks.size(); istream < nstreams; ++istream) {
                    auto dp_stream = dp+istream*maxblocks;
                    for (int i = AMREX_GPU_MAX_THREADS*blockIdx.x+threadIdx.x, stride = AMREX_GPU_MAX_THREADS*gridDim.x;
                         i < nblocks[istream]; i += stride) {
                        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, dp_stream[i]);
                    }
                }
                Reduce::detail::for_each_parallel<0, ReduceTuple, Ps...>(dst, r);
                if (threadIdx.x == 0) { *hp = dst; }
            });
#endif
            Gpu::streamSynchronize();
        }
 
        m_result_is_ready = true;
        reduce_data.markValueCalled();
        return *hp;
    }
 
private:
    template <typename... T> friend class ReduceData;
    bool m_result_is_ready = false;
    void resetResultReadiness () { m_result_is_ready = false; }
};
 
namespace Reduce {
 
template <typename T, std::integral N>
T Sum (N n, T const* v, T init_val)
{
    ReduceOps<ReduceOpSum> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {v[i]}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return amrex::get<0>(hv) + init_val;
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Sum (N n, F const& f, T init_val)
{
    ReduceOps<ReduceOpSum> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {f(i)}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return amrex::get<0>(hv) + init_val;
}
 
template <typename T, std::integral N>
T Min (N n, T const* v, T init_val)
{
    ReduceOps<ReduceOpMin> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {v[i]}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return std::min(amrex::get<0>(hv),init_val);
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Min (N n, F const& f, T init_val)
{
    ReduceOps<ReduceOpMin> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {f(i)}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return std::min(amrex::get<0>(hv),init_val);
}
 
template <typename T, std::integral N>
T Max (N n, T const* v, T init_val)
{
    ReduceOps<ReduceOpMax> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {v[i]}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return std::max(amrex::get<0>(hv),init_val);
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Max (N n, F const& f, T init_val)
{
    ReduceOps<ReduceOpMax> reduce_op;
    ReduceData<T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple { return {f(i)}; });
    ReduceTuple hv = reduce_data.value(reduce_op);
    return std::max(amrex::get<0>(hv),init_val);
}
 
template <typename T, std::integral N>
std::pair<T,T> MinMax (N n, T const* v)
{
    ReduceOps<ReduceOpMin,ReduceOpMax> reduce_op;
    ReduceData<T,T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple {
            return {v[i],v[i]};
        });
    auto hv = reduce_data.value(reduce_op);
    return std::make_pair(amrex::get<0>(hv), amrex::get<1>(hv));
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
std::pair<T,T> MinMax (N n, F const& f)
{
    ReduceOps<ReduceOpMin,ReduceOpMax> reduce_op;
    ReduceData<T,T> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;
    reduce_op.eval(n, reduce_data, [=] AMREX_GPU_DEVICE (N i) -> ReduceTuple {
            T tmp = f(i);
            return {tmp,tmp};
        });
    auto hv = reduce_data.value(reduce_op);
    return std::make_pair(amrex::get<0>(hv), amrex::get<1>(hv));
}
 
template <typename T, std::integral N, typename P>
bool AnyOf (N n, T const* v, P const& pred)
{
    Gpu::LaunchSafeGuard lsg(true);
    Gpu::DeviceScalar<int> ds(0);
    int* dp = ds.dataPtr();
    auto ec = Gpu::ExecutionConfig(n);
    ec.numBlocks.x = std::min(ec.numBlocks.x, Gpu::Device::maxBlocksPerLaunch());
 
#ifdef AMREX_USE_SYCL
    const int num_ints = std::max(Gpu::Device::warp_size, int(ec.numThreads.x)/Gpu::Device::warp_size) + 1;
    const std::size_t shared_mem_bytes = num_ints*sizeof(int);
    amrex::launch<AMREX_GPU_MAX_THREADS>(ec.numBlocks.x, shared_mem_bytes, Gpu::gpuStream(),
    [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept {
        int* has_any = &(static_cast<int*>(gh.sharedMemory())[num_ints-1]);
        if (gh.threadIdx() == 0) { *has_any = *dp; }
        gh.syncThreads();
 
        if (!(*has_any))
        {
            int r = false;
            for (Long i = Long(AMREX_GPU_MAX_THREADS)*gh.blockIdx()+gh.threadIdx(),
                 stride = Long(AMREX_GPU_MAX_THREADS)*gh.gridDim();
                 i < Long(n) && !r; i += stride)
            {
                r = pred(v[i]) ? 1 : 0;
            }
 
            r = Gpu::blockReduce<Gpu::Device::warp_size>
                (r, Gpu::warpReduce<Gpu::Device::warp_size,int,amrex::Plus<int> >(), 0, gh);
            if (gh.threadIdx() == 0 && r) { *dp = 1; }
        }
    });
#else
    amrex::launch<AMREX_GPU_MAX_THREADS>(ec.numBlocks.x, 0, Gpu::gpuStream(),
    [=] AMREX_GPU_DEVICE () noexcept {
        __shared__ int has_any;
        if (threadIdx.x == 0) { has_any = *dp; }
        __syncthreads();
 
        if (!has_any)
        {
            int r = false;
            for (Long i = Long(AMREX_GPU_MAX_THREADS)*blockIdx.x+threadIdx.x,
                 stride = Long(AMREX_GPU_MAX_THREADS)*gridDim.x;
                 i < Long(n) && !r; i += stride)
            {
                r = pred(v[i]) ? 1 : 0;
            }
            r = Gpu::blockReduce<Gpu::Device::warp_size>
                (r, Gpu::warpReduce<Gpu::Device::warp_size,int,amrex::Plus<int> >(), 0);
            if (threadIdx.x == 0 && r) *dp = 1;
        }
    });
#endif
    return ds.dataValue();
}
 
template <typename P, int dim>
bool AnyOf (BoxND<dim> const& box, P const& pred)
{
    Gpu::LaunchSafeGuard lsg(true);
    Gpu::DeviceScalar<int> ds(0);
    int* dp = ds.dataPtr();
    const BoxIndexerND<dim> indexer(box);
    auto ec = Gpu::ExecutionConfig(box.numPts());
    ec.numBlocks.x = std::min(ec.numBlocks.x, Gpu::Device::maxBlocksPerLaunch());
 
#ifdef AMREX_USE_SYCL
    const int num_ints = std::max(Gpu::Device::warp_size, int(ec.numThreads.x)/Gpu::Device::warp_size) + 1;
    const std::size_t shared_mem_bytes = num_ints*sizeof(int);
    amrex::launch<AMREX_GPU_MAX_THREADS>(ec.numBlocks.x, shared_mem_bytes, Gpu::gpuStream(),
    [=] AMREX_GPU_DEVICE (Gpu::Handler const& gh) noexcept {
        int* has_any = &(static_cast<int*>(gh.sharedMemory())[num_ints-1]);
        if (gh.threadIdx() == 0) { *has_any = *dp; }
        gh.syncThreads();
 
        if (!(*has_any))
        {
            int r = false;
            for (std::uint64_t icell = std::uint64_t(AMREX_GPU_MAX_THREADS)*gh.blockIdx()+gh.threadIdx(),
                 stride = std::uint64_t(AMREX_GPU_MAX_THREADS)*gh.gridDim();
                 icell < indexer.numPts() && !r;
                 icell += stride)
            {
                auto iv = indexer.intVect(icell);
                r = amrex::detail::call_f_intvect(pred, iv) ? 1 : 0;
            }
            r = Gpu::blockReduce<Gpu::Device::warp_size>
                (r, Gpu::warpReduce<Gpu::Device::warp_size,int,amrex::Plus<int> >(), 0, gh);
            if (gh.threadIdx() == 0 && r) { *dp = 1; }
        }
    });
#else
    AMREX_LAUNCH_KERNEL(AMREX_GPU_MAX_THREADS, ec.numBlocks, ec.numThreads, 0,
                        Gpu::gpuStream(),
    [=] AMREX_GPU_DEVICE () noexcept {
        __shared__ int has_any;
        if (threadIdx.x == 0) { has_any = *dp; }
        __syncthreads();
 
        if (!has_any)
        {
            int r = false;
            for (std::uint64_t icell = std::uint64_t(AMREX_GPU_MAX_THREADS)*blockIdx.x+threadIdx.x,
                 stride = std::uint64_t(AMREX_GPU_MAX_THREADS)*gridDim.x;
                 icell < indexer.numPts() && !r;
                 icell += stride)
            {
                auto iv = indexer.intVect(icell);
                r = amrex::detail::call_f_intvect(pred, iv) ? 1 : 0;
            }
            r = Gpu::blockReduce<Gpu::Device::warp_size>
                (r, Gpu::warpReduce<Gpu::Device::warp_size,int,amrex::Plus<int> >(), 0);
            if (threadIdx.x == 0 && r) *dp = 1;
        }
    });
#endif
    return ds.dataValue();
}
 
}
 
#else
 
template <typename... Ts>
class ReduceData
{
public:
    using Type = GpuTuple<Ts...>;
 
    template <typename... Ps>
    explicit ReduceData (ReduceOps<Ps...>& reduce_op)
        : m_tuple(OpenMP::in_parallel() ? 1 : OpenMP::get_max_threads()),
          m_fn_value([&reduce_op,this] () -> Type { return this->value(reduce_op); })
    {
        reduce_op.resetResultReadiness();
        for (auto& t : m_tuple) {
            Reduce::detail::for_each_init<0, Type, Ps...>(t);
        }
    }
 
    ~ReduceData () = default;
    ReduceData (ReduceData<Ts...> const&) = delete;
    ReduceData (ReduceData<Ts...> &&) = delete;
    void operator= (ReduceData<Ts...> const&) = delete;
    void operator= (ReduceData<Ts...> &&) = delete;
 
    Type value () { return m_fn_value(); }
 
    template <typename... Ps>
    Type value (ReduceOps<Ps...>& reduce_op)
    {
        return reduce_op.value(*this);
    }
 
    Vector<Type>& reference () { return m_tuple; }
 
    Type& reference (int tid)
    {
        if (m_tuple.size() == 1) {
            // No OpenMP or already inside OpenMP parallel when reduce_data is constructed
            return m_tuple[0];
        } else {
            return m_tuple[tid];
        }
    }
 
private:
    Vector<Type> m_tuple;
    std::function<Type()> m_fn_value;
};
 
namespace Reduce::detail {
 
    // call_f_intvect
 
    template <typename F, int dim>
    AMREX_FORCE_INLINE
    auto call_f_intvect (F const& f, IntVectND<dim> iv) noexcept ->
        decltype(amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv))
    {
        return amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv);
    }
 
    // call_f_intvect_n
 
    template <typename F, typename T, int dim>
    AMREX_FORCE_INLINE
    auto call_f_intvect_n (F const& f, IntVectND<dim> iv, T n) noexcept ->
        decltype(amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv, n))
    {
        return amrex::detail::call_f_intvect_inner(std::make_index_sequence<dim>(), f, iv, n);
    }
}
 
template <typename... Ps>
class ReduceOps
{
private:
 
    // call_f_box
 
    template <typename D, typename F, int dim>
    requires (std::same_as<std::decay_t<decltype(
                  Reduce::detail::call_f_intvect(std::declval<F const&>(), IntVectND<dim>()))>,
              typename D::Type>)
    AMREX_FORCE_INLINE
    static void call_f_box (BoxND<dim> const& box, typename D::Type & r, F const& f) noexcept
    {
        using ReduceTuple = typename D::Type;
        ReduceTuple rem;
        Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rem);
        Reduce::detail::Lanes<ReduceTuple> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        int const ilo = box.smallEnd(0);
        int const ihi = box.bigEnd(0);
        BoxND<dim> rows = box;
        rows.setBig(0, ilo);
        For(rows, [&] (IntVectND<dim> iv) {
            Reduce::detail::reduce_row<Ps...>(ilo, ihi, acc, rem, [&] (int i) {
                iv[0] = i;
                return Reduce::detail::call_f_intvect(f, iv);
            });
        });
        Reduce::detail::lanes_fold<0, ReduceTuple, Ps...>(acc, r);
        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, rem);
    }
 
    template <typename D, typename F, int dim>
    requires (std::same_as<std::decay_t<decltype(
                  std::declval<F const&>()(std::declval<BoxND<dim> const&>()))>,
              typename D::Type>)
    AMREX_FORCE_INLINE
    static void call_f_box (BoxND<dim> const& box, typename D::Type & r, F const& f) noexcept
    {
        using ReduceTuple = typename D::Type;
        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(r, f(box));
    }
 
public:
 
    template <FabArrayType MF, typename D, typename F>
    requires (IsCallable<F, int, int, int, int>::value)
    void eval (MF const& mf, IntVect const& nghost, D & reduce_data, F const& f)
    {
        using ReduceTuple = typename D::Type;
#ifdef AMREX_USE_OMP
#pragma omp parallel if (!system::regtest_reduction)
#endif
        {
            // The partial accumulators break the dependency chain so that
            // the i loop can vectorize. They live across the MFIter loop so
            // that they stay in vector registers.
            ReduceTuple rem;
            Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rem);
            Reduce::detail::Lanes<ReduceTuple> acc;
            Reduce::detail::lanes_set<0>(acc, rem);
            for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
                Box const& b = mfi.growntilebox(nghost);
                const int li = mfi.LocalIndex();
                const auto lo = amrex::lbound(b);
                const auto hi = amrex::ubound(b);
                for (int k = lo.z; k <= hi.z; ++k) {
                for (int j = lo.y; j <= hi.y; ++j) {
                    Reduce::detail::reduce_row<Ps...>(lo.x, hi.x, acc, rem,
                        [&] (int i) { return f(li,i,j,k); });
                }}
            }
            auto& rr = reduce_data.reference(OpenMP::get_thread_num());
            Reduce::detail::lanes_fold<0, ReduceTuple, Ps...>(acc, rr);
            Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(rr, rem);
        }
    }
 
    template <FabArrayType MF, typename D, typename F>
    requires (IsCallable<F, int, int, int, int, int>::value)
    void eval (MF const& mf, IntVect const& nghost, int ncomp, D & reduce_data, F const& f)
    {
        using ReduceTuple = typename D::Type;
#ifdef AMREX_USE_OMP
#pragma omp parallel if (!system::regtest_reduction)
#endif
        {
            ReduceTuple rem;
            Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rem);
            Reduce::detail::Lanes<ReduceTuple> acc;
            Reduce::detail::lanes_set<0>(acc, rem);
            for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
                Box const& b = mfi.growntilebox(nghost);
                const int li = mfi.LocalIndex();
                const auto lo = amrex::lbound(b);
                const auto hi = amrex::ubound(b);
                for (int n = 0; n < ncomp; ++n) {
                for (int k = lo.z; k <= hi.z; ++k) {
                for (int j = lo.y; j <= hi.y; ++j) {
                    Reduce::detail::reduce_row<Ps...>(lo.x, hi.x, acc, rem,
                        [&] (int i) { return f(li,i,j,k,n); });
                }}}
            }
            auto& rr = reduce_data.reference(OpenMP::get_thread_num());
            Reduce::detail::lanes_fold<0, ReduceTuple, Ps...>(acc, rr);
            Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(rr, rem);
        }
    }
 
    template <typename D, typename F, int dim>
    void eval (BoxND<dim> const& box, D & reduce_data, F&& f)
    {
        using ReduceTuple = typename D::Type;
        ReduceTuple rr;
        Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rr);
        call_f_box<D>(box, rr, std::forward<F>(f));
        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(
             reduce_data.reference(OpenMP::get_thread_num()), rr);
    }
 
    template <std::integral N, typename D, typename F, int dim>
    void eval (BoxND<dim> const& box, N ncomp, D & reduce_data, F const& f)
    {
        using ReduceTuple = typename D::Type;
        ReduceTuple rem;
        Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rem);
        Reduce::detail::Lanes<ReduceTuple> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        int const ilo = box.smallEnd(0);
        int const ihi = box.bigEnd(0);
        BoxND<dim> rows = box;
        rows.setBig(0, ilo);
        For(rows, ncomp, [&] (IntVectND<dim> iv, int n) {
            Reduce::detail::reduce_row<Ps...>(ilo, ihi, acc, rem, [&] (int i) {
                iv[0] = i;
                return Reduce::detail::call_f_intvect_n(f, iv, n);
            });
        });
        auto& rr = reduce_data.reference(OpenMP::get_thread_num());
        Reduce::detail::lanes_fold<0, ReduceTuple, Ps...>(acc, rr);
        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(rr, rem);
    }
 
    template <std::integral N, typename D, typename F>
    void eval (N n, D & reduce_data, F const& f)
    {
        using ReduceTuple = typename D::Type;
        ReduceTuple rem;
        Reduce::detail::for_each_init<0, ReduceTuple, Ps...>(rem);
        Reduce::detail::Lanes<ReduceTuple> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        if (n > 0) {
            Reduce::detail::reduce_row<Ps...>(N(0), static_cast<N>(n-1), acc, rem,
                [&] (N i) { return f(i); });
        }
        auto& rr = reduce_data.reference(OpenMP::get_thread_num());
        Reduce::detail::lanes_fold<0, ReduceTuple, Ps...>(acc, rr);
        Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(rr, rem);
    }
 
    template <typename D>
    typename D::Type value (D & reduce_data)
    {
        auto& rrv = reduce_data.reference();
        if (! m_result_is_ready) {
            using ReduceTuple = typename D::Type;
            if (rrv.size() > 1) {
                for (int i = 1, N = rrv.size(); i < N; ++i) {
                    Reduce::detail::for_each_local<0, ReduceTuple, Ps...>(rrv[0], rrv[i]);
                }
            }
            m_result_is_ready = true;
        }
        return rrv[0];
    }
 
private:
    template <typename... T> friend class ReduceData;
    bool m_result_is_ready = false;
    void resetResultReadiness () { m_result_is_ready = false; }
};
 
namespace Reduce {
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Sum (N n, F const& f, T init_val)
{
    auto const cn = Reduce::detail::chunks_1d(n);
    N const chunk = cn.first;
    N const nchunks = cn.second;
    T r = init_val;
#ifdef AMREX_USE_OMP
#pragma omp parallel for reduction(+:r)
#endif
    for (N ic = 0; ic < nchunks; ++ic) {
        N const ib = ic*chunk;
        N const ie = (n - ib < chunk) ? n : ib + chunk;
        using Tup = GpuTuple<T>;
        Tup rem{};
        Reduce::detail::Lanes<Tup> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        Reduce::detail::reduce_row<ReduceOpSum>(ib, static_cast<N>(ie-1), acc, rem,
            [&] (N i) { return Tup{f(i)}; });
        Reduce::detail::lanes_fold<0,Tup,ReduceOpSum>(acc, rem);
        r += amrex::get<0>(rem);
    }
    return r;
}
 
template <typename T, std::integral N>
T Sum (N n, T const* v, T init_val)
{
    return Sum(n, [=] (N i) -> T { return v[i]; }, init_val);
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Min (N n, F const& f, T init_val)
{
    auto const cn = Reduce::detail::chunks_1d(n);
    N const chunk = cn.first;
    N const nchunks = cn.second;
    T r = init_val;
#ifdef AMREX_USE_OMP
#pragma omp parallel for reduction(min:r)
#endif
    for (N ic = 0; ic < nchunks; ++ic) {
        N const ib = ic*chunk;
        N const ie = (n - ib < chunk) ? n : ib + chunk;
        using Tup = GpuTuple<T>;
        Tup rem{init_val};
        Reduce::detail::Lanes<Tup> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        Reduce::detail::reduce_row<ReduceOpMin>(ib, static_cast<N>(ie-1), acc, rem,
            [&] (N i) { return Tup{f(i)}; });
        Reduce::detail::lanes_fold<0,Tup,ReduceOpMin>(acc, rem);
        r = std::min(r, amrex::get<0>(rem));
    }
    return r;
}
 
template <typename T, std::integral N>
T Min (N n, T const* v, T init_val)
{
    return Reduce::Min(n, [=] (N i) -> T { return v[i]; }, init_val);
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
T Max (N n, F const& f, T init_val)
{
    auto const cn = Reduce::detail::chunks_1d(n);
    N const chunk = cn.first;
    N const nchunks = cn.second;
    T r = init_val;
#ifdef AMREX_USE_OMP
#pragma omp parallel for reduction(max:r)
#endif
    for (N ic = 0; ic < nchunks; ++ic) {
        N const ib = ic*chunk;
        N const ie = (n - ib < chunk) ? n : ib + chunk;
        using Tup = GpuTuple<T>;
        Tup rem{init_val};
        Reduce::detail::Lanes<Tup> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        Reduce::detail::reduce_row<ReduceOpMax>(ib, static_cast<N>(ie-1), acc, rem,
            [&] (N i) { return Tup{f(i)}; });
        Reduce::detail::lanes_fold<0,Tup,ReduceOpMax>(acc, rem);
        r = std::max(r, amrex::get<0>(rem));
    }
    return r;
}
 
template <typename T, std::integral N>
T Max (N n, T const* v, T init_val)
{
    return Reduce::Max(n, [=] (N i) -> T { return v[i]; }, init_val);
}
 
template <typename T, std::integral N, typename F>
requires (!std::same_as<T*,std::decay_t<F>>)
std::pair<T,T> MinMax (N n, F const& f)
{
    auto const cn = Reduce::detail::chunks_1d(n);
    N const chunk = cn.first;
    N const nchunks = cn.second;
    T r_min = std::numeric_limits<T>::max();
    T r_max = std::numeric_limits<T>::lowest();
#ifdef AMREX_USE_OMP
#pragma omp parallel for reduction(min:r_min) reduction(max:r_max)
#endif
    for (N ic = 0; ic < nchunks; ++ic) {
        N const ib = ic*chunk;
        N const ie = (n - ib < chunk) ? n : ib + chunk;
        using Tup = GpuTuple<T,T>;
        Tup rem{std::numeric_limits<T>::max(), std::numeric_limits<T>::lowest()};
        Reduce::detail::Lanes<Tup> acc;
        Reduce::detail::lanes_set<0>(acc, rem);
        Reduce::detail::reduce_row<ReduceOpMin,ReduceOpMax>(ib, static_cast<N>(ie-1), acc, rem,
            [&] (N i) { T const tmp = f(i); return Tup{tmp,tmp}; });
        Reduce::detail::lanes_fold<0,Tup,ReduceOpMin,ReduceOpMax>(acc, rem);
        r_min = std::min(r_min, amrex::get<0>(rem));
        r_max = std::max(r_max, amrex::get<1>(rem));
    }
    return std::make_pair(r_min,r_max);
}
 
template <typename T, std::integral N>
std::pair<T,T> MinMax (N n, T const* v)
{
    return Reduce::MinMax<T>(n, [=] (N i) -> T { return v[i]; });
}
 
template <typename T, std::integral N, typename P>
bool AnyOf (N n, T const* v, P const& pred)
{
    return std::any_of(v, v+n, pred);
}
 
template <typename P, int dim>
bool AnyOf (BoxND<dim> const& box, P const& pred)
{
    for (auto iv : box.iterator()) { // NOLINT(readability-use-anyofallof)
        if (Reduce::detail::call_f_intvect(pred, iv)) { return true; }
    }
    return false;
}
 
}
 
#endif

template <typename... Ts, typename... Ps>
AMREX_GPU_HOST_DEVICE
constexpr GpuTuple<Ts...>
IdentityTuple (GpuTuple<Ts...>, ReduceOps<Ps...>) noexcept
{
    GpuTuple<Ts...> r{};
    Reduce::detail::for_each_init<0, decltype(r), Ps...>(r);
    return r;
}

template <typename... Ts, typename... Ps>
AMREX_GPU_HOST_DEVICE
constexpr GpuTuple<Ts...>
IdentityTuple (GpuTuple<Ts...>, TypeList<Ps...>) noexcept
{
    GpuTuple<Ts...> r{};
    Reduce::detail::for_each_init<0, decltype(r), Ps...>(r);
    return r;
}

template <typename Ops, typename Ts>
class ReducerImpl;
 
template <typename... Ops, typename... Ts>
class ReducerImpl<TypeList<Ops...>, TypeList<Ts...>>
{
public:
    static_assert(sizeof...(Ops) > 0);
    static_assert(sizeof...(Ts) > 0);
    static_assert(sizeof...(Ops) == sizeof...(Ts));
 
    ReducerImpl ()
        : m_reduce_data(m_reduce_op)
        {}
 
protected:
    using Result_t = GpuTuple<Ts...>;
    ReduceOps<Ops...> m_reduce_op;
    ReduceData<Ts...> m_reduce_data;
};

template <typename Ops, typename Ts>
class Reducer
    : public ReducerImpl<ToTypeList_t<Ops>, ToTypeList_t<Ts>>
{
    using Base = ReducerImpl<ToTypeList_t<Ops>, ToTypeList_t<Ts>>;
public:
 
    using Result_t = typename Base::Result_t; 
    static constexpr int size = GpuTupleSize<Result_t>::value; 
 
    Reducer () = default;
    ~Reducer () = default;

    Reducer (Reducer const&) = delete;
    Reducer (Reducer &&) = delete;
    void operator= (Reducer const&) = delete;
    void operator= (Reducer &&) = delete;

    template <typename F, int dim>
    requires (IsCallable<F, int, int, int>::value ||
              IsCallable<F, IntVectND<dim>>::value)
    void eval (BoxND<dim> const& box, F&& f)
    {
        this->m_reduce_op.eval(box, this->m_reduce_data, std::forward<F>(f));
    }

    template <typename F, int dim>
    requires (IsCallable<F, int, int, int, int>::value ||
              IsCallable<F, IntVectND<dim>, int>::value)
    void eval (BoxND<dim> const& box, int ncomp, F&& f)
    {
        this->m_reduce_op.eval(box, ncomp, this->m_reduce_data, std::forward<F>(f));
    }

    template <FabArrayType MF, typename F>
    requires (IsCallable<F, int, int, int, int>::value)
    void eval (MF const& mf, IntVect const& nghost, F && f)
    {
        this->m_reduce_op.eval(mf, nghost, this->m_reduce_data, std::forward<F>(f));
    }

    template <FabArrayType MF, typename F>
    requires (IsCallable<F, int, int, int, int, int>::value)
    void eval (MF const& mf, IntVect const& nghost, int ncomp, F && f)
    {
        this->m_reduce_op.eval(mf, nghost, ncomp, this->m_reduce_data, std::forward<F>(f));
    }

    template <typename N, typename F>
    requires (IsCallable<F, N>::value)
    void eval (N n, F && f)
    {
        this->m_reduce_op.eval(n, this->m_reduce_data, std::forward<F>(f));
    }

    [[nodiscard]] Result_t getResult ()
    {
        return this->m_reduce_data.value(this->m_reduce_op);
    }
};
 
}
 
#endif