AMReX_FFT_Helper.H

Go to the documentation of this file.

#ifndef AMREX_FFT_HELPER_H_
#define AMREX_FFT_HELPER_H_
#include <AMReX_Config.H>
 
#include <AMReX.H>
#include <AMReX_BLProfiler.H>
#include <AMReX_DataAllocator.H>
#include <AMReX_DistributionMapping.H>
#include <AMReX_Enum.H>
#include <AMReX_FabArray.H>
#include <AMReX_Gpu.H>
#include <AMReX_GpuComplex.H>
#include <AMReX_Geometry.H>
#include <AMReX_Math.H>
#include <AMReX_NonLocalBC.H>
#include <AMReX_Periodicity.H>

 
#if defined(AMREX_USE_CUDA)
#  include <cufft.h>
#  include <cuComplex.h>
#elif defined(AMREX_USE_HIP)
#  if __has_include(<rocfft/rocfft.h>)  // ROCm 5.3+
#    include <rocfft/rocfft.h>
#  else
#    include <rocfft.h>
#  endif
#  include <hip/hip_complex.h>
#elif defined(AMREX_USE_SYCL)
#  if __has_include(<oneapi/mkl/dft.hpp>) // oneAPI 2025.0
#    include <oneapi/mkl/dft.hpp>
#else
#    define AMREX_USE_MKL_DFTI_2024 1
#    include <oneapi/mkl/dfti.hpp>
#  endif
#else
#  include <fftw3.h>
#endif
 
#include <algorithm>
#include <complex>
#include <initializer_list>
#include <limits>
#include <memory>
#include <tuple>
#include <utility>
#include <variant>
 
namespace amrex::FFT
{
 
enum struct Direction { forward, backward, both, none };
 
enum struct DomainStrategy { automatic, slab, pencil };

AMREX_ENUM( CommStrategy, automatic, pointtopoint, alltoallv, alltoall );
 
AMREX_ENUM( Boundary, periodic, even, odd );
 
enum struct Kind { none, r2c_f, r2c_b, c2c_f, c2c_b, r2r_ee_f, r2r_ee_b,
                   r2r_oo_f, r2r_oo_b, r2r_eo, r2r_oe };

[[nodiscard]] constexpr int
FastNumPrimeFactors () noexcept
{
#if defined(AMREX_USE_CUDA)
    return 5;
#else
    return 6;
#endif
}
 
struct Info
{
    DomainStrategy domain_strategy = DomainStrategy::automatic;

    int pencil_threshold = 4;

    bool twod_mode = false;

    bool oned_mode = false;

    int batch_size = 1;

    int nprocs = std::numeric_limits<int>::max();

    CommStrategy comm_strategy = CommStrategy::automatic;

    bool openbc_padding = true;

    int openbc_padding_nfactors = FastNumPrimeFactors();

    Info& setDomainStrategy (DomainStrategy s) { domain_strategy = s; return *this; }
    Info& setPencilThreshold (int t) { pencil_threshold = t; return *this; }
    Info& setTwoDMode (bool x) { twod_mode = x; return *this; }
    Info& setOneDMode (bool x) { oned_mode = x; return *this; }
    Info& setBatchSize (int bsize) { batch_size = bsize; return *this; }
    Info& setNumProcs (int n) { nprocs = n; return *this; }
    Info& setCommStrategy (CommStrategy s) { comm_strategy = s; return *this; }
    Info& setOpenBCPadding (bool x) { openbc_padding = x; return *this; }
    Info& setOpenBCPaddingNumPrimeFactors (int n)
    {
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            3 <= n && 6 >= n,
            "amrex::FFT::Info::setOpenBCPaddingNumPrimeFactors: nfactors must be between 3 and 6");
        openbc_padding_nfactors = n;
        return *this;
    }
};

namespace detail
{
inline Geometry shift_geom (Geometry const& geom)
{
    auto const& dom = geom.Domain();
    auto const& dlo = dom.smallEnd();
    if (dlo == 0) {
        return geom;
    } else {
        return Geometry(amrex::shift(dom,-dlo), geom.ProbDomain(), geom.CoordInt(),
                        geom.isPeriodic());
    }
}
 
template <typename MF>
void shift_mf (IntVect const& domain_lo, MF const& mf, MF& new_mf)
{
    BoxArray ba = mf.boxArray();
    ba.shift(-domain_lo);
    new_mf.define(ba, mf.DistributionMap(), mf.nComp(), mf.nGrowVect(),
                  MFInfo().SetAlloc(false));
    for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
        typename MF::fab_type fab(mf[mfi], amrex::make_alias, 0, mf.nComp());
        fab.shift(-domain_lo);
        new_mf.setFab(mfi, std::move(fab));
    }
}
 
template <typename MF>
void shift_mfs (IntVect const& domain_lo, MF const& mf1, MF const& mf2,
                MF& new_mf1, MF& new_mf2)
{
    AMREX_ALWAYS_ASSERT(mf1.boxArray() == mf2.boxArray() &&
                        mf1.DistributionMap() == mf2.DistributionMap());
    detail::shift_mf(domain_lo, mf1, new_mf1);
    detail::shift_mf(domain_lo, mf2, new_mf2);
}
 
inline void next_fast_len_doit (Long target, Long n, int const* factors,
                                int nfactors, int ifactor, Long& best)
{
    if (n >= target) {
        best = std::min(best, n);
        return;
    }
 
    if (ifactor == nfactors) {
        // Factors 2 and 3 are handled here: try each n * 3^b, multiply it
        // by the smallest power of 2 that reaches target, and keep the best.
        Long const max_m3 = std::numeric_limits<int>::max() / 3;
        Long const max_m2 = std::numeric_limits<int>::max() / 2;
        for (Long m3 = n; m3 < best; ) {
            Long k = m3;
            bool overflow = false;
            while (k < target) {
                if (k > max_m2) {
                    overflow = true;
                    break;
                }
                k *= 2;
            }
            if (!overflow && k < best) {
                best = k;
            }
            if (m3 > max_m3) {
                break;
            }
            m3 *= 3;
        }
        return;
    }
 
    int const factor = factors[ifactor];
    Long const max_len = std::numeric_limits<int>::max();
    Long const max_m = max_len / factor;
    for (Long m = n; m < best; ) {
        next_fast_len_doit(target, m, factors, nfactors, ifactor+1, best);
        if (m > max_m) {
            break;
        }
        m *= factor;
    }
}
}

[[nodiscard]] inline int
nextFastLen (int target, int nfactors = FastNumPrimeFactors())
{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(target > 0,
                                     "amrex::FFT::nextFastLen: target must be positive");
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(3 <= nfactors && 6 >= nfactors,
                                     "amrex::FFT::nextFastLen: nfactors must be between 3 and 6");
 
    int constexpr factors[] = {5, 7, 11, 13, 17}; // The extra '17' is to make the buggy GCC -Warray-bounds quiet.
 
    Long const max_len = std::numeric_limits<int>::max();
    Long best = max_len + Long(1);
 
    // Factors 2 and 3 are handled directly in the base case.
    detail::next_fast_len_doit(target, Long(1), factors, nfactors-2, 0, best);
 
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(best <= max_len,
                                     "amrex::FFT::nextFastLen: result exceeds int range");
    return static_cast<int>(best);
}
 
#ifdef AMREX_USE_HIP
namespace detail { void hip_execute (rocfft_plan plan, void **in, void **out); }
#endif
 
#ifdef AMREX_USE_SYCL
namespace detail
{
inline void assert_no_external_stream ()
{
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        !amrex::Gpu::Device::usingExternalStream(),
        "SYCL FFT does not support external GPU streams.");
}
 
template <typename T, Direction direction, typename P, typename TI, typename TO>
void sycl_execute (P* plan, TI* in, TO* out)
{
    assert_no_external_stream();
#ifndef AMREX_USE_MKL_DFTI_2024
    std::int64_t workspaceSize = 0;
#else
    std::size_t workspaceSize = 0;
#endif
    plan->get_value(oneapi::mkl::dft::config_param::WORKSPACE_BYTES,
                    &workspaceSize);
    auto* buffer = (T*)amrex::The_Arena()->alloc(workspaceSize);
    plan->set_workspace(buffer);
    sycl::event r;
    if (std::is_same_v<TI,TO>) {
        amrex::ignore_unused(in);
        if constexpr (direction == Direction::forward) {
            r = oneapi::mkl::dft::compute_forward(*plan, out);
        } else {
            r = oneapi::mkl::dft::compute_backward(*plan, out);
        }
    } else {
        if constexpr (direction == Direction::forward) {
            r = oneapi::mkl::dft::compute_forward(*plan, in, out);
        } else {
            r = oneapi::mkl::dft::compute_backward(*plan, in, out);
        }
    }
    r.wait();
    amrex::The_Arena()->free(buffer);
}
}
#endif
 
template <typename T>
struct Plan
{
#if defined(AMREX_USE_CUDA)
    using VendorPlan = cufftHandle;
    using VendorComplex = std::conditional_t<std::is_same_v<float,T>,
                                             cuComplex, cuDoubleComplex>;
#elif defined(AMREX_USE_HIP)
    using VendorPlan = rocfft_plan;
    using VendorComplex = std::conditional_t<std::is_same_v<float,T>,
                                             float2, double2>;
#elif defined(AMREX_USE_SYCL)
    using mkl_desc_r = oneapi::mkl::dft::descriptor<std::is_same_v<float,T>
                                     ? oneapi::mkl::dft::precision::SINGLE
                                     : oneapi::mkl::dft::precision::DOUBLE,
                                     oneapi::mkl::dft::domain::REAL>;
    using mkl_desc_c = oneapi::mkl::dft::descriptor<std::is_same_v<float,T>
                                     ? oneapi::mkl::dft::precision::SINGLE
                                     : oneapi::mkl::dft::precision::DOUBLE,
                                     oneapi::mkl::dft::domain::COMPLEX>;
    using VendorPlan = std::variant<mkl_desc_r*,mkl_desc_c*>;
    using VendorComplex = std::complex<T>;
#else
    using VendorPlan = std::conditional_t<std::is_same_v<float,T>,
                                          fftwf_plan, fftw_plan>;
    using VendorComplex = std::conditional_t<std::is_same_v<float,T>,
                                             fftwf_complex, fftw_complex>;
#endif
 
    int n = 0;
    int howmany = 0;
    Kind kind = Kind::none;
    bool r2r_data_is_complex = false;
    bool defined = false;
    bool defined2 = false;
    VendorPlan plan{};
    VendorPlan plan2{};
    void* pf = nullptr;
    void* pb = nullptr;
 
#ifdef AMREX_USE_GPU
    void set_ptrs (void* p0, void* p1) {
        pf = p0;
        pb = p1;
    }
#endif

    void destroy ()
    {
        if (defined) {
            destroy_vendor_plan(plan);
            defined = false;
        }
#if !defined(AMREX_USE_GPU)
        if (defined2) {
            destroy_vendor_plan(plan2);
            defined2 = false;
        }
#endif
    }

    template <Direction D>
    void init_r2c (Box const& box, T* pr, VendorComplex* pc, bool is_2d_transform = false, int ncomp = 1)
    {
        static_assert(D == Direction::forward || D == Direction::backward);
 
        int rank = is_2d_transform ? 2 : 1;
 
        kind = (D == Direction::forward) ? Kind::r2c_f : Kind::r2c_b;
        defined = true;
        pf = (void*)pr;
        pb = (void*)pc;
 
        int len[2] = {};
        if (rank == 1) {
            len[0] = box.length(0);
            len[1] = box.length(0); // Not used except for HIP. Yes it's `(0)`.
        } else {
            len[0] = box.length(1); // Most FFT libraries assume row-major ordering
            len[1] = box.length(0); // except for rocfft
        }
        int nr = (rank == 1) ? len[0] : len[0]*len[1];
        n = nr;
        int nc = (rank == 1) ? (len[0]/2+1) : (len[1]/2+1)*len[0];
#if (AMREX_SPACEDIM == 1)
        howmany = 1;
#else
        howmany = (rank == 1) ? AMREX_D_TERM(1, *box.length(1), *box.length(2))
                              : AMREX_D_TERM(1, *1            , *box.length(2));
#endif
        howmany *= ncomp;
 
        amrex::ignore_unused(nc);
 
#if defined(AMREX_USE_CUDA)
 
        AMREX_CUFFT_SAFE_CALL(cufftCreate(&plan));
        AMREX_CUFFT_SAFE_CALL(cufftSetAutoAllocation(plan, 0));
        std::size_t work_size;
        if constexpr (D == Direction::forward) {
            cufftType fwd_type = std::is_same_v<float,T> ? CUFFT_R2C : CUFFT_D2Z;
            AMREX_CUFFT_SAFE_CALL
                (cufftMakePlanMany(plan, rank, len, nullptr, 1, nr, nullptr, 1, nc, fwd_type, howmany, &work_size));
        } else {
            cufftType bwd_type = std::is_same_v<float,T> ? CUFFT_C2R : CUFFT_Z2D;
            AMREX_CUFFT_SAFE_CALL
                (cufftMakePlanMany(plan, rank, len, nullptr, 1, nc, nullptr, 1, nr, bwd_type, howmany, &work_size));
        }
 
#elif defined(AMREX_USE_HIP)
 
        auto prec = std::is_same_v<float,T> ? rocfft_precision_single : rocfft_precision_double;
        // switch to column-major ordering
        std::size_t length[2] = {std::size_t(len[1]), std::size_t(len[0])};
        if constexpr (D == Direction::forward) {
            AMREX_ROCFFT_SAFE_CALL
                (rocfft_plan_create(&plan, rocfft_placement_notinplace,
                                    rocfft_transform_type_real_forward, prec, rank,
                                    length, howmany, nullptr));
        } else {
            AMREX_ROCFFT_SAFE_CALL
                (rocfft_plan_create(&plan, rocfft_placement_notinplace,
                                    rocfft_transform_type_real_inverse, prec, rank,
                                    length, howmany, nullptr));
        }
 
#elif defined(AMREX_USE_SYCL)
 
        mkl_desc_r* pp;
        if (rank == 1) {
            pp = new mkl_desc_r(len[0]);
        } else {
            pp = new mkl_desc_r({std::int64_t(len[0]), std::int64_t(len[1])});
        }
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT,
                      oneapi::mkl::dft::config_value::NOT_INPLACE);
#else
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT, DFTI_NOT_INPLACE);
#endif
        pp->set_value(oneapi::mkl::dft::config_param::NUMBER_OF_TRANSFORMS, howmany);
        pp->set_value(oneapi::mkl::dft::config_param::FWD_DISTANCE, nr);
        pp->set_value(oneapi::mkl::dft::config_param::BWD_DISTANCE, nc);
        std::vector<std::int64_t> strides;
        strides.push_back(0);
        if (rank == 2) { strides.push_back(len[1]); }
        strides.push_back(1);
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides);
        // Do not set BWD_STRIDES
#else
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides.data());
        // Do not set BWD_STRIDES
#endif
        pp->set_value(oneapi::mkl::dft::config_param::WORKSPACE,
                      oneapi::mkl::dft::config_value::WORKSPACE_EXTERNAL);
        detail::assert_no_external_stream();
        pp->commit(amrex::Gpu::Device::streamQueue());
        plan = pp;
 
#else /* FFTW */
 
        if constexpr (std::is_same_v<float,T>) {
            if constexpr (D == Direction::forward) {
                plan = fftwf_plan_many_dft_r2c
                    (rank, len, howmany, pr, nullptr, 1, nr, pc, nullptr, 1, nc,
                     FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
            } else {
                plan = fftwf_plan_many_dft_c2r
                    (rank, len, howmany, pc, nullptr, 1, nc, pr, nullptr, 1, nr,
                     FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
            }
        } else {
            if constexpr (D == Direction::forward) {
                plan = fftw_plan_many_dft_r2c
                    (rank, len, howmany, pr, nullptr, 1, nr, pc, nullptr, 1, nc,
                     FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
            } else {
                plan = fftw_plan_many_dft_c2r
                    (rank, len, howmany, pc, nullptr, 1, nc, pr, nullptr, 1, nr,
                     FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
            }
        }
#endif
    }

    template <Direction D, int M>
    void init_r2c (IntVectND<M> const& fft_size, void* pbf, void* pbb, bool cache, int ncomp = 1);

    template <Direction D>
    void init_c2c (Box const& box, VendorComplex* p, int ncomp = 1, int ndims = 1)
    {
        static_assert(D == Direction::forward || D == Direction::backward);
 
        kind = (D == Direction::forward) ? Kind::c2c_f : Kind::c2c_b;
        defined = true;
        pf = (void*)p;
        pb = (void*)p;
 
        int len[3] = {};
 
        if (ndims == 1) {
            n = box.length(0);
            howmany = AMREX_D_TERM(1, *box.length(1), *box.length(2));
            howmany *= ncomp;
            len[0] = box.length(0);
        }
#if (AMREX_SPACEDIM >= 2)
        else if (ndims == 2) {
            n = box.length(0) * box.length(1);
#if (AMREX_SPACEDIM == 2)
            howmany = ncomp;
#else
            howmany = box.length(2) * ncomp;
#endif
            len[0] = box.length(1);
            len[1] = box.length(0);
        }
#if (AMREX_SPACEDIM == 3)
        else if (ndims == 3) {
            n = box.length(0) * box.length(1) * box.length(2);
            howmany = ncomp;
            len[0] = box.length(2);
            len[1] = box.length(1);
            len[2] = box.length(0);
        }
#endif
#endif
 
#if defined(AMREX_USE_CUDA)
        AMREX_CUFFT_SAFE_CALL(cufftCreate(&plan));
        AMREX_CUFFT_SAFE_CALL(cufftSetAutoAllocation(plan, 0));
 
        cufftType t = std::is_same_v<float,T> ? CUFFT_C2C : CUFFT_Z2Z;
        std::size_t work_size;
        AMREX_CUFFT_SAFE_CALL
            (cufftMakePlanMany(plan, ndims, len, nullptr, 1, n, nullptr, 1, n, t, howmany, &work_size));
 
#elif defined(AMREX_USE_HIP)
 
        auto prec = std::is_same_v<float,T> ? rocfft_precision_single
                                            : rocfft_precision_double;
        auto dir= (D == Direction::forward) ? rocfft_transform_type_complex_forward
                                            : rocfft_transform_type_complex_inverse;
        std::size_t length[3];
        if (ndims == 1) {
            length[0] = len[0];
        } else if (ndims == 2) {
            length[0] = len[1];
            length[1] = len[0];
        } else {
            length[0] = len[2];
            length[1] = len[1];
            length[2] = len[0];
        }
        AMREX_ROCFFT_SAFE_CALL
            (rocfft_plan_create(&plan, rocfft_placement_inplace, dir, prec, ndims,
                                length, howmany, nullptr));
 
#elif defined(AMREX_USE_SYCL)
 
        mkl_desc_c* pp;
        if (ndims == 1) {
            pp = new mkl_desc_c(n);
        } else if (ndims == 2) {
            pp = new mkl_desc_c({std::int64_t(len[0]), std::int64_t(len[1])});
        } else {
            pp = new mkl_desc_c({std::int64_t(len[0]), std::int64_t(len[1]), std::int64_t(len[2])});
        }
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT,
                      oneapi::mkl::dft::config_value::INPLACE);
#else
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT, DFTI_INPLACE);
#endif
        pp->set_value(oneapi::mkl::dft::config_param::NUMBER_OF_TRANSFORMS, howmany);
        pp->set_value(oneapi::mkl::dft::config_param::FWD_DISTANCE, n);
        pp->set_value(oneapi::mkl::dft::config_param::BWD_DISTANCE, n);
        std::vector<std::int64_t> strides(ndims+1);
        strides[0] = 0;
        strides[ndims] = 1;
        for (int i = ndims-1; i >= 1; --i) {
            strides[i] = strides[i+1] * len[i];
        }
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides);
        pp->set_value(oneapi::mkl::dft::config_param::BWD_STRIDES, strides);
#else
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides.data());
        pp->set_value(oneapi::mkl::dft::config_param::BWD_STRIDES, strides.data());
#endif
        pp->set_value(oneapi::mkl::dft::config_param::WORKSPACE,
                      oneapi::mkl::dft::config_value::WORKSPACE_EXTERNAL);
        detail::assert_no_external_stream();
        pp->commit(amrex::Gpu::Device::streamQueue());
        plan = pp;
 
#else /* FFTW */
 
        if constexpr (std::is_same_v<float,T>) {
            if constexpr (D == Direction::forward) {
                plan = fftwf_plan_many_dft
                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
                     FFTW_ESTIMATE);
            } else {
                plan = fftwf_plan_many_dft
                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, +1,
                     FFTW_ESTIMATE);
            }
        } else {
            if constexpr (D == Direction::forward) {
                plan = fftw_plan_many_dft
                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
                     FFTW_ESTIMATE);
            } else {
                plan = fftw_plan_many_dft
                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, +1,
                     FFTW_ESTIMATE);
            }
        }
#endif
    }
 
#ifndef AMREX_USE_GPU
    template <Direction D>
    fftw_r2r_kind get_fftw_kind (std::pair<Boundary,Boundary> const& bc)
    {
        if (bc.first == Boundary::even && bc.second == Boundary::even)
        {
            return (D == Direction::forward) ? FFTW_REDFT10 : FFTW_REDFT01;
        }
        else if (bc.first == Boundary::even && bc.second == Boundary::odd)
        {
            return FFTW_REDFT11;
        }
        else if (bc.first == Boundary::odd && bc.second == Boundary::even)
        {
            return FFTW_RODFT11;
        }
        else if (bc.first == Boundary::odd && bc.second == Boundary::odd)
        {
            return (D == Direction::forward) ? FFTW_RODFT10 : FFTW_RODFT01;
        }
        else {
            amrex::Abort("FFT: unsupported BC");
            return fftw_r2r_kind{};
        }
 
    }
#endif

    template <Direction D>
    Kind get_r2r_kind (std::pair<Boundary,Boundary> const& bc)
    {
        if (bc.first == Boundary::even && bc.second == Boundary::even)
        {
            return (D == Direction::forward) ? Kind::r2r_ee_f : Kind::r2r_ee_b;
        }
        else if (bc.first == Boundary::even && bc.second == Boundary::odd)
        {
            return Kind::r2r_eo;
        }
        else if (bc.first == Boundary::odd && bc.second == Boundary::even)
        {
            return Kind::r2r_oe;
        }
        else if (bc.first == Boundary::odd && bc.second == Boundary::odd)
        {
            return (D == Direction::forward) ? Kind::r2r_oo_f : Kind::r2r_oo_b;
        }
        else {
            amrex::Abort("FFT: unsupported BC");
            return Kind::none;
        }
 
    }

    template <Direction D>
    void init_r2r (Box const& box, T* p, std::pair<Boundary,Boundary> const& bc,
                   int howmany_initval = 1)
    {
        static_assert(D == Direction::forward || D == Direction::backward);
 
        kind = get_r2r_kind<D>(bc);
        defined = true;
        pf = (void*)p;
        pb = (void*)p;
 
        n = box.length(0);
        howmany = AMREX_D_TERM(howmany_initval, *box.length(1), *box.length(2));
 
#if defined(AMREX_USE_GPU)
        int nex=0;
        if (bc.first == Boundary::odd && bc.second == Boundary::odd &&
            Direction::forward == D) {
            nex = 2*n;
        } else if (bc.first == Boundary::odd && bc.second == Boundary::odd &&
            Direction::backward == D) {
            nex = 4*n;
        } else if (bc.first == Boundary::even && bc.second == Boundary::even &&
            Direction::forward == D) {
            nex = 2*n;
        } else if (bc.first == Boundary::even && bc.second == Boundary::even &&
            Direction::backward == D) {
            nex = 4*n;
        } else if ((bc.first == Boundary::even && bc.second == Boundary::odd) ||
                   (bc.first == Boundary::odd && bc.second == Boundary::even)) {
            nex = 4*n;
        } else {
            amrex::Abort("FFT: unsupported BC");
        }
        int nc = (nex/2) + 1;
 
#if defined (AMREX_USE_CUDA)
 
        AMREX_CUFFT_SAFE_CALL(cufftCreate(&plan));
        AMREX_CUFFT_SAFE_CALL(cufftSetAutoAllocation(plan, 0));
        cufftType fwd_type = std::is_same_v<float,T> ? CUFFT_R2C : CUFFT_D2Z;
        std::size_t work_size;
        AMREX_CUFFT_SAFE_CALL
            (cufftMakePlanMany(plan, 1, &nex, nullptr, 1, nc*2, nullptr, 1, nc, fwd_type, howmany, &work_size));
 
#elif defined(AMREX_USE_HIP)
 
        amrex::ignore_unused(nc);
        auto prec = std::is_same_v<float,T> ? rocfft_precision_single : rocfft_precision_double;
        const std::size_t length = nex;
        AMREX_ROCFFT_SAFE_CALL
            (rocfft_plan_create(&plan, rocfft_placement_inplace,
                                rocfft_transform_type_real_forward, prec, 1,
                                &length, howmany, nullptr));
 
#elif defined(AMREX_USE_SYCL)
 
        auto* pp = new mkl_desc_r(nex);
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT,
                      oneapi::mkl::dft::config_value::INPLACE);
#else
        pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT, DFTI_INPLACE);
#endif
        pp->set_value(oneapi::mkl::dft::config_param::NUMBER_OF_TRANSFORMS, howmany);
        pp->set_value(oneapi::mkl::dft::config_param::FWD_DISTANCE, nc*2);
        pp->set_value(oneapi::mkl::dft::config_param::BWD_DISTANCE, nc);
        std::vector<std::int64_t> strides = {0,1};
#ifndef AMREX_USE_MKL_DFTI_2024
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides);
        pp->set_value(oneapi::mkl::dft::config_param::BWD_STRIDES, strides);
#else
        pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides.data());
        pp->set_value(oneapi::mkl::dft::config_param::BWD_STRIDES, strides.data());
#endif
        pp->set_value(oneapi::mkl::dft::config_param::WORKSPACE,
                      oneapi::mkl::dft::config_value::WORKSPACE_EXTERNAL);
        detail::assert_no_external_stream();
        pp->commit(amrex::Gpu::Device::streamQueue());
        plan = pp;
 
#endif
 
#else /* FFTW */
        auto fftw_kind = get_fftw_kind<D>(bc);
        if constexpr (std::is_same_v<float,T>) {
            plan = fftwf_plan_many_r2r
                (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, &fftw_kind,
                 FFTW_ESTIMATE);
        } else {
            plan = fftw_plan_many_r2r
                (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, &fftw_kind,
                 FFTW_ESTIMATE);
        }
#endif
    }

    template <Direction D>
    void init_r2r (Box const& box, VendorComplex* pc,
                   std::pair<Boundary,Boundary> const& bc)
    {
        static_assert(D == Direction::forward || D == Direction::backward);
 
        auto* p = (T*)pc;
 
#if defined(AMREX_USE_GPU)
 
        init_r2r<D>(box, p, bc, 2);
        r2r_data_is_complex = true;
 
#else
 
        kind = get_r2r_kind<D>(bc);
        defined = true;
        pf = (void*)p;
        pb = (void*)p;
 
        n = box.length(0);
        howmany = AMREX_D_TERM(1, *box.length(1), *box.length(2));
 
        defined2 = true;
        auto fftw_kind = get_fftw_kind<D>(bc);
        if constexpr (std::is_same_v<float,T>) {
            plan = fftwf_plan_many_r2r
                (1, &n, howmany, p, nullptr, 2, n*2, p, nullptr, 2, n*2, &fftw_kind,
                 FFTW_ESTIMATE);
            plan2 = fftwf_plan_many_r2r
                (1, &n, howmany, p+1, nullptr, 2, n*2, p+1, nullptr, 2, n*2, &fftw_kind,
                 FFTW_ESTIMATE);
        } else {
            plan = fftw_plan_many_r2r
                (1, &n, howmany, p, nullptr, 2, n*2, p, nullptr, 2, n*2, &fftw_kind,
                 FFTW_ESTIMATE);
            plan2 = fftw_plan_many_r2r
                (1, &n, howmany, p+1, nullptr, 2, n*2, p+1, nullptr, 2, n*2, &fftw_kind,
                 FFTW_ESTIMATE);
        }
#endif
    }

    template <Direction D>
    void compute_r2c ()
    {
        static_assert(D == Direction::forward || D == Direction::backward);
        if (!defined) { return; }
 
        using TI = std::conditional_t<(D == Direction::forward), T, VendorComplex>;
        using TO = std::conditional_t<(D == Direction::backward), T, VendorComplex>;
        auto* pi = (TI*)((D == Direction::forward) ? pf : pb);
        auto* po = (TO*)((D == Direction::forward) ? pb : pf);
 
#if defined(AMREX_USE_CUDA)
        AMREX_CUFFT_SAFE_CALL(cufftSetStream(plan, Gpu::gpuStream()));
 
        std::size_t work_size = 0;
        AMREX_CUFFT_SAFE_CALL(cufftGetSize(plan, &work_size));
 
        auto* work_area = The_Arena()->alloc(work_size);
        AMREX_CUFFT_SAFE_CALL(cufftSetWorkArea(plan, work_area));
 
        if constexpr (D == Direction::forward) {
            if constexpr (std::is_same_v<float,T>) {
                AMREX_CUFFT_SAFE_CALL(cufftExecR2C(plan, pi, po));
            } else {
                AMREX_CUFFT_SAFE_CALL(cufftExecD2Z(plan, pi, po));
            }
        } else {
            if constexpr (std::is_same_v<float,T>) {
                AMREX_CUFFT_SAFE_CALL(cufftExecC2R(plan, pi, po));
            } else {
                AMREX_CUFFT_SAFE_CALL(cufftExecZ2D(plan, pi, po));
            }
        }
        Gpu::streamSynchronize();
        The_Arena()->free(work_area);
#elif defined(AMREX_USE_HIP)
        detail::hip_execute(plan, (void**)&pi, (void**)&po);
#elif defined(AMREX_USE_SYCL)
        detail::sycl_execute<T,D>(std::get<0>(plan), pi, po);
#else
        amrex::ignore_unused(pi,po);
        if constexpr (std::is_same_v<float,T>) {
            fftwf_execute(plan);
        } else {
            fftw_execute(plan);
        }
#endif
    }

    template <Direction D>
    void compute_c2c ()
    {
        static_assert(D == Direction::forward || D == Direction::backward);
        if (!defined) { return; }
 
        auto* p = (VendorComplex*)pf;
 
#if defined(AMREX_USE_CUDA)
        AMREX_CUFFT_SAFE_CALL(cufftSetStream(plan, Gpu::gpuStream()));
 
        std::size_t work_size = 0;
        AMREX_CUFFT_SAFE_CALL(cufftGetSize(plan, &work_size));
 
        auto* work_area = The_Arena()->alloc(work_size);
        AMREX_CUFFT_SAFE_CALL(cufftSetWorkArea(plan, work_area));
 
        auto dir = (D == Direction::forward) ? CUFFT_FORWARD : CUFFT_INVERSE;
        if constexpr (std::is_same_v<float,T>) {
            AMREX_CUFFT_SAFE_CALL(cufftExecC2C(plan, p, p, dir));
        } else {
            AMREX_CUFFT_SAFE_CALL(cufftExecZ2Z(plan, p, p, dir));
        }
        Gpu::streamSynchronize();
        The_Arena()->free(work_area);
#elif defined(AMREX_USE_HIP)
        detail::hip_execute(plan, (void**)&p, (void**)&p);
#elif defined(AMREX_USE_SYCL)
        detail::sycl_execute<T,D>(std::get<1>(plan), p, p);
#else
        amrex::ignore_unused(p);
        if constexpr (std::is_same_v<float,T>) {
            fftwf_execute(plan);
        } else {
            fftw_execute(plan);
        }
#endif
    }
 
#ifdef AMREX_USE_GPU
    [[nodiscard]] void* alloc_scratch_space () const
    {
        int nc = 0;
        if (kind == Kind::r2r_oo_f || kind == Kind::r2r_ee_f) {
            nc = n + 1;
        } else if (kind == Kind::r2r_oo_b || kind == Kind::r2r_ee_b ||
                   kind == Kind::r2r_oe || kind == Kind::r2r_eo) {
            nc = 2*n+1;
        } else {
            amrex::Abort("FFT: alloc_scratch_space: unsupported kind");
        }
        return The_Arena()->alloc(sizeof(GpuComplex<T>)*nc*howmany);
    }

    static void free_scratch_space (void* p) { The_Arena()->free(p); }

    void pack_r2r_buffer (void* pbuf, T const* psrc) const
    {
        auto* pdst = (T*) pbuf;
        if (kind == Kind::r2r_oo_f || kind == Kind::r2r_ee_f) {
            T sign = (kind == Kind::r2r_oo_f) ? T(-1) : T(1);
            int ostride = (n+1)*2;
            int istride = n;
            int nex = 2*n;
            int norig = n;
            Long nelems = Long(nex)*howmany;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto* po = pdst + (2*batch+ir)*ostride + i;
                        auto const* pi = psrc + 2*batch*istride + ir;
                        if (i < norig) {
                            *po = pi[i*2];
                        } else {
                            *po = sign * pi[(2*norig-1-i)*2];
                        }
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    auto* po = pdst + batch*ostride + i;
                    auto const* pi = psrc + batch*istride;
                    if (i < norig) {
                        *po = pi[i];
                    } else {
                        *po = sign * pi[2*norig-1-i];
                    }
                });
            }
        } else if (kind == Kind::r2r_oo_b) {
            int ostride = (2*n+1)*2;
            int istride = n;
            int nex = 4*n;
            int norig = n;
            Long nelems = Long(nex)*howmany;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto* po = pdst + (2*batch+ir)*ostride + i;
                        auto const* pi = psrc + 2*batch*istride + ir;
                        if (i < norig) {
                            *po = pi[i*2];
                        } else if (i < (2*norig-1)) {
                            *po = pi[(2*norig-2-i)*2];
                        } else if (i == (2*norig-1)) {
                            *po = T(0);
                        } else if (i < (3*norig)) {
                            *po = -pi[(i-2*norig)*2];
                        } else if (i < (4*norig-1)) {
                            *po = -pi[(4*norig-2-i)*2];
                        } else {
                            *po = T(0);
                        }
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    auto* po = pdst + batch*ostride + i;
                    auto const* pi = psrc + batch*istride;
                    if (i < norig) {
                        *po = pi[i];
                    } else if (i < (2*norig-1)) {
                        *po = pi[2*norig-2-i];
                    } else if (i == (2*norig-1)) {
                        *po = T(0);
                    } else if (i < (3*norig)) {
                        *po = -pi[i-2*norig];
                    } else if (i < (4*norig-1)) {
                        *po = -pi[4*norig-2-i];
                    } else {
                        *po = T(0);
                    }
                });
            }
        } else if (kind == Kind::r2r_ee_b) {
            int ostride = (2*n+1)*2;
            int istride = n;
            int nex = 4*n;
            int norig = n;
            Long nelems = Long(nex)*howmany;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto* po = pdst + (2*batch+ir)*ostride + i;
                        auto const* pi = psrc + 2*batch*istride + ir;
                        if (i < norig) {
                            *po = pi[i*2];
                        } else if (i == norig) {
                            *po = T(0);
                        } else if (i < (2*norig+1)) {
                            *po = -pi[(2*norig-i)*2];
                        } else if (i < (3*norig)) {
                            *po = -pi[(i-2*norig)*2];
                        } else if (i == 3*norig) {
                            *po = T(0);
                        } else {
                            *po = pi[(4*norig-i)*2];
                        }
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    auto* po = pdst + batch*ostride + i;
                    auto const* pi = psrc + batch*istride;
                    if (i < norig) {
                        *po = pi[i];
                    } else if (i == norig) {
                        *po = T(0);
                    } else if (i < (2*norig+1)) {
                        *po = -pi[2*norig-i];
                    } else if (i < (3*norig)) {
                        *po = -pi[i-2*norig];
                    } else if (i == 3*norig) {
                        *po = T(0);
                    } else {
                        *po = pi[4*norig-i];
                    }
                });
            }
        } else if (kind == Kind::r2r_eo) {
            int ostride = (2*n+1)*2;
            int istride = n;
            int nex = 4*n;
            int norig = n;
            Long nelems = Long(nex)*howmany;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto* po = pdst + (2*batch+ir)*ostride + i;
                        auto const* pi = psrc + 2*batch*istride + ir;
                        if (i < norig) {
                            *po = pi[i*2];
                        } else if (i < (2*norig)) {
                            *po = -pi[(2*norig-1-i)*2];
                        } else if (i < (3*norig)) {
                            *po = -pi[(i-2*norig)*2];
                        } else {
                            *po = pi[(4*norig-1-i)*2];
                        }
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    auto* po = pdst + batch*ostride + i;
                    auto const* pi = psrc + batch*istride;
                    if (i < norig) {
                        *po = pi[i];
                    } else if (i < (2*norig)) {
                        *po = -pi[2*norig-1-i];
                    } else if (i < (3*norig)) {
                        *po = -pi[i-2*norig];
                    } else {
                        *po = pi[4*norig-1-i];
                    }
                });
            }
        } else if (kind == Kind::r2r_oe) {
            int ostride = (2*n+1)*2;
            int istride = n;
            int nex = 4*n;
            int norig = n;
            Long nelems = Long(nex)*howmany;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto* po = pdst + (2*batch+ir)*ostride + i;
                        auto const* pi = psrc + 2*batch*istride + ir;
                        if (i < norig) {
                            *po = pi[i*2];
                        } else if (i < (2*norig)) {
                            *po = pi[(2*norig-1-i)*2];
                        } else if (i < (3*norig)) {
                            *po = -pi[(i-2*norig)*2];
                        } else {
                            *po = -pi[(4*norig-1-i)*2];
                        }
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(nex);
                    auto i = int(ielem - batch*nex);
                    auto* po = pdst + batch*ostride + i;
                    auto const* pi = psrc + batch*istride;
                    if (i < norig) {
                        *po = pi[i];
                    } else if (i < (2*norig)) {
                        *po = pi[2*norig-1-i];
                    } else if (i < (3*norig)) {
                        *po = -pi[i-2*norig];
                    } else {
                        *po = -pi[4*norig-1-i];
                    }
                });
            }
        } else {
            amrex::Abort("FFT: pack_r2r_buffer: unsupported kind");
        }
    }

    void unpack_r2r_buffer (T* pdst, void const* pbuf) const
    {
        auto const* psrc = (GpuComplex<T> const*) pbuf;
        int norig = n;
        Long nelems = Long(norig)*howmany;
        int ostride = n;
 
        if (kind == Kind::r2r_oo_f) {
            int istride = n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(k+1)/T(2*norig));
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+k+1];
                        pdst[2*batch*ostride+ir+k*2] = s * yk.real() - c * yk.imag();
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(k+1)/T(2*norig));
                    auto const& yk = psrc[batch*istride+k+1];
                    pdst[batch*ostride+k] = s * yk.real() - c * yk.imag();
                });
            }
        } else if (kind == Kind::r2r_oo_b) {
            int istride = 2*n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(2*k+1)/T(2*norig));
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+2*k+1];
                        pdst[2*batch*ostride+ir+k*2] = T(0.5)*(s * yk.real() - c * yk.imag());
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(2*k+1)/T(2*norig));
                    auto const& yk = psrc[batch*istride+2*k+1];
                    pdst[batch*ostride+k] = T(0.5)*(s * yk.real() - c * yk.imag());
                });
            }
        } else if (kind == Kind::r2r_ee_f) {
            int istride = n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(k)/T(2*norig));
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+k];
                        pdst[2*batch*ostride+ir+k*2] = c * yk.real() + s * yk.imag();
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi(T(k)/T(2*norig));
                    auto const& yk = psrc[batch*istride+k];
                    pdst[batch*ostride+k] = c * yk.real() + s * yk.imag();
                });
            }
        } else if (kind == Kind::r2r_ee_b) {
            int istride = 2*n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+2*k+1];
                        pdst[2*batch*ostride+ir+k*2] = T(0.5) * yk.real();
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto const& yk = psrc[batch*istride+2*k+1];
                    pdst[batch*ostride+k] = T(0.5) * yk.real();
                });
            }
        } else if (kind == Kind::r2r_eo) {
            int istride = 2*n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi((k+T(0.5))/T(2*norig));
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+2*k+1];
                        pdst[2*batch*ostride+ir+k*2] = T(0.5) * (c * yk.real() + s * yk.imag());
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi((k+T(0.5))/T(2*norig));
                    auto const& yk = psrc[batch*istride+2*k+1];
                    pdst[batch*ostride+k] = T(0.5) * (c * yk.real() + s * yk.imag());
                });
            }
        } else if (kind == Kind::r2r_oe) {
            int istride = 2*n+1;
            if (r2r_data_is_complex) {
                ParallelForOMP(nelems/2, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi((k+T(0.5))/T(2*norig));
                    for (int ir = 0; ir < 2; ++ir) {
                        auto const& yk = psrc[(2*batch+ir)*istride+2*k+1];
                        pdst[2*batch*ostride+ir+k*2] = T(0.5) * (s * yk.real() - c * yk.imag());
                    }
                });
            } else {
                ParallelForOMP(nelems, [=] AMREX_GPU_DEVICE (Long ielem)
                {
                    auto batch = ielem / Long(norig);
                    auto k = int(ielem - batch*norig);
                    auto [s, c] = Math::sincospi((k+T(0.5))/T(2*norig));
                    auto const& yk = psrc[batch*istride+2*k+1];
                    pdst[batch*ostride+k] = T(0.5) * (s * yk.real() - c * yk.imag());
                });
            }
        } else {
            amrex::Abort("FFT: unpack_r2r_buffer: unsupported kind");
        }
    }
#endif

    template <Direction D>
    void compute_r2r ()
    {
        static_assert(D == Direction::forward || D == Direction::backward);
        if (!defined) { return; }
 
#if defined(AMREX_USE_GPU)
 
        auto* pscratch = alloc_scratch_space();
 
        pack_r2r_buffer(pscratch, (T*)((D == Direction::forward) ? pf : pb));
 
#if defined(AMREX_USE_CUDA)
 
        AMREX_CUFFT_SAFE_CALL(cufftSetStream(plan, Gpu::gpuStream()));
 
        std::size_t work_size = 0;
        AMREX_CUFFT_SAFE_CALL(cufftGetSize(plan, &work_size));
 
        auto* work_area = The_Arena()->alloc(work_size);
        AMREX_CUFFT_SAFE_CALL(cufftSetWorkArea(plan, work_area));
 
        if constexpr (std::is_same_v<float,T>) {
            AMREX_CUFFT_SAFE_CALL(cufftExecR2C(plan, (T*)pscratch, (VendorComplex*)pscratch));
        } else {
            AMREX_CUFFT_SAFE_CALL(cufftExecD2Z(plan, (T*)pscratch, (VendorComplex*)pscratch));
        }
 
#elif defined(AMREX_USE_HIP)
        detail::hip_execute(plan, (void**)&pscratch, (void**)&pscratch);
#elif defined(AMREX_USE_SYCL)
        detail::sycl_execute<T,Direction::forward>(std::get<0>(plan), (T*)pscratch, (VendorComplex*)pscratch);
#endif
 
        unpack_r2r_buffer((T*)((D == Direction::forward) ? pb : pf), pscratch);
 
        Gpu::streamSynchronize();
        free_scratch_space(pscratch);
#if defined(AMREX_USE_CUDA)
        The_Arena()->free(work_area);
#endif
 
#else /* FFTW */
 
        if constexpr (std::is_same_v<float,T>) {
            fftwf_execute(plan);
            if (defined2) { fftwf_execute(plan2); }
        } else {
            fftw_execute(plan);
            if (defined2) { fftw_execute(plan2); }
        }
 
#endif
    }

    static void destroy_vendor_plan (VendorPlan plan)
    {
#if defined(AMREX_USE_CUDA)
        AMREX_CUFFT_SAFE_CALL(cufftDestroy(plan));
#elif defined(AMREX_USE_HIP)
        AMREX_ROCFFT_SAFE_CALL(rocfft_plan_destroy(plan));
#elif defined(AMREX_USE_SYCL)
        std::visit([](auto&& p) { delete p; }, plan);
#else
        if constexpr (std::is_same_v<float,T>) {
            fftwf_destroy_plan(plan);
        } else {
            fftw_destroy_plan(plan);
        }
#endif
    }
};

 
using Key = std::tuple<IntVectND<3>,int,Direction,Kind>;
using PlanD = typename Plan<double>::VendorPlan;
using PlanF = typename Plan<float>::VendorPlan;
 
namespace detail {
    PlanD* get_vendor_plan_d (Key const& key);
    PlanF* get_vendor_plan_f (Key const& key);
 
    void add_vendor_plan_d (Key const& key, PlanD plan);
    void add_vendor_plan_f (Key const& key, PlanF plan);
}
 
template <typename T>
template <Direction D, int M>
void Plan<T>::init_r2c (IntVectND<M> const& fft_size, void* pbf, void* pbb, bool cache, int ncomp)
{
    static_assert(D == Direction::forward || D == Direction::backward);
 
    kind = (D == Direction::forward) ? Kind::r2c_f : Kind::r2c_b;
    defined = true;
    pf = pbf;
    pb = pbb;
 
    n = 1;
    for (auto s : fft_size) { n *= s; }
    howmany = ncomp;
 
#if defined(AMREX_USE_GPU)
    Key key = {fft_size.template expand<3>(), ncomp, D, kind};
    if (cache) {
        VendorPlan* cached_plan = nullptr;
        if constexpr (std::is_same_v<float,T>) {
            cached_plan = detail::get_vendor_plan_f(key);
        } else {
            cached_plan = detail::get_vendor_plan_d(key);
        }
        if (cached_plan) {
            plan = *cached_plan;
            return;
        }
    }
#else
    amrex::ignore_unused(cache);
#endif
 
    int len[M];
    for (int i = 0; i < M; ++i) {
        len[i] = fft_size[M-1-i];
    }
 
    int nc = fft_size[0]/2+1;
    for (int i = 1; i < M; ++i) {
        nc *= fft_size[i];
    }
 
#if defined(AMREX_USE_CUDA)
 
    AMREX_CUFFT_SAFE_CALL(cufftCreate(&plan));
    AMREX_CUFFT_SAFE_CALL(cufftSetAutoAllocation(plan, 0));
    cufftType type;
    int n_in, n_out;
    if constexpr (D == Direction::forward) {
        type = std::is_same_v<float,T> ? CUFFT_R2C : CUFFT_D2Z;
        n_in = n;
        n_out = nc;
    } else {
        type = std::is_same_v<float,T> ? CUFFT_C2R : CUFFT_Z2D;
        n_in = nc;
        n_out = n;
    }
    std::size_t work_size;
    AMREX_CUFFT_SAFE_CALL
        (cufftMakePlanMany(plan, M, len, nullptr, 1, n_in, nullptr, 1, n_out, type, howmany, &work_size));
 
#elif defined(AMREX_USE_HIP)
 
    auto prec = std::is_same_v<float,T> ? rocfft_precision_single : rocfft_precision_double;
    std::size_t length[M];
    for (int idim = 0; idim < M; ++idim) { length[idim] = fft_size[idim]; }
    if constexpr (D == Direction::forward) {
        AMREX_ROCFFT_SAFE_CALL
            (rocfft_plan_create(&plan, rocfft_placement_notinplace,
                                rocfft_transform_type_real_forward, prec, M,
                                length, howmany, nullptr));
    } else {
        AMREX_ROCFFT_SAFE_CALL
            (rocfft_plan_create(&plan, rocfft_placement_notinplace,
                                rocfft_transform_type_real_inverse, prec, M,
                                length, howmany, nullptr));
    }
 
#elif defined(AMREX_USE_SYCL)
 
    mkl_desc_r* pp;
    if (M == 1) {
        pp = new mkl_desc_r(fft_size[0]);
    } else {
        std::vector<std::int64_t> len64(M);
        for (int idim = 0; idim < M; ++idim) {
            len64[idim] = len[idim];
        }
        pp = new mkl_desc_r(len64);
    }
#ifndef AMREX_USE_MKL_DFTI_2024
    pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT,
                  oneapi::mkl::dft::config_value::NOT_INPLACE);
#else
    pp->set_value(oneapi::mkl::dft::config_param::PLACEMENT, DFTI_NOT_INPLACE);
#endif
    pp->set_value(oneapi::mkl::dft::config_param::NUMBER_OF_TRANSFORMS, howmany);
    pp->set_value(oneapi::mkl::dft::config_param::FWD_DISTANCE, n);
    pp->set_value(oneapi::mkl::dft::config_param::BWD_DISTANCE, nc);
    std::vector<std::int64_t> strides(M+1);
    strides[0] = 0;
    strides[M] = 1;
    for (int i = M-1; i >= 1; --i) {
        strides[i] = strides[i+1] * fft_size[M-1-i];
    }
 
#ifndef AMREX_USE_MKL_DFTI_2024
    pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides);
    // Do not set BWD_STRIDES
#else
    pp->set_value(oneapi::mkl::dft::config_param::FWD_STRIDES, strides.data());
    // Do not set BWD_STRIDES
#endif
    pp->set_value(oneapi::mkl::dft::config_param::WORKSPACE,
                  oneapi::mkl::dft::config_value::WORKSPACE_EXTERNAL);
    detail::assert_no_external_stream();
    pp->commit(amrex::Gpu::Device::streamQueue());
    plan = pp;
 
#else /* FFTW */
 
    if (pf == nullptr || pb == nullptr) {
        defined = false;
        return;
    }
 
    if constexpr (std::is_same_v<float,T>) {
        if constexpr (D == Direction::forward) {
            plan = fftwf_plan_many_dft_r2c
                (M, len, howmany, (float*)pf, nullptr, 1, n, (fftwf_complex*)pb, nullptr, 1, nc,
                 FFTW_ESTIMATE);
        } else {
            plan = fftwf_plan_many_dft_c2r
                (M, len, howmany, (fftwf_complex*)pb, nullptr, 1, nc, (float*)pf, nullptr, 1, n,
                 FFTW_ESTIMATE);
        }
    } else {
        if constexpr (D == Direction::forward) {
            plan = fftw_plan_many_dft_r2c
                (M, len, howmany, (double*)pf, nullptr, 1, n, (fftw_complex*)pb, nullptr, 1, nc,
                 FFTW_ESTIMATE);
        } else {
            plan = fftw_plan_many_dft_c2r
                (M, len, howmany, (fftw_complex*)pb, nullptr, 1, nc, (double*)pf, nullptr, 1, n,
                 FFTW_ESTIMATE);
        }
    }
#endif
 
#if defined(AMREX_USE_GPU)
    if (cache) {
        if constexpr (std::is_same_v<float,T>) {
            detail::add_vendor_plan_f(key, plan);
        } else {
            detail::add_vendor_plan_d(key, plan);
        }
    }
#endif
}

namespace detail
{
    DistributionMapping make_iota_distromap (Long n);

    CommStrategy resolve_comm_strategy (CommStrategy s, Long max_fab_bytes,
                                        std::size_t max_block_bytes, bool dense);
 
    struct TransposeMetaData;

    inline std::size_t max_block_bytes
        (std::initializer_list<std::unique_ptr<TransposeMetaData> const*> mds);

    inline bool transposes_dense
        (std::initializer_list<std::unique_ptr<TransposeMetaData> const*> mds);
 
    template <typename FA>
    Long max_fab_bytes (FA const& fa)
    {
        using T = typename FA::FABType::value_type::value_type;
        Long n = 0;
        if (!fa.empty()) {
            auto const& ba = fa.boxArray();
            for (int i = 0, nb = int(ba.size()); i < nb; ++i) {
                n = std::max(n, ba[i].numPts());
            }
        }
        return n * fa.nComp() * Long(sizeof(T));
    }

    struct TransposePattern
    {
        std::size_t block_bytes = 0;
        bool dense = true;
    };

    struct TransposeMetaData
    {
        template <typename... Args>
        explicit TransposeMetaData (TransposePattern pattern_, Args&&... args)
            : cmd(std::forward<Args>(args)...), pattern(pattern_) {}
 
        MultiBlockCommMetaData cmd;
        NonLocalBC::CommHandler handler;
        TransposePattern pattern;
    };

    int count_active_procs (std::initializer_list<FabArrayBase const*> fas);

    MPI_Comm make_active_comm (int n_active, bool& owned);

    void free_active_comm (MPI_Comm comm, bool owned);

    inline Vector<std::pair<int,int>> box_intervals (Vector<Box> const& boxes, int dir)
    {
        Vector<std::pair<int,int>> iv;
        iv.reserve(boxes.size());
        for (auto const& b : boxes) {
            iv.emplace_back(b.smallEnd(dir), b.bigEnd(dir));
        }
        std::sort(iv.begin(), iv.end());
        iv.erase(std::unique(iv.begin(), iv.end()), iv.end());
        return iv;
    }

    template <typename FA, typename DTOS>
    TransposePattern transpose_pattern (FA const& dst, Box const& dstbox, FA const& src,
                                        DTOS const& dtos, int n_active)
    {
        using T = typename FA::FABType::value_type::value_type;
        auto const& dstba = dst.boxArray();
        auto const& srcba = src.boxArray();
 
        // Destination boxes clipped and mapped into the source index space,
        // and their indices.
        Vector<Box> dboxes;
        Vector<int> dindex;
        for (int i = 0, n = int(dstba.size()); i < n; ++i) {
            Box const b = dstba[i] & dstbox;
            if (b.ok()) {
                dboxes.push_back(NonLocalBC::Image(dtos, b).setType(srcba.ixType()));
                dindex.push_back(i);
            }
        }
        Vector<Box> sboxes;
        for (int k = 0, n = int(srcba.size()); k < n; ++k) {
            sboxes.push_back(srcba[k]);
        }
 
        // Per direction, count the overlapping (dst,src) interval pairs and
        // find the longest overlap. Their products are the number of
        // overlapping box pairs and the largest overlap.
        Long npairs = 1;
        Long maxpts = 1;
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            auto const div = box_intervals(dboxes, idim);
            auto const siv = box_intervals(sboxes, idim);
            Long nover = 0;
            int maxlen = 0;
            Long j0 = 0;
            for (auto const& a : div) {
                while (j0 < siv.size() && siv[j0].second < a.first) { ++j0; }
                for (auto j = j0; j < siv.size() && siv[j].first <= a.second; ++j) {
                    ++nover;
                    maxlen = std::max(maxlen, std::min(a.second, siv[j].second)
                                            - std::max(a.first, siv[j].first) + 1);
                }
            }
            npairs *= nover;
            maxpts *= maxlen;
        }
 
        // Both layouts put box j on local rank j, so pairs (j,j) are local
        // copies, not messages.
        for (Long i = 0; i < dboxes.size(); ++i) {
            int const j = dindex[i];
            if (j < sboxes.size() && dboxes[i].intersects(sboxes[j])) { --npairs; }
        }
        if (npairs <= 0) { npairs = 0; maxpts = 0; }
 
        TransposePattern p;
        p.block_bytes = amrex::aligned_size(alignof(T), std::size_t(maxpts) * dst.nComp() * sizeof(T));
        // Dense if at least half of all pairs in the communicator exchange
        // data. Alltoall pads every pair to the largest message.
        p.dense = (2*npairs >= Long(n_active) * Long(n_active));
        return p;
    }
 
    inline std::size_t max_block_bytes
        (std::initializer_list<std::unique_ptr<TransposeMetaData> const*> mds)
    {
        std::size_t n = 0;
        for (auto const* md : mds) {
            if (*md) { n = std::max(n, (*md)->pattern.block_bytes); }
        }
        return n;
    }
 
    inline bool transposes_dense
        (std::initializer_list<std::unique_ptr<TransposeMetaData> const*> mds)
    {
        bool dense = true;
        for (auto const* md : mds) {
            if (*md) { dense = dense && (*md)->pattern.dense; }
        }
        return dense;
    }
 
    template <typename FA, typename DTOS>
    std::unique_ptr<TransposeMetaData>
    make_transpose (int n_active, FA& dst, Box const& dstbox, FA const& src,
                    IntVect const& ngrow, DTOS const& dtos)
    {
        return std::make_unique<TransposeMetaData>
            (transpose_pattern(dst, dstbox, src, dtos, n_active), dst, dstbox, src, ngrow, dtos);
    }

    template <typename FA, typename DTOS>
    void transpose (FA& dst, FA const& src, TransposeMetaData& md, int ncomp,
                    DTOS const& dtos, CommStrategy strategy, MPI_Comm comm)
    {
        BL_PROFILE("FFT::transpose");
        if (strategy == CommStrategy::pointtopoint) {
            NonLocalBC::ParallelCopy(dst, src, md.cmd, 0, 0, ncomp, dtos);
        } else if (comm == MPI_COMM_NULL) {
            return; // this process owns no data
        } else if (strategy == CommStrategy::alltoall && md.pattern.dense) {
            // With a sparse pattern the padded buffer would be much larger
            // than the data, so alltoallv is used instead.
            NonLocalBC::ParallelCopy_alltoall(dst, src, md.cmd, 0, 0, ncomp, dtos,
                                              NonLocalBC::Identity{}, md.pattern.block_bytes,
                                              &md.handler, comm);
        } else {
            NonLocalBC::ParallelCopy_alltoallv(dst, src, md.cmd, 0, 0, ncomp, dtos,
                                               NonLocalBC::Identity{}, &md.handler, comm);
        }
    }
 
    template <typename FA>
    typename FA::FABType::value_type * get_fab (FA& fa)
    {
        auto myproc = ParallelContext::MyProcSub();
        if (myproc < fa.size()) {
            return fa.fabPtr(myproc);
        } else {
            return nullptr;
        }
    }
 
    template <typename FA1, typename FA2>
    std::unique_ptr<char,DataDeleter> make_mfs_share (FA1& fa1, FA2& fa2)
    {
        bool not_same_fa = true;
        if constexpr (std::is_same_v<FA1,FA2>) {
            not_same_fa = (&fa1 != &fa2);
        }
        using FAB1 = typename FA1::FABType::value_type;
        using FAB2 = typename FA2::FABType::value_type;
        using T1 = typename FAB1::value_type;
        using T2 = typename FAB2::value_type;
        auto myproc = ParallelContext::MyProcSub();
        bool alloc_1 = (myproc < fa1.size());
        bool alloc_2 = (myproc < fa2.size()) && not_same_fa;
        void* p = nullptr;
        if (alloc_1 && alloc_2) {
            Box const& box1 = fa1.fabbox(myproc);
            Box const& box2 = fa2.fabbox(myproc);
            int ncomp1 = fa1.nComp();
            int ncomp2 = fa2.nComp();
            p = The_Arena()->alloc(std::max(sizeof(T1)*box1.numPts()*ncomp1,
                                            sizeof(T2)*box2.numPts()*ncomp2));
            fa1.setFab(myproc, FAB1(box1, ncomp1, (T1*)p));
            fa2.setFab(myproc, FAB2(box2, ncomp2, (T2*)p));
        } else if (alloc_1) {
            Box const& box1 = fa1.fabbox(myproc);
            int ncomp1 = fa1.nComp();
            p = The_Arena()->alloc(sizeof(T1)*box1.numPts()*ncomp1);
            fa1.setFab(myproc, FAB1(box1, ncomp1, (T1*)p));
        } else if (alloc_2) {
            Box const& box2 = fa2.fabbox(myproc);
            int ncomp2 = fa2.nComp();
            p = The_Arena()->alloc(sizeof(T2)*box2.numPts()*ncomp2);
            fa2.setFab(myproc, FAB2(box2, ncomp2, (T2*)p));
        } else {
            return nullptr;
        }
        return std::unique_ptr<char,DataDeleter>((char*)p, DataDeleter{The_Arena()});
    }
}

 
struct Swap01
{
    [[nodiscard]] constexpr Dim3 operator() (Dim3 i) const noexcept
    {
        return Dim3{.x = i.y, .y = i.x, .z = i.z};
    }
 
    static constexpr Dim3 Inverse (Dim3 i)
    {
        return Dim3{.x = i.y, .y = i.x, .z = i.z};
    }
 
    [[nodiscard]] constexpr IndexType operator() (IndexType it) const noexcept
    {
        return it;
    }
 
    static constexpr IndexType Inverse (IndexType it)
    {
        return it;
    }
};
 
struct Swap02
{
    [[nodiscard]] constexpr Dim3 operator() (Dim3 i) const noexcept
    {
        return Dim3{.x = i.z, .y = i.y, .z = i.x};
    }
 
    static constexpr Dim3 Inverse (Dim3 i)
    {
        return Dim3{.x = i.z, .y = i.y, .z = i.x};
    }
 
    [[nodiscard]] constexpr IndexType operator() (IndexType it) const noexcept
    {
        return it;
    }
 
    static constexpr IndexType Inverse (IndexType it)
    {
        return it;
    }
};
 
struct RotateFwd
{
    // dest -> src: (x,y,z) -> (y,z,x)
    [[nodiscard]] constexpr Dim3 operator() (Dim3 i) const noexcept
    {
        return Dim3{.x = i.y, .y = i.z, .z = i.x};
    }
 
    // src -> dest: (x,y,z) -> (z,x,y)
    static constexpr Dim3 Inverse (Dim3 i)
    {
        return Dim3{.x = i.z, .y = i.x, .z = i.y};
    }
 
    [[nodiscard]] constexpr IndexType operator() (IndexType it) const noexcept
    {
        return it;
    }
 
    static constexpr IndexType Inverse (IndexType it)
    {
        return it;
    }
};
 
struct RotateBwd
{
    // dest -> src: (x,y,z) -> (z,x,y)
    [[nodiscard]] constexpr Dim3 operator() (Dim3 i) const noexcept
    {
        return Dim3{.x = i.z, .y = i.x, .z = i.y};
    }
 
    // src -> dest: (x,y,z) -> (y,z,x)
    static constexpr Dim3 Inverse (Dim3 i)
    {
        return Dim3{.x = i.y, .y = i.z, .z = i.x};
    }
 
    [[nodiscard]] constexpr IndexType operator() (IndexType it) const noexcept
    {
        return it;
    }
 
    static constexpr IndexType Inverse (IndexType it)
    {
        return it;
    }
};


namespace detail
{
    struct SubHelper
    {
        explicit SubHelper (Box const& domain);
 
        [[nodiscard]] Box make_box (Box const& box) const;
 
        [[nodiscard]] Periodicity make_periodicity (Periodicity const& period) const;
 
        [[nodiscard]] bool ghost_safe (IntVect const& ng) const;
 
        // This rearranges the order.
        [[nodiscard]] IntVect make_iv (IntVect const& iv) const;
 
        // This keeps the order, but zero out the values in the hidden dimension.
        [[nodiscard]] IntVect make_safe_ghost (IntVect const& ng) const;
 
        [[nodiscard]] BoxArray inverse_boxarray (BoxArray const& ba) const;
 
        [[nodiscard]] IntVect inverse_order (IntVect const& order) const;
 
        template <typename T>
        [[nodiscard]] T make_array (T const& a) const
        {
#if (AMREX_SPACEDIM == 1)
            amrex::ignore_unused(this);
            return a;
#elif (AMREX_SPACEDIM == 2)
            if (m_case == case_1n) {
                return T{a[1],a[0]};
            } else {
                return a;
            }
#else
            if (m_case == case_11n) {
                return T{a[2],a[0],a[1]};
            } else if (m_case == case_1n1) {
                return T{a[1],a[0],a[2]};
            } else if (m_case == case_1nn) {
                return T{a[1],a[2],a[0]};
            } else if (m_case == case_n1n) {
                return T{a[0],a[2],a[1]};
            } else {
                return a;
            }
#endif
        }
 
        [[nodiscard]] GpuArray<int,3> xyz_order () const;
 
        template <typename FA>
        FA make_alias_mf (FA const& mf)
        {
            BoxList bl = mf.boxArray().boxList();
            for (auto& b : bl) {
                b = make_box(b);
            }
            auto const& ng = make_iv(mf.nGrowVect());
            FA submf(BoxArray(std::move(bl)), mf.DistributionMap(), mf.nComp(), ng, MFInfo{}.SetAlloc(false));
            using FAB = typename FA::fab_type;
            for (MFIter mfi(submf, MFItInfo().DisableDeviceSync()); mfi.isValid(); ++mfi) {
                submf.setFab(mfi, FAB(mfi.fabbox(), mf.nComp(), mf[mfi].dataPtr()));
            }
            return submf;
        }
 
#if (AMREX_SPACEDIM == 2)
        enum Case { case_1n, case_other };
        int m_case = case_other;
#elif (AMREX_SPACEDIM == 3)
        enum Case { case_11n, case_1n1, case_1nn, case_n1n, case_other };
        int m_case = case_other;
#endif
    };
}
 
}
 
#endif