amrex/doxygen/AMReX__ParReduce_8H_source.html

#ifndef AMREX_PAR_REDUCE_H_

#define AMREX_PAR_REDUCE_H_

#include <AMReX_Config.H>


#include <AMReX_Concepts.H>

#include <AMReX_Reduce.H>


namespace amrex {


template <typename... Ops, typename... Ts, BaseFabType FAB, typename F>

typename ReduceData<Ts...>::Type


ParReduce (TypeList<Ops...> operation_list, TypeList<Ts...> type_list,

           FabArray<FAB> const& fa, IntVect const& nghost, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    ReduceOps<Ops...> reduce_op;

    ReduceData<Ts...> reduce_data(reduce_op);

    reduce_op.eval(fa, nghost, reduce_data, std::forward<F>(f));

    return reduce_data.value(reduce_op);

}


template <typename Op, typename T, BaseFabType FAB, typename F>

T


ParReduce (TypeList<Op> operation_list, TypeList<T> type_list,

           FabArray<FAB> const& fa, IntVect const& nghost, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    ReduceOps<Op> reduce_op;

    ReduceData<T> reduce_data(reduce_op);

    reduce_op.eval(fa, nghost, reduce_data, std::forward<F>(f));

    auto const& hv = reduce_data.value(reduce_op);

    return amrex::get<0>(hv);

}


template <typename... Ops, typename... Ts, BaseFabType FAB, typename F>

typename ReduceData<Ts...>::Type


ParReduce (TypeList<Ops...> operation_list, TypeList<Ts...> type_list,

           FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    ReduceOps<Ops...> reduce_op;

    ReduceData<Ts...> reduce_data(reduce_op);

    reduce_op.eval(fa, nghost, ncomp, reduce_data, std::forward<F>(f));

    return reduce_data.value(reduce_op);

}


template <typename Op, typename T, BaseFabType FAB, typename F>

T


ParReduce (TypeList<Op> operation_list, TypeList<T> type_list,

           FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    ReduceOps<Op> reduce_op;

    ReduceData<T> reduce_data(reduce_op);

    reduce_op.eval(fa, nghost, ncomp, reduce_data, std::forward<F>(f));

    auto const& hv = reduce_data.value(reduce_op);

    return amrex::get<0>(hv);

}


template <typename... Ops, typename... Ts, BaseFabType FAB, typename F>

typename ReduceData<Ts...>::Type


ParReduce (TypeList<Ops...> operation_list, TypeList<Ts...> type_list,

           FabArray<FAB> const& fa, F&& f)

{

    return ParReduce(operation_list, type_list, fa, IntVect(0), std::forward<F>(f));

}


template <typename Op, typename T, BaseFabType FAB, typename F>

T


ParReduce (TypeList<Op> operation_list, TypeList<T> type_list,

           FabArray<FAB> const& fa, F&& f)

{

    return ParReduce(operation_list, type_list, fa, IntVect(0), std::forward<F>(f));

}


namespace Reduce::detail {


// Reduce f(box_no,i,j,k,n) separately for each component n (PerBox false) or

// for each local box and component (PerBox true, result index box_no*ncomp+n).

// A class template, so that nvcc can name the function of the device lambdas.

template <bool PerBox, typename... Ops>

struct ParReduceComp

{

    template <typename RT, BaseFabType FAB, typename F>

    static Vector<RT>

    eval (FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F const& f)

    {

        const int nlocal = fa.local_size();

        const int nkeys = PerBox ? nlocal*ncomp : ncomp;

        Vector<RT> r(std::max(nkeys,0));

        for (auto& x : r) { for_each_init<0, RT, Ops...>(x); }

        if (nlocal == 0 || ncomp <= 0) { return r; }


#ifdef AMREX_USE_GPU

        // Block b works on component n = b/nb1 of local box L and stores its

        // partial result in slot b, so each key's slots are contiguous.  Box L

        // owns blocks [hoff[L],hoff[L+1]) of every component.

        constexpr int MT = AMREX_GPU_MAX_THREADS;

        const int maxblocks = int(Gpu::Device::maxBlocksPerLaunch());

        Long ntot = 0;

        for (int L = 0; L < nlocal; ++L) {

            Box b = fa.box(fa.IndexArray()[L]);

            ntot += b.ok() ? amrex::grow(b,nghost).numPts() : Long(0);

        }

        // About one launch worth of blocks, with at least 4 cells per thread.

        const Long cells_per_block = std::max(Long(MT)*4, (ntot*ncomp + maxblocks - 1) / maxblocks);

        Gpu::PinnedVector<int> hoff(nlocal+1, 0);

        for (int L = 0; L < nlocal; ++L) {

            Box b = fa.box(fa.IndexArray()[L]);

            Long const npts = b.ok() ? amrex::grow(b,nghost).numPts() : Long(0);

            Long const nb = (npts + cells_per_block - 1) / cells_per_block;

            AMREX_ALWAYS_ASSERT(Long(hoff[L]) + nb < Long(std::numeric_limits<int>::max()));

            hoff[L+1] = hoff[L] + int(nb);

        }

        const int nb1 = hoff[nlocal];

        AMREX_ALWAYS_ASSERT(Long(nb1)*ncomp < Long(std::numeric_limits<int>::max()));

        const int nblocks = nb1*ncomp;

        if (nblocks == 0) { return r; }


        auto const& stream = Gpu::gpuStream();

        Gpu::DeviceVector<int> doff(nlocal+1);

        Gpu::copyAsync(Gpu::hostToDevice, hoff.begin(), hoff.end(), doff.begin());

        Gpu::DeviceVector<RT> dpart(nblocks);

        Gpu::DeviceVector<RT> dr(nkeys);

        int const* poff = doff.data();

        RT* ppart = dpart.data();

        RT* pr = dr.data();

        BoxIndexer const* pbx = fa.getParForInfo(nghost).getBoxes();


#ifdef AMREX_USE_SYCL

        constexpr std::size_t shared_mem_bytes = sizeof(RT)

            * std::max(Gpu::Device::warp_size, MT/Gpu::Device::warp_size);

        amrex::launch<MT>(std::min(nblocks,maxblocks), 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<MT>(std::min(nblocks,maxblocks), 0, stream,

        [=] AMREX_GPU_DEVICE () noexcept

        {

#endif

            for (int b = int(blockIdx.x); b < nblocks; b += int(gridDim.x)) {

                const int n = b / nb1;

                const int bb = b - n*nb1;

                const int L = amrex::bisect(poff, 0, nlocal, bb);

                const int nb = poff[L+1] - poff[L];

                BoxIndexer const& indexer = pbx[L];

                RT t;

                for_each_init<0, RT, Ops...>(t);

                for (auto icell = std::uint64_t(bb-poff[L])*MT + threadIdx.x;

                     icell < indexer.numPts(); icell += std::uint64_t(nb)*MT)

                {

                    auto [i, j, k] = indexer(icell);

                    for_each_local<0, RT, Ops...>(t, f(L,i,j,k,n));

                }

                RT& dst = ppart[b];

                if (threadIdx.x == 0) { for_each_init<0, RT, Ops...>(dst); }

#ifdef AMREX_USE_SYCL

                for_each_parallel<0, RT, Ops...>(dst, t, gh);

#else

                for_each_parallel<0, RT, Ops...>(dst, t);

#endif

            }

        });


        // Combine each key's partial results in a fixed order.

#ifdef AMREX_USE_SYCL

        amrex::launch<MT>(std::min(nkeys,maxblocks), 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<MT>(std::min(nkeys,maxblocks), 0, stream,

        [=] AMREX_GPU_DEVICE () noexcept

        {

#endif

            for (int key = int(blockIdx.x); key < nkeys; key += int(gridDim.x)) {

                int begin = 0;

                int end = 0;

                if (PerBox) { // not constexpr: nvcc cannot first-capture ncomp there

                    const int L = key / ncomp;

                    const int n = key - L*ncomp;

                    begin = n*nb1 + poff[L];

                    end   = n*nb1 + poff[L+1];

                } else {

                    begin = key*nb1;

                    end   = begin + nb1;

                }

                RT t;

                for_each_init<0, RT, Ops...>(t);

                for (int ib = begin + int(threadIdx.x); ib < end; ib += MT) {

                    for_each_local<0, RT, Ops...>(t, ppart[ib]);

                }

                RT& dst = pr[key];

                if (threadIdx.x == 0) { for_each_init<0, RT, Ops...>(dst); }

#ifdef AMREX_USE_SYCL

                for_each_parallel<0, RT, Ops...>(dst, t, gh);

#else

                for_each_parallel<0, RT, Ops...>(dst, t);

#endif

            }

        });


        Gpu::dtoh_memcpy_async(r.data(), pr, sizeof(RT)*nkeys);

        Gpu::streamSynchronize();

#else

        if constexpr (PerBox) {

            // With too few (box, component) pairs for the threads, split each

            // box along its last dimension into contiguous chunks of about 16K

            // cells. Items are ordered by box, component and chunk, which is

            // the memory order, and their partial results are combined in order.

            constexpr int d = AMREX_SPACEDIM-1;

            const int nthreads = OpenMP::in_parallel() ? 1 : OpenMP::get_max_threads();

            const bool split = nthreads > 1 && nkeys < 4*nthreads;

            Vector<Box> bx(nlocal);

            Vector<int> thick(nlocal), nchunk(nlocal), ioff(nlocal+1, 0);

            for (int L = 0; L < nlocal; ++L) {

                bx[L] = fa.box(fa.IndexArray()[L]);

                if (bx[L].ok()) { bx[L].grow(nghost); }

                // An empty box is one empty chunk.

                const int len = bx[L].ok() ? bx[L].length(d) : 1;

                thick[L] = (split && bx[L].ok()) ? int(std::max(Long(1), Long(16384)/(bx[L].numPts()/len))) : len;

                nchunk[L] = (len + thick[L] - 1) / thick[L];

                ioff[L+1] = ioff[L] + nchunk[L]*ncomp;

            }

            // Without splitting, item is the result index.

            const int nitems = ioff[nlocal];

            Vector<RT> part(split ? nitems : 0);

            RT* pr = split ? part.data() : r.data();

#ifdef AMREX_USE_OMP

#pragma omp parallel for schedule(dynamic) if (!system::regtest_reduction)

#endif

            for (int item = 0; item < nitems; ++item) {

                const int L = split ? amrex::bisect(ioff.data(), 0, nlocal, item) : item/ncomp;

                const int m = item - ioff[L];

                const int n = m / nchunk[L];

                const int c = m - n*nchunk[L];

                Box b = bx[L];

                b.setSmall(d, bx[L].smallEnd(d) + c*thick[L]);

                b.setBig(d, std::min(b.smallEnd(d)+thick[L]-1, bx[L].bigEnd(d)));

                const auto lo = amrex::lbound(b);

                const auto hi = amrex::ubound(b);

                RT rem;

                for_each_init<0, RT, Ops...>(rem);

                Lanes<RT> acc;

                lanes_set<0>(acc, rem);

                for (int k = lo.z; k <= hi.z; ++k) {

                for (int j = lo.y; j <= hi.y; ++j) {

                    reduce_row<Ops...>(lo.x, hi.x, acc, rem,

                        [&] (int i) { return f(L,i,j,k,n); });

                }}

                RT& x = pr[item];

                for_each_init<0, RT, Ops...>(x);

                lanes_fold<0, RT, Ops...>(acc, x);

                for_each_local<0, RT, Ops...>(x, rem);

            }

            if (split) {

                for (int L = 0; L < nlocal; ++L) {

                    for (int m = 0; m < nchunk[L]*ncomp; ++m) {

                        for_each_local<0, RT, Ops...>(r[L*ncomp+m/nchunk[L]], part[ioff[L]+m]);

                    }

                }

            }

        } else {

            const int nthreads = OpenMP::in_parallel() ? 1 : OpenMP::get_max_threads();

            Vector<RT> rt(std::size_t(nthreads)*ncomp);

            for (auto& x : rt) { for_each_init<0, RT, Ops...>(x); }

#ifdef AMREX_USE_OMP

#pragma omp parallel if (nthreads > 1 && !system::regtest_reduction)

#endif

            {

                RT* rr = rt.data() + std::size_t(nthreads > 1 ? OpenMP::get_thread_num() : 0)*ncomp;

                // Component loop outermost, so that each pass streams one component.

                for (int n = 0; n < ncomp; ++n) {

                    RT rem;

                    for_each_init<0, RT, Ops...>(rem);

                    Lanes<RT> acc;

                    lanes_set<0>(acc, rem);

                    for (MFIter mfi(fa,true); mfi.isValid(); ++mfi) {

                        Box const& b = mfi.growntilebox(nghost);

                        const int L = 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_row<Ops...>(lo.x, hi.x, acc, rem,

                                [&] (int i) { return f(L,i,j,k,n); });

                        }}

                    }

                    lanes_fold<0, RT, Ops...>(acc, rr[n]);

                    for_each_local<0, RT, Ops...>(rr[n], rem);

                }

            }

            for (int it = 0; it < nthreads; ++it) {

                for (int n = 0; n < ncomp; ++n) {

                    for_each_local<0, RT, Ops...>(r[n], rt[std::size_t(it)*ncomp+n]);

                }

            }

        }

#endif

        return r;

    }

};


}


template <typename... Ops, typename... Ts, BaseFabType FAB, typename F>

Vector<typename ReduceData<Ts...>::Type>


ParReducePerComp (TypeList<Ops...> operation_list, TypeList<Ts...> type_list,

                  FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    return Reduce::detail::ParReduceComp<false,Ops...>::template

        eval<typename ReduceData<Ts...>::Type>(fa, nghost, ncomp, std::forward<F>(f));

}


template <typename Op, typename T, BaseFabType FAB, typename F>

Vector<T>


ParReducePerComp (TypeList<Op> operation_list, TypeList<T> type_list,

                  FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    auto const& rt = Reduce::detail::ParReduceComp<false,Op>::template

        eval<typename ReduceData<T>::Type>(fa, nghost, ncomp, std::forward<F>(f));

    Vector<T> r(rt.size());

    for (int i = 0, N = int(rt.size()); i < N; ++i) { r[i] = amrex::get<0>(rt[i]); }

    return r;

}


template <typename... Ops, typename... Ts, BaseFabType FAB, typename F>

Vector<typename ReduceData<Ts...>::Type>


ParReducePerBoxComp (TypeList<Ops...> operation_list, TypeList<Ts...> type_list,

                     FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    return Reduce::detail::ParReduceComp<true,Ops...>::template

        eval<typename ReduceData<Ts...>::Type>(fa, nghost, ncomp, std::forward<F>(f));

}


template <typename Op, typename T, BaseFabType FAB, typename F>

Vector<T>


ParReducePerBoxComp (TypeList<Op> operation_list, TypeList<T> type_list,

                     FabArray<FAB> const& fa, IntVect const& nghost, int ncomp, F&& f)

{

    amrex::ignore_unused(operation_list,type_list);

    auto const& rt = Reduce::detail::ParReduceComp<true,Op>::template

        eval<typename ReduceData<T>::Type>(fa, nghost, ncomp, std::forward<F>(f));

    Vector<T> r(rt.size());

    for (int i = 0, N = int(rt.size()); i < N; ++i) { r[i] = amrex::get<0>(rt[i]); }

    return r;

}


}

#endif

AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50

AMReX_Concepts.H

AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18

AMReX_Reduce.H

amrex::BoxND::setBig
__host__ __device__ BoxND & setBig(const IntVectND< dim > &bg) noexcept
Redefine the big end of the BoxND.
Definition AMReX_Box.H:516

amrex::BoxND::ok
__host__ __device__ bool ok() const noexcept
Return true if high bounds are >= low bounds and the index type is valid.
Definition AMReX_Box.H:229

amrex::BoxND::setSmall
__host__ __device__ BoxND & setSmall(const IntVectND< dim > &sm) noexcept
Definition AMReX_Box.H:508

amrex::FabArray
An Array of FortranArrayBox(FAB)-like Objects.
Definition AMReX_FabArray.H:356

amrex::ReduceData
Definition AMReX_Reduce.H:556

amrex::ReduceData::value
Type value()
Definition AMReX_Reduce.H:591

amrex::ReduceOps
Definition AMReX_Reduce.H:715

amrex::ReduceOps::eval
void eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition AMReX_Reduce.H:852

amrex::Vector
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29

amrex::BaseFabType
Checks if a type is derived from amrex::BaseFab.
Definition AMReX_Concepts.H:13

amrex::Long
amrex_long Long
Definition AMReX_INT.H:30

amrex::ubound
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
Return the inclusive upper bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1365

amrex::lbound
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Return the inclusive lower bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1351

amrex::grow
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Return a copy of b grown uniformly by i cells in every direction.
Definition AMReX_Box.H:1328

amrex::bisect
__host__ __device__ T bisect(T lo, T hi, F f, T tol=1e-12, int max_iter=100)
Find a root of a scalar function on a bracketing interval using bisection.
Definition AMReX_Algorithm.H:170

amrex
Definition AMReX_Amr.cpp:50

amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
No-op helper that marks variables as intentionally unused.
Definition AMReX.H:273

amrex::Order::F
@ F
Definition AMReX_Order.H:7

amrex::ParReduce
ReduceData< Ts... >::Type ParReduce(TypeList< Ops... > operation_list, TypeList< Ts... > type_list, FabArray< FAB > const &fa, IntVect const &nghost, F &&f)
Parallel reduce for MultiFab/FabArray. The reduce result is local and it's the user's responsibility ...
Definition AMReX_ParReduce.H:48

amrex::ParReducePerComp
Vector< typename ReduceData< Ts... >::Type > ParReducePerComp(TypeList< Ops... > operation_list, TypeList< Ts... > type_list, FabArray< FAB > const &fa, IntVect const &nghost, int ncomp, F &&f)
Parallel reduce for MultiFab/FabArray with a separate result for each component. The reduce result is...
Definition AMReX_ParReduce.H:563

amrex::begin
__host__ __device__ Dim3 begin(BoxND< dim > const &box) noexcept
Return the iterator begin coordinate of box as Dim3.
Definition AMReX_Box.H:2241

amrex::launch
void launch(T const &n, L const &f) noexcept
Definition AMReX_GpuLaunchFunctsC.H:122

amrex::Box
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:35

amrex::Direction::x
@ x
Definition AMReX_Orientation.H:14

amrex::IntVect
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:38

amrex::BoxIndexer
BoxIndexerND< 3 > BoxIndexer
Definition AMReX_Box.H:2610

amrex::int
const int[]
Definition AMReX_BLProfiler.cpp:1665

amrex::ParReducePerBoxComp
Vector< typename ReduceData< Ts... >::Type > ParReducePerBoxComp(TypeList< Ops... > operation_list, TypeList< Ts... > type_list, FabArray< FAB > const &fa, IntVect const &nghost, int ncomp, F &&f)
Parallel reduce for MultiFab/FabArray with a separate result for each local box and component....
Definition AMReX_ParReduce.H:600

amrex::end
__host__ __device__ Dim3 end(BoxND< dim > const &box) noexcept
Return the iterator end coordinate of box as Dim3.
Definition AMReX_Box.H:2259

amrex::split
std::vector< std::string > split(std::string const &s, std::string const &sep)
Split a string using given tokens in sep.
Definition AMReX_String.cpp:44

amrex::get
__host__ __device__ constexpr int get(IntVectND< dim > const &iv) noexcept
Get the Idx-th element of IntVectND<dim>.
Definition AMReX_IntVect.H:1360

amrex::TypeList
Struct for holding types.
Definition AMReX_TypeList.H:13