diff --git a/Docs/sphinx_documentation/source/FFT.rst b/Docs/sphinx_documentation/source/FFT.rst
index 27dbb9ec2b..ac24e1c818 100644
--- a/Docs/sphinx_documentation/source/FFT.rst
+++ b/Docs/sphinx_documentation/source/FFT.rst
@@ -93,6 +93,14 @@ in an :cpp:`FFT::Info` object passed to the constructor of
 
     r2c.backward(cmf, mf);
 
+.. _sec:FFT:c2c:
+
+FFT::C2C Class
+==============
+
+:cpp:`FFT::C2C` is a class template that supports complex to complex Fourier
+transforms. It has a similar interface as :cpp:`FFT::R2C`.
+
 .. _sec:FFT:localr2c:
 
 FFT::LocalR2C Class
diff --git a/Src/FFT/AMReX_FFT_Helper.H b/Src/FFT/AMReX_FFT_Helper.H
index 1c41cbffff..c8d6f2babe 100644
--- a/Src/FFT/AMReX_FFT_Helper.H
+++ b/Src/FFT/AMReX_FFT_Helper.H
@@ -61,7 +61,7 @@ struct Info
 
     //! For automatic strategy, this is the size per process below which we
     //! switch from slab to pencil.
-    int pencil_threshold = 8;
+    int pencil_threshold = 4;
 
     //! Supported only in 3D. When twod_mode is true, FFT is performed on
     //! the first two dimensions only and the third dimension size is the
@@ -310,7 +310,7 @@ struct Plan
     void init_r2c (IntVectND<M> const& fft_size, void*, void*, bool cache, int ncomp = 1);
 
     template <Direction D>
-    void init_c2c (Box const& box, VendorComplex* p, int ncomp = 1)
+    void init_c2c (Box const& box, VendorComplex* p, int ncomp = 1, int ndims = 1)
     {
         static_assert(D == Direction::forward || D == Direction::backward);
 
@@ -319,9 +319,35 @@ struct Plan
         pf = (void*)p;
         pb = (void*)p;
 
-        n = box.length(0);
-        howmany = AMREX_D_TERM(1, *box.length(1), *box.length(2));
-        howmany *= ncomp;
+        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));
@@ -330,7 +356,7 @@ struct Plan
         cufftType t = std::is_same_v<float,T> ? CUFFT_C2C : CUFFT_Z2Z;
         std::size_t work_size;
         AMREX_CUFFT_SAFE_CALL
-            (cufftMakePlanMany(plan, 1, &n, nullptr, 1, n, nullptr, 1, n, t, howmany, &work_size));
+            (cufftMakePlanMany(plan, ndims, len, nullptr, 1, n, nullptr, 1, n, t, howmany, &work_size));
 
 #elif defined(AMREX_USE_HIP)
 
@@ -338,14 +364,31 @@ struct Plan
                                             : rocfft_precision_double;
         auto dir= (D == Direction::forward) ? rocfft_transform_type_complex_forward
                                             : rocfft_transform_type_complex_inverse;
-        const std::size_t length = n;
+        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, 1,
-                                &length, howmany, nullptr));
+            (rocfft_plan_create(&plan, rocfft_placement_inplace, dir, prec, ndims,
+                                length, howmany, nullptr));
 
 #elif defined(AMREX_USE_SYCL)
 
-        auto* pp = new mkl_desc_c(n);
+        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);
@@ -355,7 +398,12 @@ struct Plan
         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 = {0,1};
+        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[ndims-1-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);
@@ -373,21 +421,21 @@ struct Plan
         if constexpr (std::is_same_v<float,T>) {
             if constexpr (D == Direction::forward) {
                 plan = fftwf_plan_many_dft
-                    (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
+                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
                      FFTW_ESTIMATE);
             } else {
                 plan = fftwf_plan_many_dft
-                    (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, +1,
+                    (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
-                    (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
+                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, -1,
                      FFTW_ESTIMATE);
             } else {
                 plan = fftw_plan_many_dft
-                    (1, &n, howmany, p, nullptr, 1, n, p, nullptr, 1, n, +1,
+                    (ndims, len, howmany, p, nullptr, 1, n, p, nullptr, 1, n, +1,
                      FFTW_ESTIMATE);
             }
         }
diff --git a/Src/FFT/AMReX_FFT_R2C.H b/Src/FFT/AMReX_FFT_R2C.H
index 8d1837e971..36db327d90 100644
--- a/Src/FFT/AMReX_FFT_R2C.H
+++ b/Src/FFT/AMReX_FFT_R2C.H
@@ -29,13 +29,14 @@ template <typename T> class PoissonHybrid;
  * For more details, we refer the users to
  * https://amrex-codes.github.io/amrex/docs_html/FFT_Chapter.html.
  */
-template <typename T = Real, FFT::Direction D = FFT::Direction::both>
+template <typename T = Real, FFT::Direction D = FFT::Direction::both, bool C = false>
 class R2C
 {
 public:
-    using MF = std::conditional_t<std::is_same_v<T,Real>,
-                                  MultiFab, FabArray<BaseFab<T> > >;
     using cMF = FabArray<BaseFab<GpuComplex<T> > >;
+    using MF = std::conditional_t
+        <C, cMF, std::conditional_t<std::is_same_v<T,Real>,
+                                    MultiFab, FabArray<BaseFab<T> > >>;
 
     template <typename U> friend class OpenBCSolver;
     template <typename U> friend class Poisson;
@@ -183,7 +184,52 @@ private:
                         Periodicity const& period = Periodicity::NonPeriodic(),
                         int incomp = 0, int outcomp = 0);
 
-    std::pair<Plan<T>,Plan<T>> make_c2c_plans (cMF& inout) const;
+    std::pair<Plan<T>,Plan<T>> make_c2c_plans (cMF& inout, int ndims) const;
+
+    static Box make_domain_x (Box const& domain)
+    {
+        if constexpr (C) {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(0)-1,
+                                                        domain.length(1)-1,
+                                                        domain.length(2)-1)),
+                       domain.ixType());
+        } else {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(0)/2,
+                                                        domain.length(1)-1,
+                                                        domain.length(2)-1)),
+                       domain.ixType());
+        }
+    }
+
+    static Box make_domain_y (Box const& domain)
+    {
+        if constexpr (C) {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(1)-1,
+                                                        domain.length(0)-1,
+                                                        domain.length(2)-1)),
+                       domain.ixType());
+        } else {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(1)-1,
+                                                        domain.length(0)/2,
+                                                        domain.length(2)-1)),
+                       domain.ixType());
+        }
+    }
+
+    static Box make_domain_z (Box const& domain)
+    {
+        if constexpr (C) {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(2)-1,
+                                                        domain.length(0)-1,
+                                                        domain.length(1)-1)),
+                       domain.ixType());
+        } else {
+            return Box(IntVect(0), IntVect(AMREX_D_DECL(domain.length(2)-1,
+                                                        domain.length(0)/2,
+                                                        domain.length(1)-1)),
+                       domain.ixType());
+        }
+    }
 
     Plan<T> m_fft_fwd_x{};
     Plan<T> m_fft_bwd_x{};
@@ -225,7 +271,7 @@ private:
     Box m_spectral_domain_y;
     Box m_spectral_domain_z;
 
-    std::unique_ptr<R2C<T,D>> m_r2c_sub;
+    std::unique_ptr<R2C<T,D,C>> m_r2c_sub;
     detail::SubHelper m_sub_helper;
 
     Info m_info;
@@ -235,23 +281,14 @@ private:
     bool m_openbc_half = false;
 };
 
-template <typename T, Direction D>
-R2C<T,D>::R2C (Box const& domain, Info const& info)
+template <typename T, Direction D, bool C>
+R2C<T,D,C>::R2C (Box const& domain, Info const& info)
     : m_real_domain(domain),
-      m_spectral_domain_x(IntVect(0), IntVect(AMREX_D_DECL(domain.length(0)/2,
-                                                           domain.length(1)-1,
-                                                           domain.length(2)-1)),
-                          domain.ixType()),
+      m_spectral_domain_x(make_domain_x(domain)),
 #if (AMREX_SPACEDIM >= 2)
-      m_spectral_domain_y(IntVect(0), IntVect(AMREX_D_DECL(domain.length(1)-1,
-                                                           domain.length(0)/2,
-                                                           domain.length(2)-1)),
-                          domain.ixType()),
+      m_spectral_domain_y(make_domain_y(domain)),
 #if (AMREX_SPACEDIM == 3)
-      m_spectral_domain_z(IntVect(0), IntVect(AMREX_D_DECL(domain.length(2)-1,
-                                                           domain.length(0)/2,
-                                                           domain.length(1)-1)),
-                          domain.ixType()),
+      m_spectral_domain_z(make_domain_z(domain)),
 #endif
 #endif
       m_sub_helper(domain),
@@ -275,7 +312,7 @@ R2C<T,D>::R2C (Box const& domain, Info const& info)
     {
         Box subbox = m_sub_helper.make_box(m_real_domain);
         if (subbox.size() != m_real_domain.size()) {
-            m_r2c_sub = std::make_unique<R2C<T,D>>(subbox, m_info);
+            m_r2c_sub = std::make_unique<R2C<T,D,C>>(subbox, m_info);
             return;
         }
     }
@@ -285,11 +322,19 @@ R2C<T,D>::R2C (Box const& domain, Info const& info)
 
 #if (AMREX_SPACEDIM == 3)
     if (m_info.domain_strategy == DomainStrategy::automatic) {
-        int shortside = m_real_domain.shortside();
-        if (shortside <= m_info.pencil_threshold*nprocs) {
-            m_info.domain_strategy = DomainStrategy::pencil;
+        if (m_info.twod_mode) {
+            if (m_real_domain.length(2) < nprocs) {
+                m_info.domain_strategy = DomainStrategy::pencil;
+            } else {
+                m_info.domain_strategy = DomainStrategy::slab;
+            }
         } else {
-            m_info.domain_strategy = DomainStrategy::slab;
+            int shortside = m_real_domain.shortside();
+            if (shortside < m_info.pencil_threshold*nprocs) {
+                m_info.domain_strategy = DomainStrategy::pencil;
+            } else {
+                m_info.domain_strategy = DomainStrategy::slab;
+            }
         }
     }
     if (m_info.domain_strategy == DomainStrategy::slab && (m_real_domain.length(1) > 1)) {
@@ -364,12 +409,28 @@ R2C<T,D>::R2C (Box const& domain, Info const& info)
         }
 #endif
 
-        if (m_slab_decomp) {
-            m_data_1 = detail::make_mfs_share(m_rx, m_cz);
-            m_data_2 = detail::make_mfs_share(m_cx, m_cx);
+        if constexpr (C) {
+            if (m_slab_decomp) {
+                m_data_1 = detail::make_mfs_share(m_rx, m_cx);
+                m_data_2 = detail::make_mfs_share(m_cz, m_cz);
+            } else {
+                m_data_1 = detail::make_mfs_share(m_rx, m_cz);
+                m_data_2 = detail::make_mfs_share(m_cy, m_cy);
+                // make m_cx an alias to m_rx
+                if (myproc < m_cx.size()) {
+                    Box const& box = m_cx.fabbox(myproc);
+                    using FAB = typename cMF::FABType::value_type;
+                    m_cx.setFab(myproc, FAB(box, ncomp, m_rx[myproc].dataPtr()));
+                }
+            }
         } else {
-            m_data_1 = detail::make_mfs_share(m_rx, m_cy);
-            m_data_2 = detail::make_mfs_share(m_cx, m_cz);
+            if (m_slab_decomp) {
+                m_data_1 = detail::make_mfs_share(m_rx, m_cz);
+                m_data_2 = detail::make_mfs_share(m_cx, m_cx);
+            } else {
+                m_data_1 = detail::make_mfs_share(m_rx, m_cy);
+                m_data_2 = detail::make_mfs_share(m_cx, m_cz);
+            }
         }
 
         //
@@ -409,57 +470,67 @@ R2C<T,D>::R2C (Box const& domain, Info const& info)
 
         if (myproc < m_rx.size())
         {
-            Box const& box = m_rx.box(myproc);
-            auto* pr = m_rx[myproc].dataPtr();
-            auto* pc = (typename Plan<T>::VendorComplex *)m_cx[myproc].dataPtr();
+            if constexpr (C) {
+                int ndims = m_slab_decomp ? 2 : 1;
+                std::tie(m_fft_fwd_x, m_fft_bwd_x) = make_c2c_plans(m_cx, ndims);
+            } else {
+                Box const& box = m_rx.box(myproc);
+                auto* pr = m_rx[myproc].dataPtr();
+                auto* pc = (typename Plan<T>::VendorComplex *)m_cx[myproc].dataPtr();
 #ifdef AMREX_USE_SYCL
-            m_fft_fwd_x.template init_r2c<Direction::forward>(box, pr, pc, m_slab_decomp, ncomp);
-            m_fft_bwd_x = m_fft_fwd_x;
-#else
-            if constexpr (D == Direction::both || D == Direction::forward) {
                 m_fft_fwd_x.template init_r2c<Direction::forward>(box, pr, pc, m_slab_decomp, ncomp);
-            }
-            if constexpr (D == Direction::both || D == Direction::backward) {
-                m_fft_bwd_x.template init_r2c<Direction::backward>(box, pr, pc, m_slab_decomp, ncomp);
-            }
+                m_fft_bwd_x = m_fft_fwd_x;
+#else
+                if constexpr (D == Direction::both || D == Direction::forward) {
+                    m_fft_fwd_x.template init_r2c<Direction::forward>(box, pr, pc, m_slab_decomp, ncomp);
+                }
+                if constexpr (D == Direction::both || D == Direction::backward) {
+                    m_fft_bwd_x.template init_r2c<Direction::backward>(box, pr, pc, m_slab_decomp, ncomp);
+                }
 #endif
+            }
         }
 
 #if (AMREX_SPACEDIM >= 2)
         if (! m_cy.empty()) {
-            std::tie(m_fft_fwd_y, m_fft_bwd_y) = make_c2c_plans(m_cy);
+            std::tie(m_fft_fwd_y, m_fft_bwd_y) = make_c2c_plans(m_cy,1);
         }
 #endif
 #if (AMREX_SPACEDIM == 3)
         if (! m_cz.empty()) {
-            std::tie(m_fft_fwd_z, m_fft_bwd_z) = make_c2c_plans(m_cz);
+            std::tie(m_fft_fwd_z, m_fft_bwd_z) = make_c2c_plans(m_cz,1);
         }
 #endif
     }
     else // do fft in all dimensions at the same time
     {
-        m_data_1 = detail::make_mfs_share(m_rx, m_rx);
-        m_data_2 = detail::make_mfs_share(m_cx, m_cx);
+        if constexpr (C) {
+            m_data_1 = detail::make_mfs_share(m_rx, m_cx);
+            std::tie(m_fft_fwd_x, m_fft_bwd_x) = make_c2c_plans(m_cx,AMREX_SPACEDIM);
+        } else {
+            m_data_1 = detail::make_mfs_share(m_rx, m_rx);
+            m_data_2 = detail::make_mfs_share(m_cx, m_cx);
 
-        auto const& len = m_real_domain.length();
-        auto* pr = (void*)m_rx[0].dataPtr();
-        auto* pc = (void*)m_cx[0].dataPtr();
+            auto const& len = m_real_domain.length();
+            auto* pr = (void*)m_rx[0].dataPtr();
+            auto* pc = (void*)m_cx[0].dataPtr();
 #ifdef AMREX_USE_SYCL
-        m_fft_fwd_x.template init_r2c<Direction::forward>(len, pr, pc, false, ncomp);
-        m_fft_bwd_x = m_fft_fwd_x;
-#else
-        if constexpr (D == Direction::both || D == Direction::forward) {
             m_fft_fwd_x.template init_r2c<Direction::forward>(len, pr, pc, false, ncomp);
-        }
-        if constexpr (D == Direction::both || D == Direction::backward) {
-            m_fft_bwd_x.template init_r2c<Direction::backward>(len, pr, pc, false, ncomp);
-        }
+            m_fft_bwd_x = m_fft_fwd_x;
+#else
+            if constexpr (D == Direction::both || D == Direction::forward) {
+                m_fft_fwd_x.template init_r2c<Direction::forward>(len, pr, pc, false, ncomp);
+            }
+            if constexpr (D == Direction::both || D == Direction::backward) {
+                m_fft_bwd_x.template init_r2c<Direction::backward>(len, pr, pc, false, ncomp);
+            }
 #endif
+        }
     }
 }
 
-template <typename T, Direction D>
-R2C<T,D>::~R2C<T,D> ()
+template <typename T, Direction D, bool C>
+R2C<T,D,C>::~R2C<T,D,C> ()
 {
     if (m_fft_bwd_x.plan != m_fft_fwd_x.plan) {
         m_fft_bwd_x.destroy();
@@ -479,10 +550,10 @@ R2C<T,D>::~R2C<T,D> ()
     m_fft_fwd_x_half.destroy();
 }
 
-template <typename T, Direction D>
-void R2C<T,D>::prepare_openbc ()
+template <typename T, Direction D, bool C>
+void R2C<T,D,C>::prepare_openbc ()
 {
-    if (m_r2c_sub) { amrex::Abort("R2C: OpenBC not supported with reduced dimensions"); }
+    if (C || m_r2c_sub) { amrex::Abort("R2C: OpenBC not supported with reduced dimensions or complex inputs"); }
 
 #if (AMREX_SPACEDIM == 3)
     if (m_do_alld_fft) { return; }
@@ -535,10 +606,10 @@ void R2C<T,D>::prepare_openbc ()
 #endif
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <Direction DIR, std::enable_if_t<DIR == Direction::forward ||
                                           DIR == Direction::both, int> >
-void R2C<T,D>::forward (MF const& inmf, int incomp)
+void R2C<T,D,C>::forward (MF const& inmf, int incomp)
 {
     BL_PROFILE("FFT::R2C::forward(in)");
 
@@ -560,12 +631,20 @@ void R2C<T,D>::forward (MF const& inmf, int incomp)
     }
 
     if (m_do_alld_fft) {
-        m_fft_fwd_x.template compute_r2c<Direction::forward>();
+        if constexpr (C) {
+            m_fft_fwd_x.template compute_c2c<Direction::forward>();
+        } else {
+            m_fft_fwd_x.template compute_r2c<Direction::forward>();
+        }
         return;
     }
 
     auto& fft_x = m_openbc_half ? m_fft_fwd_x_half : m_fft_fwd_x;
-    fft_x.template compute_r2c<Direction::forward>();
+    if constexpr (C) {
+        fft_x.template compute_c2c<Direction::forward>();
+    } else {
+        fft_x.template compute_r2c<Direction::forward>();
+    }
 
     if (                          m_cmd_x2y) {
         ParallelCopy(m_cy, m_cx, *m_cmd_x2y, 0, 0, ncomp, m_dtos_x2y);
@@ -599,15 +678,15 @@ void R2C<T,D>::forward (MF const& inmf, int incomp)
     m_fft_fwd_z.template compute_c2c<Direction::forward>();
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <Direction DIR, std::enable_if_t<DIR == Direction::both, int> >
-void R2C<T,D>::backward (MF& outmf, int outcomp)
+void R2C<T,D,C>::backward (MF& outmf, int outcomp)
 {
     backward_doit(outmf, IntVect(0), Periodicity::NonPeriodic(), outcomp);
 }
 
-template <typename T, Direction D>
-void R2C<T,D>::backward_doit (MF& outmf, IntVect const& ngout,
+template <typename T, Direction D, bool C>
+void R2C<T,D,C>::backward_doit (MF& outmf, IntVect const& ngout,
                               Periodicity const& period, int outcomp)
 {
     BL_PROFILE("FFT::R2C::backward(out)");
@@ -630,7 +709,11 @@ void R2C<T,D>::backward_doit (MF& outmf, IntVect const& ngout,
     }
 
     if (m_do_alld_fft) {
-        m_fft_bwd_x.template compute_r2c<Direction::backward>();
+        if constexpr (C) {
+            m_fft_bwd_x.template compute_c2c<Direction::backward>();
+        } else {
+            m_fft_bwd_x.template compute_r2c<Direction::backward>();
+        }
         outmf.ParallelCopy(m_rx, 0, outcomp, ncomp, IntVect(0),
                            amrex::elemwiseMin(ngout,outmf.nGrowVect()), period);
         return;
@@ -653,14 +736,18 @@ void R2C<T,D>::backward_doit (MF& outmf, IntVect const& ngout,
     }
 
     auto& fft_x = m_openbc_half ? m_fft_bwd_x_half : m_fft_bwd_x;
-    fft_x.template compute_r2c<Direction::backward>();
+    if constexpr (C) {
+        fft_x.template compute_c2c<Direction::backward>();
+    } else {
+        fft_x.template compute_r2c<Direction::backward>();
+    }
     outmf.ParallelCopy(m_rx, 0, outcomp, ncomp, IntVect(0),
                        amrex::elemwiseMin(ngout,outmf.nGrowVect()), period);
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 std::pair<Plan<T>, Plan<T>>
-R2C<T,D>::make_c2c_plans (cMF& inout) const
+R2C<T,D,C>::make_c2c_plans (cMF& inout, int ndims) const
 {
     Plan<T> fwd;
     Plan<T> bwd;
@@ -674,23 +761,23 @@ R2C<T,D>::make_c2c_plans (cMF& inout) const
     auto const ncomp = m_info.batch_size;
 
 #ifdef AMREX_USE_SYCL
-    fwd.template init_c2c<Direction::forward>(box, pio, ncomp);
+    fwd.template init_c2c<Direction::forward>(box, pio, ncomp, ndims);
     bwd = fwd;
 #else
     if constexpr (D == Direction::both || D == Direction::forward) {
-        fwd.template init_c2c<Direction::forward>(box, pio, ncomp);
+        fwd.template init_c2c<Direction::forward>(box, pio, ncomp, ndims);
     }
     if constexpr (D == Direction::both || D == Direction::backward) {
-        bwd.template init_c2c<Direction::backward>(box, pio, ncomp);
+        bwd.template init_c2c<Direction::backward>(box, pio, ncomp, ndims);
     }
 #endif
 
     return {fwd, bwd};
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <typename F>
-void R2C<T,D>::post_forward_doit_0 (F const& post_forward)
+void R2C<T,D,C>::post_forward_doit_0 (F const& post_forward)
 {
     if (m_info.twod_mode || m_info.batch_size > 1) {
         amrex::Abort("xxxxx todo: post_forward");
@@ -743,9 +830,9 @@ void R2C<T,D>::post_forward_doit_0 (F const& post_forward)
     }
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <typename F>
-void R2C<T,D>::post_forward_doit_1 (F const& post_forward)
+void R2C<T,D,C>::post_forward_doit_1 (F const& post_forward)
 {
     if (m_info.twod_mode || m_info.batch_size > 1) {
         amrex::Abort("xxxxx todo: post_forward");
@@ -786,8 +873,8 @@ void R2C<T,D>::post_forward_doit_1 (F const& post_forward)
     }
 }
 
-template <typename T, Direction D>
-T R2C<T,D>::scalingFactor () const
+template <typename T, Direction D, bool C>
+T R2C<T,D,C>::scalingFactor () const
 {
 #if (AMREX_SPACEDIM == 3)
     if (m_info.twod_mode) {
@@ -804,11 +891,11 @@ T R2C<T,D>::scalingFactor () const
     }
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <Direction DIR, std::enable_if_t<DIR == Direction::forward ||
                                           DIR == Direction::both, int> >
-std::pair<typename R2C<T,D>::cMF *, IntVect>
-R2C<T,D>::getSpectralData ()
+std::pair<typename R2C<T,D,C>::cMF *, IntVect>
+R2C<T,D,C>::getSpectralData ()
 {
 #if (AMREX_SPACEDIM > 1)
     if (m_r2c_sub) {
@@ -825,10 +912,10 @@ R2C<T,D>::getSpectralData ()
     }
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <Direction DIR, std::enable_if_t<DIR == Direction::forward ||
                                           DIR == Direction::both, int> >
-void R2C<T,D>::forward (MF const& inmf, cMF& outmf, int incomp, int outcomp)
+void R2C<T,D,C>::forward (MF const& inmf, cMF& outmf, int incomp, int outcomp)
 {
     BL_PROFILE("FFT::R2C::forward(inout)");
 
@@ -885,16 +972,16 @@ void R2C<T,D>::forward (MF const& inmf, cMF& outmf, int incomp, int outcomp)
     }
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 template <Direction DIR, std::enable_if_t<DIR == Direction::backward ||
                                           DIR == Direction::both, int> >
-void R2C<T,D>::backward (cMF const& inmf, MF& outmf, int incomp, int outcomp)
+void R2C<T,D,C>::backward (cMF const& inmf, MF& outmf, int incomp, int outcomp)
 {
     backward_doit(inmf, outmf, IntVect(0), Periodicity::NonPeriodic(), incomp, outcomp);
 }
 
-template <typename T, Direction D>
-void R2C<T,D>::backward_doit (cMF const& inmf, MF& outmf, IntVect const& ngout,
+template <typename T, Direction D, bool C>
+void R2C<T,D,C>::backward_doit (cMF const& inmf, MF& outmf, IntVect const& ngout,
                               Periodicity const& period, int incomp, int outcomp)
 {
     BL_PROFILE("FFT::R2C::backward(inout)");
@@ -955,9 +1042,9 @@ void R2C<T,D>::backward_doit (cMF const& inmf, MF& outmf, IntVect const& ngout,
     }
 }
 
-template <typename T, Direction D>
+template <typename T, Direction D, bool C>
 std::pair<BoxArray,DistributionMapping>
-R2C<T,D>::getSpectralDataLayout () const
+R2C<T,D,C>::getSpectralDataLayout () const
 {
 #if (AMREX_SPACEDIM > 1)
     if (m_r2c_sub) {
@@ -1001,6 +1088,10 @@ R2C<T,D>::getSpectralDataLayout () const
     }
 }
 
+//! FFT between complex data.
+template <typename T = Real, FFT::Direction D = FFT::Direction::both>
+using C2C = R2C<T, D, true>;
+
 }
 
 #endif
diff --git a/Tests/FFT/C2C/CMakeLists.txt b/Tests/FFT/C2C/CMakeLists.txt
new file mode 100644
index 0000000000..21a9d3b268
--- /dev/null
+++ b/Tests/FFT/C2C/CMakeLists.txt
@@ -0,0 +1,10 @@
+foreach(D IN LISTS AMReX_SPACEDIM)
+    set(_sources  main.cpp)
+
+    set(_input_files)
+
+    setup_test(${D} _sources _input_files)
+
+    unset(_sources)
+    unset(_input_files)
+endforeach()
diff --git a/Tests/FFT/C2C/GNUmakefile b/Tests/FFT/C2C/GNUmakefile
new file mode 100644
index 0000000000..93376f4485
--- /dev/null
+++ b/Tests/FFT/C2C/GNUmakefile
@@ -0,0 +1,26 @@
+AMREX_HOME := ../../..
+
+DEBUG	= FALSE
+
+DIM	= 3
+
+COMP    = gcc
+
+USE_MPI   = TRUE
+USE_OMP   = FALSE
+USE_CUDA  = FALSE
+USE_HIP   = FALSE
+USE_SYCL  = FALSE
+
+USE_FFT = TRUE
+
+BL_NO_FORT = TRUE
+
+TINY_PROFILE = FALSE
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ./Make.package
+include $(AMREX_HOME)/Src/Base/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
diff --git a/Tests/FFT/C2C/Make.package b/Tests/FFT/C2C/Make.package
new file mode 100644
index 0000000000..6b4b865e8f
--- /dev/null
+++ b/Tests/FFT/C2C/Make.package
@@ -0,0 +1 @@
+CEXE_sources += main.cpp
diff --git a/Tests/FFT/C2C/main.cpp b/Tests/FFT/C2C/main.cpp
new file mode 100644
index 0000000000..df87d4e641
--- /dev/null
+++ b/Tests/FFT/C2C/main.cpp
@@ -0,0 +1,134 @@
+#include <AMReX_FFT.H> // Put this at the top for testing
+
+#include <AMReX.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_PlotFileUtil.H>
+
+using namespace amrex;
+
+int main (int argc, char* argv[])
+{
+    amrex::Initialize(argc, argv);
+    {
+        BL_PROFILE("main");
+
+        AMREX_D_TERM(int n_cell_x = 128;,
+                     int n_cell_y = 32;,
+                     int n_cell_z = 64);
+
+        AMREX_D_TERM(int max_grid_size_x = 16;,
+                     int max_grid_size_y = 32;,
+                     int max_grid_size_z = 32);
+
+        AMREX_D_TERM(Real prob_lo_x = 0.;,
+                     Real prob_lo_y = 0.;,
+                     Real prob_lo_z = 0.);
+        AMREX_D_TERM(Real prob_hi_x = 1.;,
+                     Real prob_hi_y = 1.;,
+                     Real prob_hi_z = 1.);
+
+        {
+            ParmParse pp;
+            AMREX_D_TERM(pp.query("n_cell_x", n_cell_x);,
+                         pp.query("n_cell_y", n_cell_y);,
+                         pp.query("n_cell_z", n_cell_z));
+            AMREX_D_TERM(pp.query("max_grid_size_x", max_grid_size_x);,
+                         pp.query("max_grid_size_y", max_grid_size_y);,
+                         pp.query("max_grid_size_z", max_grid_size_z));
+        }
+
+        Box domain(IntVect(0),IntVect(AMREX_D_DECL(n_cell_x-1,n_cell_y-1,n_cell_z-1)));
+        BoxArray ba(domain);
+        ba.maxSize(IntVect(AMREX_D_DECL(max_grid_size_x,
+                                        max_grid_size_y,
+                                        max_grid_size_z)));
+        DistributionMapping dm(ba);
+
+        Geometry geom;
+        {
+            geom.define(domain,
+                        RealBox(AMREX_D_DECL(prob_lo_x,prob_lo_y,prob_lo_z),
+                                AMREX_D_DECL(prob_hi_x,prob_hi_y,prob_hi_z)),
+                        CoordSys::cartesian, {AMREX_D_DECL(1,1,1)});
+        }
+        auto const& dx = geom.CellSizeArray();
+
+        cMultiFab mf(ba,dm,1,0);
+        auto const& ma = mf.arrays();
+        ParallelFor(mf, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+        {
+            AMREX_D_TERM(Real x = (i+0.5_rt) * dx[0] - 0.5_rt;,
+                         Real y = (j+0.5_rt) * dx[1] - 0.5_rt;,
+                         Real z = (k+0.5_rt) * dx[2] - 0.5_rt);
+            auto tmp = std::exp(-10._rt*
+                (AMREX_D_TERM(x*x*1.05_rt, + y*y*0.90_rt, + z*z)));
+            ma[b](i,j,k) = GpuComplex<Real>{tmp, 1._rt-tmp};
+        });
+
+        cMultiFab mf2(ba,dm,1,0);
+
+        auto scaling = Real(1) / Real(geom.Domain().d_numPts());
+
+        {
+            FFT::Info info{};
+            info.setDomainStrategy(FFT::DomainStrategy::slab);
+            FFT::C2C<Real,FFT::Direction::both> c2c(geom.Domain(), info);
+            c2c.forward(mf);
+            c2c.backward(mf2);
+        }
+
+        {
+            MultiFab errmf(ba,dm,1,0);
+            auto const& errma = errmf.arrays();
+
+            auto const& ma2 = mf2.const_arrays();
+            ParallelFor(mf2, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+            {
+                auto err = ma[b](i,j,k) - ma2[b](i,j,k)*scaling;
+                errma[b](i,j,k) = amrex::norm(err);
+            });
+
+            auto error = errmf.norminf();
+            amrex::Print() << "  Expected to be close to zero: " << error << "\n";
+#ifdef AMREX_USE_FLOAT
+            auto eps = 1.e-6f;
+#else
+            auto eps = 1.e-13;
+#endif
+            AMREX_ALWAYS_ASSERT(error < eps);
+        }
+
+        mf2.setVal(0);
+
+        {
+            FFT::Info info{};
+            info.setDomainStrategy(FFT::DomainStrategy::pencil);
+            FFT::C2C<Real,FFT::Direction::both> c2c(geom.Domain(), info);
+            c2c.forward(mf);
+            c2c.backward(mf2);
+        }
+
+        {
+            MultiFab errmf(ba,dm,1,0);
+            auto const& errma = errmf.arrays();
+
+            auto const& ma2 = mf2.const_arrays();
+            ParallelFor(mf2, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
+            {
+                auto err = ma[b](i,j,k) - ma2[b](i,j,k)*scaling;
+                errma[b](i,j,k) = amrex::norm(err);
+            });
+
+            auto error = errmf.norminf();
+            amrex::Print() << "  Expected to be close to zero: " << error << "\n";
+#ifdef AMREX_USE_FLOAT
+            auto eps = 1.e-6f;
+#else
+            auto eps = 1.e-13;
+#endif
+            AMREX_ALWAYS_ASSERT(error < eps);
+        }
+    }
+    amrex::Finalize();
+}