Add DUCC FFT support

unsupported/Eigen/FFT

@@ -37,12 +37,13 @@
  * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
  * - MKL (https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html) : fastest, free -- may be
  * incompatible with Eigen in GPL form.
- * - pocketfft (https://gitlab.mpcdf.mpg.de/mtr/pocketfft) : faster than kissfft, BSD 3-clause.
+ * - PocketFFT/DUCC (https://gitlab.mpcdf.mpg.de/mtr/pocketfft, https://gitlab.mpcdf.mpg.de/mtr/ducc) : faster than kissfft, BSD 3-clause.
  *   It is a heavily modified implementation of FFTPack, with the following advantages:
  *   1.strictly C++11 compliant
  *   2.more accurate twiddle factor computation
  *   3.very fast plan generation
  *   4.worst case complexity for transform sizes with large prime factors is N*log(N), because Bluestein's algorithm is
+ *   According to the author, DUCC contains the "evolution" of pocketfft, though the interface is very similar.
  * used for these cases
  *
  * \section FFTDesign Design
@@ -85,7 +86,7 @@
 #ifdef EIGEN_FFTW_DEFAULT
 // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
 #include <fftw3.h>
-#include "src/FFT/ei_fftw_impl.h"
+#include "src/FFT/fftw_impl.h"
 namespace Eigen {
 // template <typename T> typedef struct internal::fftw_impl  default_fft_impl; this does not work
 template <typename T>
@@ -93,7 +94,7 @@ struct default_fft_impl : public internal::fftw_impl<T> {};
 }  // namespace Eigen
 #elif defined EIGEN_MKL_DEFAULT
 // intel Math Kernel Library: fastest, free -- may be incompatible with Eigen in GPL form
-#include "src/FFT/ei_imklfft_impl.h"
+#include "src/FFT/imklfft_impl.h"
 namespace Eigen {
 template <typename T>
 struct default_fft_impl : public internal::imklfft::imklfft_impl<T> {};
@@ -101,14 +102,24 @@ struct default_fft_impl : public internal::imklfft::imklfft_impl<T> {};
 #elif defined EIGEN_POCKETFFT_DEFAULT
 // internal::pocketfft_impl: a heavily modified implementation of FFTPack, with many advantages.
 #include <pocketfft_hdronly.h>
-#include "src/FFT/ei_pocketfft_impl.h"
+#include "src/FFT/pocketfft_impl.h"
 namespace Eigen {
 template <typename T>
 struct default_fft_impl : public internal::pocketfft_impl<T> {};
 }  // namespace Eigen
+#elif defined EIGEN_DUCCFFT_DEFAULT
+#include <ducc0/fft/fft.h>
+#include <ducc0/infra/string_utils.h>
+#include <ducc0/fft/fft.h>
+#include <ducc0/fft/fftnd_impl.h>
+#include "src/FFT/duccfft_impl.h"
+namespace Eigen {
+template <typename T>
+struct default_fft_impl : public internal::duccfft_impl<T> {};
+}  // namespace Eigen
 #else
 // internal::kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
-#include "src/FFT/ei_kissfft_impl.h"
+#include "src/FFT/kissfft_impl.h"
 namespace Eigen {
 template <typename T>
 struct default_fft_impl : public internal::kissfft_impl<T> {};
@@ -204,7 +215,8 @@ class FFT {
 
   inline void fwd(Complex* dst, const Complex* src, Index nfft) { m_impl.fwd(dst, src, static_cast<int>(nfft)); }
 
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+    defined EIGEN_MKL_DEFAULT
   inline void fwd2(Complex* dst, const Complex* src, int n0, int n1) { m_impl.fwd2(dst, src, n0, n1); }
 #endif
 
@@ -366,7 +378,8 @@ class FFT {
     inv(&dst[0], &src[0], nfft);
   }
 
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+    defined EIGEN_MKL_DEFAULT
   inline void inv2(Complex* dst, const Complex* src, int n0, int n1) {
     m_impl.inv2(dst, src, n0, n1);
     if (HasFlag(Unscaled) == false) scale(dst, 1. / (n0 * n1), n0 * n1);
@@ -385,7 +398,6 @@ class FFT {
       Matrix<T_Data, Dynamic, 1>::Map(x, nx) *= s;
     else
       Matrix<T_Data, Dynamic, 1>::MapAligned(x, nx) *= s;
-      // Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
 #endif
   }
 

unsupported/Eigen/src/FFT/duccfft_impl.h

@@ -0,0 +1,71 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+namespace Eigen {
+
+namespace internal {
+
+template <typename _Scalar>
+struct duccfft_impl {
+  using Scalar = _Scalar;
+  using Complex = std::complex<Scalar>;
+  using shape_t = ducc0::fmav_info::shape_t;
+  using stride_t = ducc0::fmav_info::stride_t;
+
+  inline void clear() {}
+
+  inline void fwd(Complex* dst, const Scalar* src, int nfft) {
+    const shape_t axes{0};
+    ducc0::cfmav<Scalar> m_in(src, shape_t{static_cast<size_t>(nfft)});
+    ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft) / 2 + 1});
+    ducc0::r2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+  }
+
+  inline void fwd(Complex* dst, const Complex* src, int nfft) {
+    const shape_t axes{0};
+    ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft)});
+    ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+    ducc0::c2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+  }
+
+  inline void inv(Scalar* dst, const Complex* src, int nfft) {
+    const shape_t axes{0};
+    ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft) / 2 + 1});
+    ducc0::vfmav<Scalar> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+    ducc0::c2r(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+  }
+
+  inline void inv(Complex* dst, const Complex* src, int nfft) {
+    const shape_t axes{0};
+    ducc0::cfmav<Complex> m_in(src, shape_t{static_cast<size_t>(nfft)});
+    ducc0::vfmav<Complex> m_out(dst, shape_t{static_cast<size_t>(nfft)});
+    ducc0::c2c(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+  }
+
+  inline void fwd2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
+    const shape_t axes{0, 1};
+    const shape_t in_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+    const shape_t out_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+    const stride_t stride{static_cast<ptrdiff_t>(nfft1), static_cast<ptrdiff_t>(1)};
+    ducc0::cfmav<Complex> m_in(src, in_shape, stride);
+    ducc0::vfmav<Complex> m_out(dst, out_shape, stride);
+    ducc0::c2c(m_in, m_out, axes, /*forward=*/true, /*scale=*/static_cast<Scalar>(1));
+  }
+
+  inline void inv2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
+    const shape_t axes{0, 1};
+    const shape_t in_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+    const shape_t out_shape{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
+    const stride_t stride{static_cast<ptrdiff_t>(nfft1), static_cast<ptrdiff_t>(1)};
+    ducc0::cfmav<Complex> m_in(src, in_shape, stride);
+    ducc0::vfmav<Complex> m_out(dst, out_shape, stride);
+    ducc0::c2c(m_in, m_out, axes, /*forward=*/false, /*scale=*/static_cast<Scalar>(1));
+  }
+};
+
+}  // namespace internal
+}  // namespace Eigen

unsupported/Eigen/src/FFT/pocketfft_impl.h

@@ -5,17 +5,16 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-using namespace pocketfft;
-using namespace pocketfft::detail;
-
 namespace Eigen {
 
 namespace internal {
 
 template <typename _Scalar>
 struct pocketfft_impl {
-  typedef _Scalar Scalar;
+  using Scalar = _Scalar;
-  typedef std::complex<Scalar> Complex;
+  using Complex = std::complex<Scalar>;
+  using shape_t = pocketfft::shape_t;
+  using stride_t = pocketfft::stride_t;
 
   inline void clear() {}
 
@@ -24,14 +23,14 @@ struct pocketfft_impl {
     const shape_t axes_{0};
     const stride_t stride_in{sizeof(Scalar)};
     const stride_t stride_out{sizeof(Complex)};
-    r2c(shape_, stride_in, stride_out, axes_, FORWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::r2c(shape_, stride_in, stride_out, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
   }
 
   inline void fwd(Complex* dst, const Complex* src, int nfft) {
     const shape_t shape_{static_cast<size_t>(nfft)};
     const shape_t axes_{0};
     const stride_t stride_{sizeof(Complex)};
-    c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
   }
 
   inline void inv(Scalar* dst, const Complex* src, int nfft) {
@@ -39,28 +38,28 @@ struct pocketfft_impl {
     const shape_t axes_{0};
     const stride_t stride_in{sizeof(Complex)};
     const stride_t stride_out{sizeof(Scalar)};
-    c2r(shape_, stride_in, stride_out, axes_, BACKWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::c2r(shape_, stride_in, stride_out, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
   }
 
   inline void inv(Complex* dst, const Complex* src, int nfft) {
     const shape_t shape_{static_cast<size_t>(nfft)};
     const shape_t axes_{0};
     const stride_t stride_{sizeof(Complex)};
-    c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
   }
 
   inline void fwd2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
     const shape_t shape_{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
     const shape_t axes_{0, 1};
     const stride_t stride_{static_cast<ptrdiff_t>(sizeof(Complex) * nfft1), static_cast<ptrdiff_t>(sizeof(Complex))};
-    c2c(shape_, stride_, stride_, axes_, FORWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::FORWARD, src, dst, static_cast<Scalar>(1));
   }
 
   inline void inv2(Complex* dst, const Complex* src, int nfft0, int nfft1) {
     const shape_t shape_{static_cast<size_t>(nfft0), static_cast<size_t>(nfft1)};
     const shape_t axes_{0, 1};
     const stride_t stride_{static_cast<ptrdiff_t>(sizeof(Complex) * nfft1), static_cast<ptrdiff_t>(sizeof(Complex))};
-    c2c(shape_, stride_, stride_, axes_, BACKWARD, src, dst, static_cast<Scalar>(1));
+    pocketfft::c2c(shape_, stride_, stride_, axes_, pocketfft::BACKWARD, src, dst, static_cast<Scalar>(1));
   }
 };
 

unsupported/test/CMakeLists.txt

@@ -88,6 +88,25 @@ else()
   ei_add_property(EIGEN_MISSING_BACKENDS "pocketfft, ")
 endif()
 
+if( NOT DUCC_ROOT AND ENV{DUCC_ROOT} )
+  set( DUCC_ROOT $ENV{DUCC_ROOT} )
+endif()
+find_path(DUCCFFT
+  NAMES "src/ducc0/fft/fft.h"
+  PATHS ${DUCC_ROOT})
+message(INFO " ${DUCC_ROOT} ${DUCCFFT}")
+if(DUCCFFT)
+  ei_add_property(EIGEN_TESTED_BACKENDS "duccfft, ")
+  include_directories( "${DUCCFFT}/src" )
+  add_library(ducc_lib "${DUCCFFT}/src/ducc0/infra/string_utils.cc" "${DUCCFFT}/src/ducc0/infra/threading.cc")
+  target_compile_definitions(ducc_lib PUBLIC "DUCC0_NO_THREADING=1")
+  ei_add_test(duccfft "-DEIGEN_DUCCFFT_DEFAULT -DDUCC0_NO_THREADING=1" "ducc_lib" )
+  set_target_properties(ducc_lib duccfft PROPERTIES CXX_STANDARD 17)
+else()
+  ei_add_property(EIGEN_MISSING_BACKENDS "duccfft, ")
+endif()
+
+
 option(EIGEN_TEST_OPENGL "Enable OpenGL support in unit tests" OFF)
 if(EIGEN_TEST_OPENGL)
   find_package(OpenGL)

unsupported/test/duccfft.cpp

@@ -0,0 +1,4 @@
+#define EIGEN_DUCCFFT_DEFAULT 1
+#include <ducc0/fft/fft.h>         // Needs to be included before main.h
+#include <ducc0/fft/fftnd_impl.h>  // Same requirement
+#include "fft_test_shared.h"

unsupported/test/fft_test_shared.h

@@ -272,7 +272,7 @@ EIGEN_DECLARE_TEST(FFTW) {
   CALL_SUBTEST(test_scalar<float>(2 * 3 * 4 * 5 * 7));
   CALL_SUBTEST(test_scalar<double>(2 * 3 * 4 * 5 * 7));
 
-#if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT
+#if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT
   CALL_SUBTEST(test_complex<long double>(32));
   CALL_SUBTEST(test_complex<long double>(256));
   CALL_SUBTEST(test_complex<long double>(3 * 8));
@@ -294,13 +294,15 @@ EIGEN_DECLARE_TEST(FFTW) {
 // fail to build since Eigen limit the stack allocation size,too big here.
 // CALL_SUBTEST( ( test_complex2d<long double, 256, 256> () ) );
 #endif
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+    defined EIGEN_MKL_DEFAULT
   CALL_SUBTEST((test_complex2d<float, 24, 24>()));
   CALL_SUBTEST((test_complex2d<float, 60, 60>()));
   CALL_SUBTEST((test_complex2d<float, 24, 60>()));
   CALL_SUBTEST((test_complex2d<float, 60, 24>()));
 #endif
-#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
+#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_DUCCFFT_DEFAULT || \
+    defined EIGEN_MKL_DEFAULT
   CALL_SUBTEST((test_complex2d<double, 24, 24>()));
   CALL_SUBTEST((test_complex2d<double, 60, 60>()));
   CALL_SUBTEST((test_complex2d<double, 24, 60>()));