AMReX_BaseFab.H

Go to the documentation of this file.

#ifndef AMREX_BASEFAB_H_
#define AMREX_BASEFAB_H_
#include <AMReX_Config.H>
 
#include <AMReX_Algorithm.H>
#include <AMReX_Extension.H>
#include <AMReX_BLassert.H>
#include <AMReX_Array.H>
#include <AMReX_Box.H>
#include <AMReX_Loop.H>
#include <AMReX_BoxList.H>
#include <AMReX_BArena.H>
#include <AMReX_CArena.H>
#include <AMReX_DataAllocator.H>
#include <AMReX_REAL.H>
#include <AMReX_BLProfiler.H>
#include <AMReX_BoxIterator.H>
#include <AMReX_MakeType.H>
#include <AMReX_Utility.H>
#include <AMReX_Reduce.H>
#include <AMReX_Gpu.H>
#include <AMReX_Scan.H>
#include <AMReX_Math.H>
#include <AMReX_OpenMP.H>
#include <AMReX_MemPool.H>
 
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <limits>
#include <climits>
#include <array>
#include <type_traits>
#include <memory>
#include <atomic>
#include <utility>
 
namespace amrex
{
 
extern std::atomic<Long> atomic_total_bytes_allocated_in_fabs;
extern std::atomic<Long> atomic_total_bytes_allocated_in_fabs_hwm;
extern std::atomic<Long> atomic_total_cells_allocated_in_fabs;
extern std::atomic<Long> atomic_total_cells_allocated_in_fabs_hwm;
extern Long private_total_bytes_allocated_in_fabs;     
extern Long private_total_bytes_allocated_in_fabs_hwm; 
extern Long private_total_cells_allocated_in_fabs;     
extern Long private_total_cells_allocated_in_fabs_hwm; 
#ifdef AMREX_USE_OMP
#pragma omp threadprivate(private_total_bytes_allocated_in_fabs)
#pragma omp threadprivate(private_total_bytes_allocated_in_fabs_hwm)
#pragma omp threadprivate(private_total_cells_allocated_in_fabs)
#pragma omp threadprivate(private_total_cells_allocated_in_fabs_hwm)
#endif
 
Long TotalBytesAllocatedInFabs () noexcept;
Long TotalBytesAllocatedInFabsHWM () noexcept;
Long TotalCellsAllocatedInFabs () noexcept;
Long TotalCellsAllocatedInFabsHWM () noexcept;
void ResetTotalBytesAllocatedInFabsHWM () noexcept;
void update_fab_stats (Long n, Long s, std::size_t szt) noexcept;
 
void BaseFab_Initialize ();
void BaseFab_Finalize ();
 
struct SrcComp {
    AMREX_GPU_HOST_DEVICE
    explicit SrcComp (int ai) noexcept : i(ai) {}
    int i;
};
 
struct DestComp {
    AMREX_GPU_HOST_DEVICE
    explicit DestComp (int ai) noexcept : i(ai) {}
    int i;
};
 
struct NumComps {
    AMREX_GPU_HOST_DEVICE
    explicit NumComps (int an) noexcept : n(an) {}
    int n;
};
 
template <typename T>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Array4<T>
makeArray4 (T* p, Box const& bx, int ncomp) noexcept
{
    return Array4<T>{p, amrex::begin(bx), amrex::end(bx), ncomp};
}
 
template <typename T>
std::enable_if_t<std::is_arithmetic_v<T>>
placementNew (T* const /*ptr*/, Long /*n*/)
{}
 
template <typename T>
std::enable_if_t<std::is_trivially_default_constructible_v<T>
                 && !std::is_arithmetic_v<T>>
placementNew (T* const ptr, Long n)
{
    for (Long i = 0; i < n; ++i) {
        new (ptr+i) T;
    }
}
 
template <typename T>
std::enable_if_t<!std::is_trivially_default_constructible_v<T>>
placementNew (T* const ptr, Long n)
{
    AMREX_HOST_DEVICE_FOR_1D ( n, i,
    {
        new (ptr+i) T;
    });
}
 
template <typename T>
std::enable_if_t<std::is_trivially_destructible_v<T>>
placementDelete (T* const /*ptr*/, Long /*n*/)
{}
 
template <typename T>
std::enable_if_t<!std::is_trivially_destructible_v<T>>
placementDelete (T* const ptr, Long n)
{
    AMREX_HOST_DEVICE_FOR_1D (n, i,
    {
        (ptr+i)->~T();
    });
}
 
template <class T>
class BaseFab
    : public DataAllocator
{
public:
 
    template <class U> friend class BaseFab;
 
    using value_type = T;
 
    BaseFab () noexcept = default;
 
    explicit BaseFab (Arena* ar) noexcept;
 
    BaseFab (const Box& bx, int n, Arena* ar);
 
    explicit BaseFab (const Box& bx, int n = 1, bool alloc = true,
                      bool shared = false, Arena* ar = nullptr);
 
    BaseFab (const BaseFab<T>& rhs, MakeType make_type, int scomp, int ncomp);
 
    BaseFab (const Box& bx, int ncomp, T* p);
    BaseFab (const Box& bx, int ncomp, T const* p);
 
    explicit BaseFab (Array4<T> const& a) noexcept;
 
    explicit BaseFab (Array4<T> const& a, IndexType t) noexcept;
 
    explicit BaseFab (Array4<T const> const& a) noexcept;
 
    explicit BaseFab (Array4<T const> const& a, IndexType t) noexcept;
 
    virtual ~BaseFab () noexcept;
 
    BaseFab (const BaseFab<T>& rhs) = delete;
    BaseFab<T>& operator= (const BaseFab<T>& rhs) = delete;
 
    BaseFab (BaseFab<T>&& rhs) noexcept;
    BaseFab<T>& operator= (BaseFab<T>&& rhs) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab& operator= (T const&) noexcept;
 
    static void Initialize();
    static void Finalize();
 
    void resize (const Box& b, int N = 1, Arena* ar = nullptr);
 
    template <class U=T, std::enable_if_t<std::is_trivially_destructible_v<U>,int> = 0>
    [[nodiscard]] Elixir elixir () noexcept;
 
    void clear () noexcept;
 
    [[nodiscard]] std::unique_ptr<T,DataDeleter> release () noexcept;
 
    [[nodiscard]] std::size_t nBytes () const noexcept { return this->truesize*sizeof(T); }
 
    [[nodiscard]] std::size_t nBytesOwned () const noexcept {
        return (this->ptr_owner) ? nBytes() : 0;
    }
 
    [[nodiscard]] std::size_t nBytes (const Box& bx, int ncomps) const noexcept
        { return bx.numPts() * sizeof(T) * ncomps; }
 
    [[nodiscard]] int nComp () const noexcept { return this->nvar; }
 
    [[nodiscard]] const int* nCompPtr() const noexcept {
        return &(this->nvar);
    }
 
    [[nodiscard]] Long numPts () const noexcept { return this->domain.numPts(); }
 
    [[nodiscard]] Long size () const noexcept { return this->nvar*this->domain.numPts(); }
 
    [[nodiscard]] const Box& box () const noexcept { return this->domain; }
 
    [[nodiscard]] IntVect length () const noexcept { return this->domain.length(); }
 
    [[nodiscard]] const IntVect& smallEnd () const noexcept { return this->domain.smallEnd(); }
 
    [[nodiscard]] const IntVect& bigEnd () const noexcept { return this->domain.bigEnd(); }
 
    [[nodiscard]] const int* loVect () const noexcept { return this->domain.loVect(); }
 
    [[nodiscard]] const int* hiVect () const noexcept { return this->domain.hiVect(); }
 
    [[nodiscard]] bool contains (const BaseFab<T>& fab) const noexcept
    {
        return box().contains(fab.box()) && this->nvar <= fab.nvar;
    }
 
    [[nodiscard]] bool contains (const Box& bx) const noexcept { return box().contains(bx); }
 
    [[nodiscard]] T* dataPtr (int n = 0) noexcept {
        if (this->dptr) {
            return &(this->dptr[n*this->domain.numPts()]);
        } else {
            return nullptr;
        }
    }
 
    [[nodiscard]] const T* dataPtr (int n = 0) const noexcept {
        if (this->dptr) {
            return &(this->dptr[n*this->domain.numPts()]);
        } else {
            return nullptr;
        }
    }
 
    [[nodiscard]] T* dataPtr (const IntVect& p, int n = 0) noexcept;
 
    [[nodiscard]] const T* dataPtr (const IntVect& p, int n = 0) const noexcept;
 
    void setPtr (T* p, Long sz) noexcept { AMREX_ASSERT(this->dptr == nullptr && this->truesize == 0); this->dptr = p; this->truesize = sz; }
 
    void prefetchToHost () const noexcept;
    void prefetchToDevice () const noexcept;
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> array () const noexcept
    {
        return makeArray4<T const>(this->dptr, this->domain, this->nvar);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> array (int start_comp) const noexcept
    {
        return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar),start_comp);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> array (int start_comp, int num_comps) const noexcept
    {
        return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T> array () noexcept
    {
        return makeArray4<T>(this->dptr, this->domain, this->nvar);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T> array (int start_comp) noexcept
    {
        return Array4<T>(makeArray4<T>(this->dptr, this->domain, this->nvar),start_comp);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T> array (int start_comp, int num_comps) noexcept
    {
        return Array4<T>(makeArray4<T>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> const_array () const noexcept
    {
        return makeArray4<T const>(this->dptr, this->domain, this->nvar);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> const_array (int start_comp) const noexcept
    {
        return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar),start_comp);
    }
 
    [[nodiscard]] AMREX_FORCE_INLINE
    Array4<T const> const_array (int start_comp, int num_comps) const noexcept
    {
        return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
    }
 
    [[nodiscard]] bool isAllocated () const noexcept { return this->dptr != nullptr; }
 
    [[nodiscard]] T& operator() (const IntVect& p, int N) noexcept;
 
    [[nodiscard]] T& operator() (const IntVect& p) noexcept;
 
    [[nodiscard]] const T& operator() (const IntVect& p, int N) const noexcept;
 
    [[nodiscard]] const T& operator() (const IntVect& p) const noexcept;
 
    void getVal (T* data, const IntVect& pos, int N, int numcomp) const noexcept;
    void getVal (T* data, const IntVect& pos) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on,
#else
    template <RunOn run_on=RunOn::Host,
#endif
              class U=T, std::enable_if_t<std::is_same_v<U,float> || std::is_same_v<U,double>,int> FOO = 0>
    void fill_snan () noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setVal (T const& x, const Box& bx, int dcomp, int ncomp) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setVal (T const& x, const Box& bx, int N = 0) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setVal (T const& x, int N) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setValIfNot (T const& val, const Box& bx, const BaseFab<int>& mask, int nstart, int num) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setComplement (T const& x, const Box& b, int ns, int num) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& copy (const BaseFab<T>& src, const Box& srcbox, int srccomp,
                      const Box& destbox, int destcomp, int numcomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& copy (const BaseFab<T>& src, int srccomp, int destcomp,
                      int numcomp = 1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& copy (const BaseFab<T>& src, const Box& destbox) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    std::size_t copyToMem (const Box& srcbox, int srccomp,
                           int numcomp, void* dst) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on, typename BUF = T>
#else
    template <RunOn run_on=RunOn::Host, typename BUF = T>
#endif
    std::size_t copyFromMem (const Box& dstbox, int dstcomp,
                             int numcomp, const void* src) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on, typename BUF = T>
#else
    template <RunOn run_on=RunOn::Host, typename BUF = T>
#endif
    std::size_t addFromMem (const Box& dstbox, int dstcomp,
                            int numcomp, const void* src) noexcept;
 
    BaseFab<T>& shift (const IntVect& v) noexcept;
    BaseFab<T>& shift (int idir, int n_cell) noexcept;
    BaseFab<T>& shiftHalf (int dir, int n_cell) noexcept;
    BaseFab<T>& shiftHalf (const IntVect& v) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] Real norminfmask (const Box& subbox, const BaseFab<int>& mask, int scomp=0, int ncomp=1) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] Real norm (int p, int scomp = 0, int numcomp = 1) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] Real norm (const Box& subbox, int p, int scomp = 0, int numcomp = 1) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void abs () noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void abs (int comp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void abs (const Box& subbox, int comp = 0, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T min (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T min (const Box& subbox, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T max (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T max (const Box& subbox, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] std::pair<T,T> minmax (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] std::pair<T,T> minmax (const Box& subbox, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T maxabs (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T maxabs (const Box& subbox, int comp = 0) const noexcept;
 
    /*(
    * \return location of given value
    */
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] IntVect indexFromValue (const Box& subbox, int comp, T const& value) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] IntVect minIndex (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] IntVect minIndex (const Box& subbox, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void  minIndex (const Box& subbox, Real& min_val, IntVect& min_idx, int comp = 0) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] IntVect maxIndex (int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] IntVect maxIndex (const Box& subbox, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void  maxIndex (const Box& subbox, Real& max_value, IntVect& max_idx, int comp = 0) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    int maskLT (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    int maskLE (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    int maskEQ (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    int maskGT (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    int maskGE (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T sum (int comp, int numcomp = 1) const noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T sum (const Box& subbox, int comp, int numcomp = 1) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& invert (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& invert (T const& r, int comp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& negate (const Box& b, int comp=0, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& negate (int comp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (T const& r, int comp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                      int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& atomicAdd (const BaseFab<T>& x) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& atomicAdd (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& atomicAdd (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                           int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& atomicAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                           int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& lockAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                         int srccomp, int destcomp, int numcomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& saxpy (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
                       int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& saxpy (T a, const BaseFab<T>& x) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& xpay (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
                      int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& addproduct (const Box& destbox, int destcomp, int numcomp,
                            const BaseFab<T>& src1, int comp1,
                            const BaseFab<T>& src2, int comp2) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                       int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                       int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (T const& r, int comp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                      int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                      int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (T const& r, int comp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                        int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                        int srccomp, int destcomp, int numcomp=1) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& protected_divide (const BaseFab<T>& src) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& protected_divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& protected_divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                                  int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& protected_divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                                  int srccomp, int destcomp, int numcomp=1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& linInterp (const BaseFab<T>& f1, const Box& b1, int comp1,
                           const BaseFab<T>& f2, const Box& b2, int comp2,
                           Real t1, Real t2, Real t,
                           const Box& b, int comp, int numcomp = 1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& linInterp (const BaseFab<T>& f1, int comp1,
                           const BaseFab<T>& f2, int comp2,
                           Real t1, Real t2, Real t,
                           const Box& b, int comp, int numcomp = 1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& linComb (const BaseFab<T>& f1, const Box& b1, int comp1,
                         const BaseFab<T>& f2, const Box& b2, int comp2,
                         Real alpha, Real beta, const Box& b,
                         int comp, int numcomp = 1) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dot (const Box& xbx, int xcomp, const BaseFab<T>& y, const Box& ybx, int ycomp,
           int numcomp = 1) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dotmask (const BaseFab<int>& mask, const Box& xbx, int xcomp,
               const BaseFab<T>& y, const Box& ybx, int ycomp,
               int numcomp) const noexcept;
 
    void SetBoxType (const IndexType& typ) noexcept { this->domain.setType(typ); }
 
    //
    // New interfaces
    //
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setVal (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setVal (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setValIf (T const& val, const BaseFab<int>& mask) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setValIf (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setValIfNot (T const& val, const BaseFab<int>& mask) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setValIfNot (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    void setComplement (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& copy (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& copy (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator+= (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator+= (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& plus (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator-= (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator-= (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& minus (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator*= (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator*= (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& mult (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator/= (T const& val) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& operator/= (const BaseFab<T>& src) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& divide (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& negate () noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& negate (const Box& bx, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& invert (T const& r) noexcept;
    //
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    BaseFab<T>& invert (T const& r, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T sum (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dot (const BaseFab<T>& src, const Box& bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dot (const Box& bx, int destcomp, int numcomp) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dot (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept;
 
#if defined(AMREX_USE_GPU)
    template <RunOn run_on>
#else
    template <RunOn run_on=RunOn::Host>
#endif
    [[nodiscard]] T dotmask (const BaseFab<T>& src, const Box& bx, const BaseFab<int>& mask,
               SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept;
 
protected:
    void define ();
 
    T*   dptr = nullptr;        
    Box  domain;                
    int  nvar = 0;              
    Long truesize = 0L;         
    bool ptr_owner = false;     
    bool shared_memory = false; 
#ifdef AMREX_USE_GPU
    gpuStream_t alloc_stream{};
#endif
};
 
template <class T>
AMREX_FORCE_INLINE
T*
BaseFab<T>::dataPtr (const IntVect& p, int n) noexcept
{
    AMREX_ASSERT(n >= 0);
    AMREX_ASSERT(n < this->nvar);
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr + (this->domain.index(p)+n*this->domain.numPts());
}
 
template <class T>
AMREX_FORCE_INLINE
const T*
BaseFab<T>::dataPtr (const IntVect& p, int n) const noexcept
{
    AMREX_ASSERT(n >= 0);
    AMREX_ASSERT(n < this->nvar);
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr + (this->domain.index(p)+n*this->domain.numPts());
}
 
template <class T>
void
BaseFab<T>::prefetchToHost () const noexcept
{
#ifdef AMREX_USE_GPU
    if (this->arena()->isManaged()) {
#if defined(AMREX_USE_SYCL)
        // xxxxx SYCL todo: prefetchToHost
        // std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
        // auto& q = Gpu::Device::streamQueue();
        // q.submit([&] (sycl::handler& h) { h.prefetch(this->dptr, s); });
#elif defined(AMREX_USE_CUDA) && !defined(_WIN32)
        if (Gpu::Device::devicePropMajor() >= 6) {
            std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
            AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s,
                                                      cudaCpuDeviceId,
                                                      Gpu::gpuStream()));
        }
#elif defined(AMREX_USE_HIP)
        // xxxxx HIP FIX HERE after managed memory is supported
#endif
    }
#endif
}
 
template <class T>
void
BaseFab<T>::prefetchToDevice () const noexcept
{
#ifdef AMREX_USE_GPU
    if (this->arena()->isManaged()) {
#if defined(AMREX_USE_SYCL)
        std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
        auto& q = Gpu::Device::streamQueue();
        q.submit([&] (sycl::handler& h) { h.prefetch(this->dptr, s); });
#elif defined(AMREX_USE_CUDA) && !defined(_WIN32)
        if (Gpu::Device::devicePropMajor() >= 6) {
            std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
            AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s,
                                                      Gpu::Device::deviceId(),
                                                      Gpu::gpuStream()));
        }
#elif defined(AMREX_USE_HIP)
        // xxxxx HIP FIX HERE after managed memory is supported
#endif
    }
#endif
}
 
template <class T>
AMREX_FORCE_INLINE
T&
BaseFab<T>::operator() (const IntVect& p, int n) noexcept
{
    AMREX_ASSERT(n >= 0);
    AMREX_ASSERT(n < this->nvar);
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr[this->domain.index(p)+n*this->domain.numPts()];
}
 
template <class T>
AMREX_FORCE_INLINE
T&
BaseFab<T>::operator() (const IntVect& p) noexcept
{
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr[this->domain.index(p)];
}
 
template <class T>
AMREX_FORCE_INLINE
const T&
BaseFab<T>::operator() (const IntVect& p, int n) const noexcept
{
    AMREX_ASSERT(n >= 0);
    AMREX_ASSERT(n < this->nvar);
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr[this->domain.index(p)+n*this->domain.numPts()];
}
 
template <class T>
AMREX_FORCE_INLINE
const T&
BaseFab<T>::operator() (const IntVect& p) const noexcept
{
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(this->domain.contains(p));
 
    return this->dptr[this->domain.index(p)];
}
 
template <class T>
void
BaseFab<T>::getVal  (T*             data,
                     const IntVect& pos,
                     int            n,
                     int            numcomp) const noexcept
{
    const int loc      = this->domain.index(pos);
    const Long sz    = this->domain.numPts();
 
    AMREX_ASSERT(!(this->dptr == nullptr));
    AMREX_ASSERT(n >= 0 && n + numcomp <= this->nvar);
 
    for (int k = 0; k < numcomp; k++) {
        data[k] = this->dptr[loc+(n+k)*sz];
    }
}
 
template <class T>
void
BaseFab<T>::getVal (T*             data,
                    const IntVect& pos) const noexcept
{
    getVal(data,pos,0,this->nvar);
}
 
template <class T>
BaseFab<T>&
BaseFab<T>::shift (const IntVect& v) noexcept
{
    this->domain += v;
    return *this;
}
 
template <class T>
BaseFab<T>&
BaseFab<T>::shift (int idir, int n_cell) noexcept
{
    this->domain.shift(idir,n_cell);
    return *this;
}
 
template <class T>
BaseFab<T> &
BaseFab<T>::shiftHalf (const IntVect& v) noexcept
{
    this->domain.shiftHalf(v);
    return *this;
}
 
template <class T>
BaseFab<T> &
BaseFab<T>::shiftHalf (int idir, int n_cell) noexcept
{
    this->domain.shiftHalf(idir,n_cell);
    return *this;
}
 
template <class T>
template <RunOn run_on, class U,
          std::enable_if_t<std::is_same_v<U,float> || std::is_same_v<U,double>, int> FOO>
void
BaseFab<T>::fill_snan () noexcept
{
    amrex::fill_snan<run_on>(this->dptr, this->truesize);
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setVal (T const& x, const Box& bx, int n) noexcept
{
    this->setVal<run_on>(x, bx, DestComp{n}, NumComps{1});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setVal (T const& x, int n) noexcept
{
    this->setVal<run_on>(x, this->domain, DestComp{n}, NumComps{1});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setVal (T const& x, const Box& bx, int dcomp, int ncomp) noexcept
{
    this->setVal<run_on>(x, bx, DestComp{dcomp}, NumComps{ncomp});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setValIfNot (T const& val, const Box& bx, const BaseFab<int>& mask, int ns, int num) noexcept
{
    this->setValIfNot<run_on>(val, bx, mask, DestComp{ns}, NumComps{num});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src, const Box& srcbox, int srccomp,
                  const Box& destbox, int destcomp, int numcomp) noexcept
{
    AMREX_ASSERT(destbox.ok());
    AMREX_ASSERT(srcbox.sameSize(destbox));
    AMREX_ASSERT(src.box().contains(srcbox));
    AMREX_ASSERT(this->domain.contains(destbox));
    AMREX_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    AMREX_ASSERT(destcomp >= 0 && destcomp+numcomp <= this->nvar);
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
 
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) = s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src, const Box& destbox) noexcept
{
    return this->copy<run_on>(src, destbox, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    return copy<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
void
BaseFab<T>::define ()
{
    AMREX_ASSERT(this->dptr == nullptr);
    AMREX_ASSERT(this->domain.numPts() > 0);
    AMREX_ASSERT(this->nvar >= 0);
    if (this->nvar == 0) { return; }
    AMREX_ASSERT(std::numeric_limits<Long>::max()/this->nvar > this->domain.numPts());
 
    this->truesize  = this->nvar*this->domain.numPts();
    this->ptr_owner = true;
    this->dptr = static_cast<T*>(this->alloc(this->truesize*sizeof(T)));
#ifdef AMREX_USE_GPU
    this->alloc_stream = Gpu::gpuStream();
#endif
 
    placementNew(this->dptr, this->truesize);
 
    amrex::update_fab_stats(this->domain.numPts(), this->truesize, sizeof(T));
 
    if constexpr (std::is_same_v<T,float> || std::is_same_v<T,double>) {
        if (amrex::InitSNaN() && this->truesize > 0) {
#ifdef AMREX_USE_GPU
            if (Gpu::inLaunchRegion() && arena()->isDeviceAccessible()) {
                this->template fill_snan<RunOn::Device>();
                Gpu::streamSynchronize();
            } else
#endif
            {
                this->template fill_snan<RunOn::Host>();
            }
        }
    }
}
 
template <class T>
BaseFab<T>::BaseFab (Arena* ar) noexcept
    : DataAllocator{ar}
{}
 
template <class T>
BaseFab<T>::BaseFab (const Box& bx, int n, Arena* ar)
    : DataAllocator{ar}, domain(bx), nvar(n)
{
    define();
}
 
template <class T>
BaseFab<T>::BaseFab (const Box& bx, int n, bool alloc, bool shared, Arena* ar)
    : DataAllocator{ar}, domain(bx), nvar(n), shared_memory(shared)
{
    if (!this->shared_memory && alloc) { define(); }
}
 
template <class T>
BaseFab<T>::BaseFab (const BaseFab<T>& rhs, MakeType make_type, int scomp, int ncomp)
    : DataAllocator{rhs.arena()},
      dptr(const_cast<T*>(rhs.dataPtr(scomp))),
      domain(rhs.domain), nvar(ncomp),
      truesize(ncomp*rhs.domain.numPts())
{
    AMREX_ASSERT(scomp+ncomp <= rhs.nComp());
    if (make_type == amrex::make_deep_copy)
    {
        this->dptr = nullptr;
        define();
        this->copy<RunOn::Device>(rhs, this->domain, scomp, this->domain, 0, ncomp);
    } else if (make_type == amrex::make_alias) {
        ; // nothing to do
    } else {
        amrex::Abort("BaseFab: unknown MakeType");
    }
}
 
template<class T>
BaseFab<T>::BaseFab (const Box& bx, int ncomp, T* p)
    : dptr(p), domain(bx), nvar(ncomp), truesize(bx.numPts()*ncomp)
{
}
 
template<class T>
BaseFab<T>::BaseFab (const Box& bx, int ncomp, T const* p)
    : dptr(const_cast<T*>(p)), domain(bx), nvar(ncomp), truesize(bx.numPts()*ncomp)
{
}
 
template<class T>
BaseFab<T>::BaseFab (Array4<T> const& a) noexcept
    : dptr(a.p),
      domain(IntVect(AMREX_D_DECL(a.begin.x,a.begin.y,a.begin.z)),
             IntVect(AMREX_D_DECL(a.end.x-1,a.end.y-1,a.end.z-1))),
      nvar(a.ncomp), truesize(a.ncomp*a.nstride)
{}
 
template<class T>
BaseFab<T>::BaseFab (Array4<T> const& a, IndexType t) noexcept
    : dptr(a.p),
      domain(IntVect(AMREX_D_DECL(a.begin.x,a.begin.y,a.begin.z)),
             IntVect(AMREX_D_DECL(a.end.x-1,a.end.y-1,a.end.z-1)), t),
      nvar(a.ncomp), truesize(a.ncomp*a.nstride)
{}
 
template<class T>
BaseFab<T>::BaseFab (Array4<T const> const& a) noexcept
    : dptr(const_cast<T*>(a.p)),
      domain(IntVect(AMREX_D_DECL(a.begin.x,a.begin.y,a.begin.z)),
             IntVect(AMREX_D_DECL(a.end.x-1,a.end.y-1,a.end.z-1))),
      nvar(a.ncomp), truesize(a.ncomp*a.nstride)
{}
 
template<class T>
BaseFab<T>::BaseFab (Array4<T const> const& a, IndexType t) noexcept
    : dptr(const_cast<T*>(a.p)),
      domain(IntVect(AMREX_D_DECL(a.begin.x,a.begin.y,a.begin.z)),
             IntVect(AMREX_D_DECL(a.end.x-1,a.end.y-1,a.end.z-1)), t),
      nvar(a.ncomp), truesize(a.ncomp*a.nstride)
{}
 
template <class T>
BaseFab<T>::~BaseFab () noexcept
{
    clear();
}
 
template <class T>
BaseFab<T>::BaseFab (BaseFab<T>&& rhs) noexcept
    : DataAllocator{rhs.arena()},
      dptr(rhs.dptr), domain(rhs.domain),
      nvar(rhs.nvar), truesize(rhs.truesize),
      ptr_owner(rhs.ptr_owner), shared_memory(rhs.shared_memory)
#ifdef AMREX_USE_GPU
      , alloc_stream(rhs.alloc_stream)
#endif
{
    rhs.dptr = nullptr;
    rhs.ptr_owner = false;
}
 
template <class T>
BaseFab<T>&
BaseFab<T>::operator= (BaseFab<T>&& rhs) noexcept
{
    if (this != &rhs) {
        clear();
        DataAllocator::operator=(rhs);
        dptr = rhs.dptr;
        domain = rhs.domain;
        nvar = rhs.nvar;
        truesize = rhs.truesize;
        ptr_owner = rhs.ptr_owner;
        shared_memory = rhs.shared_memory;
#ifdef AMREX_USE_GPU
        alloc_stream = rhs.alloc_stream;
#endif
 
        rhs.dptr = nullptr;
        rhs.ptr_owner = false;
    }
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator= (T const& t) noexcept
{
    setVal<run_on>(t);
    return *this;
}
 
template <class T>
void
BaseFab<T>::resize (const Box& b, int n, Arena* ar)
{
    this->nvar   = n;
    this->domain = b;
 
    if (ar == nullptr) {
        ar = m_arena;
    }
 
    if (arena() != DataAllocator(ar).arena()) {
        clear();
        m_arena = ar;
        define();
    }
    else if (this->dptr == nullptr || !this->ptr_owner)
    {
        if (this->shared_memory) {
            amrex::Abort("BaseFab::resize: BaseFab in shared memory cannot increase size");
        }
 
        this->dptr = nullptr;
        define();
    }
    else if (this->nvar*this->domain.numPts() > this->truesize
#ifdef AMREX_USE_GPU
             || (arena()->isStreamOrderedArena() && alloc_stream != Gpu::gpuStream())
#endif
             )
    {
        if (this->shared_memory) {
            amrex::Abort("BaseFab::resize: BaseFab in shared memory cannot increase size");
        }
 
        clear();
 
        define();
    }
}
 
template <class T>
template <class U, std::enable_if_t<std::is_trivially_destructible_v<U>,int>>
Elixir
BaseFab<T>::elixir () noexcept
{
    bool o;
    if (Gpu::inLaunchRegion()) {
        o = this->ptr_owner;
        this->ptr_owner = false;
        if (o && this->dptr) {
            if (this->nvar > 1) {
                amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
            } else {
                amrex::update_fab_stats(0, -this->truesize, sizeof(T));
            }
        }
    } else {
        o = false;
    }
    return Elixir((o ? this->dptr : nullptr), this->arena());
}
 
template <class T>
void
BaseFab<T>::clear () noexcept
{
    if (this->dptr)
    {
        //
        // Call T::~T() on the to-be-destroyed memory.
        //
        if (this->ptr_owner)
        {
            if (this->shared_memory)
            {
                amrex::Abort("BaseFab::clear: BaseFab cannot be owner of shared memory");
            }
 
            placementDelete(this->dptr, this->truesize);
 
#ifdef AMREX_USE_GPU
            auto current_stream = Gpu::Device::gpuStream();
            Gpu::Device::setStream(alloc_stream);
#endif
            this->free(this->dptr);
#ifdef AMREX_USE_GPU
            Gpu::Device::setStream(current_stream);
#endif
 
            if (this->nvar > 1) {
                amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
            } else {
                amrex::update_fab_stats(0, -this->truesize, sizeof(T));
            }
        }
 
        this->dptr = nullptr;
        this->truesize = 0;
    }
}
 
template <class T>
std::unique_ptr<T,DataDeleter>
BaseFab<T>::release () noexcept
{
    std::unique_ptr<T,DataDeleter> r(nullptr, DataDeleter{this->arena()});
    if (this->dptr && this->ptr_owner) {
        r.reset(this->dptr);
        this->ptr_owner = false;
        if (this->nvar > 1) {
            amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
        } else {
            amrex::update_fab_stats(0, -this->truesize, sizeof(T));
        }
    }
    return r;
}
 
template <class T>
template <RunOn run_on>
std::size_t
BaseFab<T>::copyToMem (const Box& srcbox,
                       int        srccomp,
                       int        numcomp,
                       void*      dst) const noexcept
{
    BL_ASSERT(box().contains(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= nComp());
 
    if (srcbox.ok())
    {
        Array4<T> d(static_cast<T*>(dst),amrex::begin(srcbox),amrex::end(srcbox),numcomp);
        Array4<T const> const& s = this->const_array();
        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, srcbox, numcomp, i, j, k, n,
        {
            d(i,j,k,n) = s(i,j,k,n+srccomp);
        });
        return sizeof(T)*d.size();
    }
    else
    {
        return 0;
    }
}
 
template <class T>
template <RunOn run_on, typename BUF>
std::size_t
BaseFab<T>::copyFromMem (const Box&  dstbox,
                         int         dstcomp,
                         int         numcomp,
                         const void* src) noexcept
{
    BL_ASSERT(box().contains(dstbox));
    BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
 
    if (dstbox.ok())
    {
        Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
                            amrex::end(dstbox), numcomp);
        Array4<T> const& d = this->array();
        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
        {
            d(i,j,k,n+dstcomp) = static_cast<T>(s(i,j,k,n));
        });
        return sizeof(BUF)*s.size();
    }
    else
    {
        return 0;
    }
}
 
template <class T>
template <RunOn run_on, typename BUF>
std::size_t
BaseFab<T>::addFromMem (const Box&  dstbox,
                        int         dstcomp,
                        int         numcomp,
                        const void* src) noexcept
{
    BL_ASSERT(box().contains(dstbox));
    BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
 
    if (dstbox.ok())
    {
        Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
                            amrex::end(dstbox), numcomp);
        Array4<T> const& d = this->array();
        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
        {
            d(i,j,k,n+dstcomp) += static_cast<T>(s(i,j,k,n));
        });
        return sizeof(BUF)*s.size();
    }
    else
    {
        return 0;
    }
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setComplement (T const& x, const Box& b, int ns, int num) noexcept
{
    this->setComplement<run_on>(x, b, DestComp{ns}, NumComps{num});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::abs () noexcept
{
    this->abs<run_on>(this->domain,0,this->nvar);
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::abs (int comp, int numcomp) noexcept
{
    this->abs<run_on>(this->domain,comp,numcomp);
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::abs (const Box& subbox, int comp, int numcomp) noexcept
{
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, subbox, numcomp, i, j, k, n,
    {
        a(i,j,k,n+comp) = std::abs(a(i,j,k,n+comp));
    });
}
 
template <class T>
template <RunOn run_on>
Real
BaseFab<T>::norminfmask (const Box& subbox, const BaseFab<int>& mask,
                         int scomp, int ncomp) const noexcept
{
    BL_ASSERT(this->domain.contains(subbox));
    BL_ASSERT(scomp >= 0 && scomp + ncomp <= this->nvar);
 
    Array4<T const> const& a = this->const_array();
    Array4<int const> const& m = mask.const_array();
    Real r = 0.0;
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpMax> reduce_op;
        ReduceData<Real> reduce_data(reduce_op);
        using ReduceTuple = ReduceData<Real>::Type;
        reduce_op.eval(subbox, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            Real t = 0.0;
            if (m(i,j,k)) {
                for (int n = 0; n < ncomp; ++n) {
                    t = amrex::max(t,static_cast<Real>(std::abs(a(i,j,k,n+scomp))));
                }
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        amrex::LoopOnCpu(subbox, ncomp, [=,&r] (int i, int j, int k, int n) noexcept
        {
            if (m(i,j,k)) {
                Real t = static_cast<Real>(std::abs(a(i,j,k,n+scomp)));
                r = amrex::max(r,t);
            }
        });
    }
    return r;
}
 
template <class T>
template <RunOn run_on>
Real
BaseFab<T>::norm (int p, int comp, int numcomp) const noexcept
{
    return norm<run_on>(this->domain,p,comp,numcomp);
}
 
template <class T>
template <RunOn run_on>
Real
BaseFab<T>::norm (const Box& subbox, int p, int comp, int numcomp) const noexcept
{
    BL_ASSERT(this->domain.contains(subbox));
    BL_ASSERT(comp >= 0 && comp + numcomp <= this->nvar);
 
    Array4<T const> const& a = this->const_array();
    Real nrm = 0.;
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        if (p == 0) {
            ReduceOps<ReduceOpMax> reduce_op;
            ReduceData<Real> reduce_data(reduce_op);
            using ReduceTuple = ReduceData<Real>::Type;
            reduce_op.eval(subbox, reduce_data,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
            {
                Real t = 0.0;
                for (int n = 0; n < numcomp; ++n) {
                    t = amrex::max(t, static_cast<Real>(std::abs(a(i,j,k,n+comp))));
                }
                return {t};
            });
            ReduceTuple hv = reduce_data.value(reduce_op);
            nrm = amrex::get<0>(hv);
        } else if (p == 1) {
            ReduceOps<ReduceOpSum> reduce_op;
            ReduceData<Real> reduce_data(reduce_op);
            using ReduceTuple = ReduceData<Real>::Type;
            reduce_op.eval(subbox, reduce_data,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
            {
                Real t = 0.0;
                for (int n = 0; n < numcomp; ++n) {
                    t += static_cast<Real>(std::abs(a(i,j,k,n+comp)));
                }
                return {t};
            });
            ReduceTuple hv = reduce_data.value(reduce_op);
            nrm = amrex::get<0>(hv);
        } else if (p == 2) {
            ReduceOps<ReduceOpSum> reduce_op;
            ReduceData<Real> reduce_data(reduce_op);
            using ReduceTuple = ReduceData<Real>::Type;
            reduce_op.eval(subbox, reduce_data,
            [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
            {
                Real t = 0.0;
                for (int n = 0; n < numcomp; ++n) {
                    t += static_cast<Real>(a(i,j,k,n+comp)*a(i,j,k,n+comp));
                }
                return {t};
            });
            ReduceTuple hv = reduce_data.value(reduce_op);
            nrm = amrex::get<0>(hv);
        } else {
            amrex::Error("BaseFab<T>::norm: wrong p");
        }
    } else
#endif
    {
        if (p == 0) {
            amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
            {
                Real t = static_cast<Real>(std::abs(a(i,j,k,n+comp)));
                nrm = amrex::max(nrm,t);
            });
        } else if (p == 1) {
            amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
            {
                nrm += std::abs(a(i,j,k,n+comp));
            });
        } else if (p == 2) {
            amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
            {
                nrm += a(i,j,k,n+comp)*a(i,j,k,n+comp);
            });
        } else {
            amrex::Error("BaseFab<T>::norm: wrong p");
        }
    }
 
    return nrm;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::min (int comp) const noexcept
{
    return this->min<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::min (const Box& subbox, int comp) const noexcept
{
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpMin> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(subbox, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            return { a(i,j,k) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        T r = std::numeric_limits<T>::max();
        amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
        {
            r = amrex::min(r, a(i,j,k));
        });
        return r;
    }
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::max (int comp) const noexcept
{
    return this->max<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::max (const Box& subbox, int comp) const noexcept
{
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpMax> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(subbox, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            return { a(i,j,k) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        T r = std::numeric_limits<T>::lowest();
        amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
        {
            r = amrex::max(r, a(i,j,k));
        });
        return r;
    }
}
 
template <class T>
template <RunOn run_on>
std::pair<T,T>
BaseFab<T>::minmax (int comp) const noexcept
{
    return this->minmax<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
std::pair<T,T>
BaseFab<T>::minmax (const Box& subbox, int comp) const noexcept
{
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpMin,ReduceOpMax> reduce_op;
        ReduceData<T,T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(subbox, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            auto const x = a(i,j,k);
            return { x, x };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return std::make_pair(amrex::get<0>(hv), amrex::get<1>(hv));
    } else
#endif
    {
        T rmax = std::numeric_limits<T>::lowest();
        T rmin = std::numeric_limits<T>::max();
        amrex::LoopOnCpu(subbox, [=,&rmin,&rmax] (int i, int j, int k) noexcept
        {
            auto const x = a(i,j,k);
            rmin = amrex::min(rmin, x);
            rmax = amrex::max(rmax, x);
        });
        return std::make_pair(rmin,rmax);
    }
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::maxabs (int comp) const noexcept
{
    return this->maxabs<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::maxabs (const Box& subbox, int comp) const noexcept
{
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpMax> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(subbox, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            return { std::abs(a(i,j,k)) };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        return amrex::get<0>(hv);
    } else
#endif
    {
        T r = 0;
        amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
        {
            r = amrex::max(r, std::abs(a(i,j,k)));
        });
        return r;
    }
}
 
 
template <class T>
template <RunOn run_on>
IntVect
BaseFab<T>::indexFromValue (Box const& subbox, int comp, T const& value) const noexcept
{
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        Array<int,1+AMREX_SPACEDIM> ha{0,AMREX_D_DECL(std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest(),
                                                      std::numeric_limits<int>::lowest())};
        Gpu::DeviceVector<int> dv(1+AMREX_SPACEDIM);
        int* p = dv.data();
        Gpu::htod_memcpy_async(p, ha.data(), sizeof(int)*ha.size());
        amrex::ParallelFor(subbox, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            int* flag = p;
            if ((*flag == 0) && (a(i,j,k) == value)) {
                if (Gpu::Atomic::Exch(flag,1) == 0) {
                    AMREX_D_TERM(p[1] = i;,
                                 p[2] = j;,
                                 p[3] = k;);
                }
            }
        });
        Gpu::dtoh_memcpy_async(ha.data(), p, sizeof(int)*ha.size());
        Gpu::streamSynchronize();
        return IntVect(AMREX_D_DECL(ha[1],ha[2],ha[2]));
    } else
#endif
    {
        AMREX_LOOP_3D(subbox, i, j, k,
        {
            if (a(i,j,k) == value) { return IntVect(AMREX_D_DECL(i,j,k)); }
        });
        return IntVect::TheMinVector();
    }
}
 
template <class T>
template <RunOn run_on>
IntVect
BaseFab<T>::minIndex (int comp) const noexcept
{
    return this->minIndex<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
IntVect
BaseFab<T>::minIndex (const Box& subbox, int comp) const noexcept
{
    T min_val = this->min<run_on>(subbox, comp);
    return this->indexFromValue<run_on>(subbox, comp, min_val);
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::minIndex (const Box& subbox, Real& min_val, IntVect& min_idx, int comp) const noexcept
{
    min_val = this->min<run_on>(subbox, comp);
    min_idx = this->indexFromValue<run_on>(subbox, comp, min_val);
}
 
template <class T>
template <RunOn run_on>
IntVect
BaseFab<T>::maxIndex (int comp) const noexcept
{
    return this->maxIndex<run_on>(this->domain,comp);
}
 
template <class T>
template <RunOn run_on>
IntVect
BaseFab<T>::maxIndex (const Box& subbox, int comp) const noexcept
{
    T max_val = this->max<run_on>(subbox, comp);
    return this->indexFromValue<run_on>(subbox, comp, max_val);
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::maxIndex (const Box& subbox, Real& max_val, IntVect& max_idx, int comp) const noexcept
{
    max_val = this->max<run_on>(subbox, comp);
    max_idx = this->indexFromValue<run_on>(subbox, comp, max_val);
}
 
template <class T>
template <RunOn run_on>
int
BaseFab<T>::maskLT (BaseFab<int>& mask, T const& val, int comp) const noexcept
{
    mask.resize(this->domain,1);
    int cnt = 0;
    Array4<int> const& m = mask.array();
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(this->domain, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int t;
            if (a(i,j,k) < val) {
                m(i,j,k) = 1;
                t = 1;
            } else {
                m(i,j,k) = 0;
                t = 0;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        cnt = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_3D(this->domain, i, j, k,
        {
            if (a(i,j,k) < val) {
                m(i,j,k) = 1;
                ++cnt;
            } else {
                m(i,j,k) = 0;
            }
        });
    }
 
    return cnt;
}
 
template <class T>
template <RunOn run_on>
int
BaseFab<T>::maskLE (BaseFab<int>& mask, T const& val, int comp) const noexcept
{
    mask.resize(this->domain,1);
    int cnt = 0;
    Array4<int> const& m = mask.array();
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(this->domain, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int t;
            if (a(i,j,k) <= val) {
                m(i,j,k) = 1;
                t = 1;
            } else {
                m(i,j,k) = 0;
                t = 0;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        cnt = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_3D(this->domain, i, j, k,
        {
            if (a(i,j,k) <= val) {
                m(i,j,k) = 1;
                ++cnt;
            } else {
                m(i,j,k) = 0;
            }
        });
    }
 
    return cnt;
}
 
template <class T>
template <RunOn run_on>
int
BaseFab<T>::maskEQ (BaseFab<int>& mask, T const& val, int comp) const noexcept
{
    mask.resize(this->domain,1);
    int cnt = 0;
    Array4<int> const& m = mask.array();
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(this->domain, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int t;
            if (a(i,j,k) == val) {
                m(i,j,k) = 1;
                t = 1;
            } else {
                m(i,j,k) = 0;
                t = 0;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        cnt = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_3D(this->domain, i, j, k,
        {
            if (a(i,j,k) == val) {
                m(i,j,k) = 1;
                ++cnt;
            } else {
                m(i,j,k) = 0;
            }
        });
    }
 
    return cnt;
}
 
template <class T>
template <RunOn run_on>
int
BaseFab<T>::maskGT (BaseFab<int>& mask, T const& val, int comp) const noexcept
{
    mask.resize(this->domain,1);
    int cnt = 0;
    Array4<int> const& m = mask.array();
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(this->domain, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int t;
            if (a(i,j,k) > val) {
                m(i,j,k) = 1;
                t = 1;
            } else {
                m(i,j,k) = 0;
                t = 0;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        cnt = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_3D(this->domain, i, j, k,
        {
            if (a(i,j,k) > val) {
                m(i,j,k) = 1;
                ++cnt;
            } else {
                m(i,j,k) = 0;
            }
        });
    }
 
    return cnt;
}
 
template <class T>
template <RunOn run_on>
int
BaseFab<T>::maskGE (BaseFab<int>& mask, T const& val, int comp) const noexcept
{
    mask.resize(this->domain,1);
    int cnt = 0;
    Array4<int> const& m = mask.array();
    Array4<T const> const& a = this->const_array(comp);
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<int> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(this->domain, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int t;
            if (a(i,j,k) >= val) {
                m(i,j,k) = 1;
                t = 1;
            } else {
                m(i,j,k) = 0;
                t = 0;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        cnt = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_3D(this->domain, i, j, k,
        {
            if (a(i,j,k) >= val) {
                m(i,j,k) = 1;
                ++cnt;
            } else {
                m(i,j,k) = 0;
            }
        });
    }
 
    return cnt;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::atomicAdd (const BaseFab<T>& x) noexcept
{
    Box ovlp(this->domain);
    ovlp &= x.domain;
    return ovlp.ok() ? this->atomicAdd<run_on>(x,ovlp,ovlp,0,0,this->nvar) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::saxpy (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
                   int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(srcbox.ok());
    BL_ASSERT(x.box().contains(srcbox));
    BL_ASSERT(destbox.ok());
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT( srccomp >= 0 &&  srccomp+numcomp <= x.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <=   nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = x.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) += a * s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::saxpy (T a, const BaseFab<T>& x) noexcept
{
    Box ovlp(this->domain);
    ovlp &= x.domain;
    return ovlp.ok() ? saxpy<run_on>(a,x,ovlp,ovlp,0,0,this->nvar) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::xpay (T a, const BaseFab<T>& x,
                  const Box& srcbox, const Box& destbox,
                  int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(srcbox.ok());
    BL_ASSERT(x.box().contains(srcbox));
    BL_ASSERT(destbox.ok());
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT( srccomp >= 0 &&  srccomp+numcomp <= x.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <=   nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = x.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) = s(i+offset.x,j+offset.y,k+offset.z,n+srccomp) + a*d(i,j,k,n+destcomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::addproduct (const Box& destbox, int destcomp, int numcomp,
                        const BaseFab<T>& src1, int comp1,
                        const BaseFab<T>& src2, int comp2) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(   comp1 >= 0 &&    comp1+numcomp <= src1.nComp());
    BL_ASSERT(   comp2 >= 0 &&    comp2+numcomp <= src2.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <=      nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s1 = src1.const_array();
    Array4<T const> const& s2 = src2.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) += s1(i,j,k,n+comp1) * s2(i,j,k,n+comp2);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::linComb (const BaseFab<T>& f1, const Box& b1, int comp1,
                     const BaseFab<T>& f2, const Box& b2, int comp2,
                     Real alpha, Real beta, const Box& b,
                     int comp, int numcomp) noexcept
{
    BL_ASSERT(b1.ok());
    BL_ASSERT(f1.box().contains(b1));
    BL_ASSERT(b2.ok());
    BL_ASSERT(f2.box().contains(b2));
    BL_ASSERT(b.ok());
    BL_ASSERT(box().contains(b));
    BL_ASSERT(b.sameSize(b1));
    BL_ASSERT(b.sameSize(b2));
    BL_ASSERT(comp1 >= 0 && comp1+numcomp <= f1.nComp());
    BL_ASSERT(comp2 >= 0 && comp2+numcomp <= f2.nComp());
    BL_ASSERT(comp  >= 0 && comp +numcomp <=    nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s1 = f1.const_array();
    Array4<T const> const& s2 = f2.const_array();
    const auto dlo = amrex::lbound(b);
    const auto slo1 = amrex::lbound(b1);
    const auto slo2 = amrex::lbound(b2);
    const Dim3 off1{slo1.x-dlo.x,slo1.y-dlo.y,slo1.z-dlo.z};
    const Dim3 off2{slo2.x-dlo.x,slo2.y-dlo.y,slo2.z-dlo.z};
 
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, b, numcomp, i, j, k, n,
    {
        d(i,j,k,n+comp) = alpha*s1(i+off1.x,j+off1.y,k+off1.z,n+comp1)
                         + beta*s2(i+off2.x,j+off2.y,k+off2.z,n+comp2);
    });
    return *this;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dot (const Box& xbx, int xcomp,
                 const BaseFab<T>& y, const Box& ybx, int ycomp,
                 int numcomp) const noexcept
{
    BL_ASSERT(xbx.ok());
    BL_ASSERT(box().contains(xbx));
    BL_ASSERT(y.box().contains(ybx));
    BL_ASSERT(xbx.sameSize(ybx));
    BL_ASSERT(xcomp >= 0 && xcomp+numcomp <=   nComp());
    BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
 
    T r = 0;
 
    const auto xlo = amrex::lbound(xbx);
    const auto ylo = amrex::lbound(ybx);
    const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
    Array4<T const> const& xa = this->const_array();
    Array4<T const> const& ya = y.const_array();
 
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(xbx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            T t = 0;
            for (int n = 0; n < numcomp; ++n) {
                t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
        {
            r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
        });
    }
 
    return r;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dotmask (const BaseFab<int>& mask, const Box& xbx, int xcomp,
                     const BaseFab<T>& y, const Box& ybx, int ycomp,
                     int numcomp) const noexcept
{
    BL_ASSERT(xbx.ok());
    BL_ASSERT(box().contains(xbx));
    BL_ASSERT(y.box().contains(ybx));
    BL_ASSERT(xbx.sameSize(ybx));
    BL_ASSERT(xcomp >= 0 && xcomp+numcomp <=   nComp());
    BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
 
    T r = 0;
 
    const auto xlo = amrex::lbound(xbx);
    const auto ylo = amrex::lbound(ybx);
    const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
 
    Array4<T const> const& xa = this->const_array();
    Array4<T const> const& ya = y.const_array();
    Array4<int const> const& ma = mask.const_array();
 
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(xbx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
            T t = 0;
            for (int n = 0; n < numcomp; ++n) {
                t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
            }
            return {t};
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
        {
            int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
            r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
        });
    }
 
    return r;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::sum (int comp, int numcomp) const noexcept
{
    return this->sum<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::sum (const Box& subbox, int comp, int numcomp) const noexcept
{
    return this->sum<run_on>(subbox, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::negate (int comp, int numcomp) noexcept
{
    return this->negate<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::negate (const Box& b, int comp, int numcomp) noexcept
{
    return this->negate<run_on>(b, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::invert (T const& r, int comp, int numcomp) noexcept
{
    return this->invert<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::invert (T const& r, const Box& b, int comp, int numcomp) noexcept
{
    return this->invert<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (T const& r, int comp, int numcomp) noexcept
{
    return this->plus<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (T const& r, const Box& b, int comp, int numcomp) noexcept
{
    return this->plus<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->plus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::atomicAdd (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    Box ovlp(this->domain);
    ovlp &= src.domain;
    return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                  int numcomp) noexcept
{
    return this->plus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                       int numcomp) noexcept
{
    Box ovlp(this->domain);
    ovlp &= src.domain;
    ovlp &= subbox;
    return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                  int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) += s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                       int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        T* p = d.ptr(i,j,k,n+destcomp);
        HostDevice::Atomic::Add(p, s(i+offset.x,j+offset.y,k+offset.z,n+srccomp));
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::lockAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                     int srccomp, int destcomp, int numcomp) noexcept
{
#if defined(AMREX_USE_OMP) && (AMREX_SPACEDIM > 1)
#if defined(AMREX_USE_GPU)
    if (run_on == RunOn::Host || Gpu::notInLaunchRegion()) {
#endif
        BL_ASSERT(destbox.ok());
        BL_ASSERT(src.box().contains(srcbox));
        BL_ASSERT(box().contains(destbox));
        BL_ASSERT(destbox.sameSize(srcbox));
        BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
        BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
        Array4<T> const& d = this->array();
        Array4<T const> const& s = src.const_array();
        auto const& dlo = amrex::lbound(destbox);
        auto const& dhi = amrex::ubound(destbox);
        auto const& len = amrex::length(destbox);
        auto const& slo = amrex::lbound(srcbox);
        Dim3 const offset{slo.x-dlo.x, slo.y-dlo.y, slo.z-dlo.z};
 
        int planedim;
        int nplanes;
        int plo;
        if (len.z == 1) {
            planedim = 1;
            nplanes = len.y;
            plo = dlo.y;
        } else {
            planedim = 2;
            nplanes = len.z;
            plo = dlo.z;
        }
 
        auto* mask = (bool*) amrex_mempool_alloc(sizeof(bool)*nplanes);
        for (int ip = 0; ip < nplanes; ++ip) {
            mask[ip] = false;
        }
 
        int mm = 0;
        int planes_left = nplanes;
        while (planes_left > 0) {
            AMREX_ASSERT(mm < nplanes);
            auto const m = mm + plo;
            auto* lock = OpenMP::get_lock(m);
            if (omp_test_lock(lock))
            {
                auto lo = dlo;
                auto hi = dhi;
                if (planedim == 1) {
                    lo.y = m;
                    hi.y = m;
                } else {
                    lo.z = m;
                    hi.z = m;
                }
 
                for (int n = 0; n < numcomp; ++n) {
                    for     (int k = lo.z; k <= hi.z; ++k) {
                        for (int j = lo.y; j <= hi.y; ++j) {
                            auto      * pdst = d.ptr(dlo.x,j         ,k         ,n+destcomp);
                            auto const* psrc = s.ptr(slo.x,j+offset.y,k+offset.z,n+ srccomp);
#pragma omp simd
                            for (int ii = 0; ii < len.x; ++ii) {
                                pdst[ii] += psrc[ii];
                            }
                        }
                    }
                }
 
                mask[mm] = true;
                --planes_left;
                omp_unset_lock(lock);
                if (planes_left == 0) { break; }
            }
 
            ++mm;
            for (int ip = 0; ip < nplanes; ++ip) {
                int new_mm = (mm+ip) % nplanes;
                if ( ! mask[new_mm] ) {
                    mm = new_mm;
                    break;
                }
            }
        }
 
        amrex_mempool_free(mask);
 
        return *this;
 
#if defined(AMREX_USE_GPU)
    } else {
        return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
    }
#endif
#else
    return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
#endif
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->minus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->minus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                   int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) -= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (T const& r, int comp, int numcomp) noexcept
{
    return this->mult<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (T const& r, const Box& b, int comp, int numcomp) noexcept
{
    return this->mult<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->mult<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->mult<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                  int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) *= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (T const& r, int comp, int numcomp) noexcept
{
    return this->divide<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (T const& r, const Box& b, int comp, int numcomp) noexcept
{
    return this->divide<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->divide<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
{
    return this->divide<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                    int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::protected_divide (const BaseFab<T>& src) noexcept
{
    Box ovlp(this->domain);
    ovlp &= src.domain;
    return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,0,0,this->nvar) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::protected_divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
{
    Box ovlp(this->domain);
    ovlp &= src.domain;
    return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
                              int numcomp) noexcept
{
    Box ovlp(this->domain);
    ovlp &= src.domain;
    ovlp &= subbox;
    return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
                              int srccomp, int destcomp, int numcomp) noexcept
{
    BL_ASSERT(destbox.ok());
    BL_ASSERT(src.box().contains(srcbox));
    BL_ASSERT(box().contains(destbox));
    BL_ASSERT(destbox.sameSize(srcbox));
    BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
    BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    const auto dlo = amrex::lbound(destbox);
    const auto slo = amrex::lbound(srcbox);
    const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
    {
        if (s(i+offset.x,j+offset.y,k+offset.z,n+srccomp) != 0) {
            d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
        }
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::linInterp (const BaseFab<T>& f1, const Box& b1, int comp1,
                       const BaseFab<T>& f2, const Box& b2, int comp2,
                       Real t1, Real t2, Real t,
                       const Box& b, int comp, int numcomp) noexcept
{
    if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
        return copy<run_on>(f1,b1,comp1,b,comp,numcomp);
    } else if (amrex::almostEqual(t2,t)) {
        return copy<run_on>(f2,b2,comp2,b,comp,numcomp);
    } else {
        Real alpha = (t2-t)/(t2-t1);
        Real beta = (t-t1)/(t2-t1);
        return linComb<run_on>(f1,b1,comp1,f2,b2,comp2,alpha,beta,b,comp,numcomp);
    }
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::linInterp (const BaseFab<T>& f1, int comp1,
                       const BaseFab<T>& f2, int comp2,
                       Real t1, Real t2, Real t,
                       const Box& b, int comp, int numcomp) noexcept
{
    if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
        return copy<run_on>(f1,b,comp1,b,comp,numcomp);
    } else if (amrex::almostEqual(t2,t)) {
        return copy<run_on>(f2,b,comp2,b,comp,numcomp);
    } else {
        Real alpha = (t2-t)/(t2-t1);
        Real beta = (t-t1)/(t2-t1);
        return linComb<run_on>(f1,b,comp1,f2,b,comp2,alpha,beta,b,comp,numcomp);
    }
}
 
//
// New interfaces
//
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setVal (T const& val) noexcept
{
    this->setVal<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setVal (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) = x;
    });
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setValIf (T const& val, const BaseFab<int>& mask) noexcept
{
    this->setValIf<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setValIf (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
    Array4<T> const& a = this->array();
    Array4<int const> const& m = mask.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
    {
        if (m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
    });
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setValIfNot (T const& val, const BaseFab<int>& mask) noexcept
{
    this->setValIfNot<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setValIfNot (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
    Array4<T> const& a = this->array();
    Array4<int const> const& m = mask.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
    {
        if (!m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
    });
}
 
template <class T>
template <RunOn run_on>
void
BaseFab<T>::setComplement (T const& x, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
{
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        Array4<T> const& a = this->array();
        amrex::ParallelFor(this->domain, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            if (! bx.contains(IntVect(AMREX_D_DECL(i,j,k)))) {
                for (int n = dcomp.i; n < ncomp.n+dcomp.i; ++n) {
                    a(i,j,k,n) = x;
                }
            }
        });
    } else
#endif
    {
        const BoxList b_lst = amrex::boxDiff(this->domain,bx);
        for (auto const& b : b_lst) {
            this->setVal<RunOn::Host>(x, b, dcomp, ncomp);
        }
    }
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src) noexcept
{
    this->copy<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::copy (const BaseFab<T>& src, Box bx,
                  SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    bx &= src.domain;
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        d(i,j,k,n+dcomp.i) = s(i,j,k,n+scomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (T const& val) noexcept
{
    return this->plus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator+= (T const& val) noexcept
{
    return this->plus<run_on>(val, this->domain,  DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) += val;
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (const BaseFab<T>& src) noexcept
{
    return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator+= (const BaseFab<T>& src) noexcept
{
    return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::plus (const BaseFab<T>& src, Box bx,
                  SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    bx &= src.domain;
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        d(i,j,k,n+dcomp.i) += s(i,j,k,n+scomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (T const& val) noexcept
{
    return this->minus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator-= (T const& val) noexcept
{
    return this->minus<run_on>(val, this->domain,  DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) -= val;
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (const BaseFab<T>& src) noexcept
{
    return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator-= (const BaseFab<T>& src) noexcept
{
    return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::minus (const BaseFab<T>& src, Box bx,
                   SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    bx &= src.domain;
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        d(i,j,k,n+dcomp.i) -= s(i,j,k,n+scomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (T const& val) noexcept
{
    return this->mult<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator*= (T const& val) noexcept
{
    return this->mult<run_on>(val, this->domain,  DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) *= val;
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (const BaseFab<T>& src) noexcept
{
    return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator*= (const BaseFab<T>& src) noexcept
{
    return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::mult (const BaseFab<T>& src, Box bx,
                  SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    bx &= src.domain;
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        d(i,j,k,n+dcomp.i) *= s(i,j,k,n+scomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (T const& val) noexcept
{
    return this->divide<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator/= (T const& val) noexcept
{
    return this->divide<run_on>(val, this->domain,  DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) /= val;
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (const BaseFab<T>& src) noexcept
{
    return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::operator/= (const BaseFab<T>& src) noexcept
{
    return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::divide (const BaseFab<T>& src, Box bx,
                    SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    bx &= src.domain;
 
    Array4<T> const& d = this->array();
    Array4<T const> const& s = src.const_array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        d(i,j,k,n+dcomp.i) /= s(i,j,k,n+scomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::negate () noexcept
{
    return this->negate<run_on>(this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::negate (const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) = -a(i,j,k,n+dcomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::invert (T const& r) noexcept
{
    return this->invert<run_on>(r, this->domain, DestComp{0}, NumComps{this->nvar});
}
 
template <class T>
template <RunOn run_on>
BaseFab<T>&
BaseFab<T>::invert (T const& r, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
{
    BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
 
    Array4<T> const& a = this->array();
    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
    {
        a(i,j,k,n+dcomp.i) = r / a(i,j,k,n+dcomp.i);
    });
 
    return *this;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::sum (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
{
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    T r = 0;
    Array4<T const> const& a = this->const_array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            T t = 0;
            for (int n = 0; n < ncomp.n; ++n) {
                t += a(i,j,k,n+dcomp.i);
            }
            return { t };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
        {
            r += a(i,j,k,n+dcomp.i);
        });
    }
 
    return r;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dot (const BaseFab<T>& src, const Box& bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    T r = 0;
    Array4<T const> const& d = this->const_array();
    Array4<T const> const& s = src.const_array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            T t = 0;
            for (int n = 0; n < ncomp.n; ++n) {
                t += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
            }
            return { t };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
        {
            r += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
        });
    }
 
    return r;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dot (const Box& bx, int destcomp, int numcomp) const noexcept
{
    return dot<run_on>(bx, DestComp{destcomp}, NumComps{numcomp});
}
 
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dot (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
{
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    T r = 0;
    Array4<T const> const& a = this->const_array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            T t = 0;
            for (int n = 0; n < ncomp.n; ++n) {
                t += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
            }
            return { t };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
        {
            r += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
        });
    }
 
    return r;
}
 
template <class T>
template <RunOn run_on>
T
BaseFab<T>::dotmask (const BaseFab<T>& src, const Box& bx, const BaseFab<int>& mask,
                     SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
{
    AMREX_ASSERT(this->domain.sameType(src.domain));
    AMREX_ASSERT(this->domain.sameType(mask.domain));
    AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
    AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
 
    T r = 0;
    Array4<T const> const& d = this->const_array();
    Array4<T const> const& s = src.const_array();
    Array4<int const> const& m = mask.const_array();
#ifdef AMREX_USE_GPU
    if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
        ReduceOps<ReduceOpSum> reduce_op;
        ReduceData<T> reduce_data(reduce_op);
        using ReduceTuple = typename decltype(reduce_data)::Type;
        reduce_op.eval(bx, reduce_data,
        [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
        {
            T t = 0;
            T mi = static_cast<T>(static_cast<int>(static_cast<bool>(m(i,j,k))));
            for (int n = 0; n < ncomp.n; ++n) {
                t += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
            }
            return { t };
        });
        ReduceTuple hv = reduce_data.value(reduce_op);
        r = amrex::get<0>(hv);
    } else
#endif
    {
        amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
        {
            int mi = static_cast<int>(static_cast<bool>(m(i,j,k)));
            r += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
        });
    }
 
    return r;
}
 
}
 
#endif /*BL_BASEFAB_H*/