From cb72d6026521c760b1344f948e262252d85afc14 Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sat, 3 Aug 2019 17:17:44 +0100
Subject: [PATCH 1/7] Use alias templates for vtype

---
 pocketfft_hdronly.h | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index e6bfd73..b1d5a4d 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2812,6 +2812,8 @@ template<> struct VTYPE<long double>
   {
   using type = long double __attribute__ ((vector_size (VLEN<long double>::val*sizeof(long double))));
   };
+
+template <typename T> using vtype_t = typename VTYPE<T>::type;
 #endif
 
 template<typename T> arr<char> alloc_tmp(const shape_t &shape,
@@ -2866,9 +2868,8 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
     if (vlen>1)
       while (it.remaining()>=vlen)
         {
-        using vtype = typename VTYPE<T>::type;
         it.advance(vlen);
-        auto tdatav = reinterpret_cast<cmplx<vtype> *>(storage.data());
+        auto tdatav = reinterpret_cast<cmplx<vtype_t<T>> *>(storage.data());
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             {
@@ -2924,9 +2925,8 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
     if (vlen>1)
       while (it.remaining()>=vlen)
         {
-        using vtype = typename VTYPE<T>::type;
         it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype *>(storage.data());
+        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             tdatav[i][j] = tin[it.iofs(j,i)];
@@ -2992,9 +2992,8 @@ template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
     if (vlen>1)
       while (it.remaining()>=vlen)
         {
-        using vtype = typename VTYPE<T>::type;
         it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype *>(storage.data());
+        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             tdatav[i][j] = tin[it.iofs(j,i)];
@@ -3040,9 +3039,8 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
   if (vlen>1)
     while (it.remaining()>=vlen)
       {
-      using vtype = typename VTYPE<T>::type;
       it.advance(vlen);
-      auto tdatav = reinterpret_cast<vtype *>(storage.data());
+      auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
       for (size_t i=0; i<len; ++i)
         for (size_t j=0; j<vlen; ++j)
           tdatav[i][j] = in[it.iofs(j,i)];
@@ -3100,9 +3098,8 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
   if (vlen>1)
     while (it.remaining()>=vlen)
       {
-      using vtype = typename VTYPE<T>::type;
       it.advance(vlen);
-      auto tdatav = reinterpret_cast<vtype *>(storage.data());
+      auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
       for (size_t j=0; j<vlen; ++j)
         tdatav[0][j]=in[it.iofs(j,0)].r;
       {
@@ -3172,9 +3169,8 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
     if (vlen>1)
       while (it.remaining()>=vlen)
         {
-        using vtype = typename VTYPE<T>::type;
         it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype *>(storage.data());
+        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             tdatav[i][j] = tin[it.iofs(j,i)];
-- 
GitLab


From 343ff57eee80bec3978b972dd2720523032b30ff Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sat, 3 Aug 2019 17:57:40 +0100
Subject: [PATCH 2/7] Isolate copying into its own function

---
 pocketfft_hdronly.h | 178 +++++++++++++++++++++++++-------------------
 1 file changed, 102 insertions(+), 76 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index b1d5a4d..305fa4d 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2844,6 +2844,57 @@ template<typename T> arr<char> alloc_tmp(const shape_t &shape,
 #define POCKETFFT_NTHREADS
 #endif
 
+template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
+  const cndarr<cmplx<T>> &src, cmplx<vtype_t<T>> *POCKETFFT_RESTRICT dst)
+  {
+  for (size_t i=0; i<it.length_in(); ++i)
+    for (size_t j=0; j<vlen; ++j)
+      {
+      dst[i].r[j] = src[it.iofs(j,i)].r;
+      dst[i].i[j] = src[it.iofs(j,i)].i;
+      }
+  }
+
+template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
+  const cndarr<T> &src, vtype_t<T> *POCKETFFT_RESTRICT dst)
+  {
+  for (size_t i=0; i<it.length_in(); ++i)
+    for (size_t j=0; j<vlen; ++j)
+      dst[i][j] = src[it.iofs(j,i)];
+  }
+
+template <typename T, size_t vlen> void copy_input(const multi_iter<vlen> &it,
+  const cndarr<T> &src, T *POCKETFFT_RESTRICT dst)
+  {
+  if (dst == &src[it.iofs(0)]) return;  // in-place
+  for (size_t i=0; i<it.length_in(); ++i)
+    dst[i] = src[it.iofs(i)];
+  }
+
+template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
+  const cmplx<vtype_t<T>> *POCKETFFT_RESTRICT src, ndarr<cmplx<T>> &dst)
+  {
+  for (size_t i=0; i<it.length_out(); ++i)
+    for (size_t j=0; j<vlen; ++j)
+      dst[it.oofs(j,i)].Set(src[i].r[j],src[i].i[j]);
+  }
+
+template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
+  const vtype_t<T> *POCKETFFT_RESTRICT src, ndarr<T> &dst)
+  {
+  for (size_t i=0; i<it.length_out(); ++i)
+    for (size_t j=0; j<vlen; ++j)
+      dst[it.oofs(j,i)] = src[i][j];
+  }
+
+template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
+  const T *POCKETFFT_RESTRICT src, ndarr<T> &dst)
+  {
+  if (src == &dst[it.oofs(0)]) return;  // in-place
+  for (size_t i=0; i<it.length_out(); ++i)
+    dst[it.oofs(i)] = src[i];
+  }
+
 template<typename T> POCKETFFT_NOINLINE void general_c(
   const cndarr<cmplx<T>> &in, ndarr<cmplx<T>> &out,
   const shape_t &axes, bool forward, T fct, size_t POCKETFFT_NTHREADS)
@@ -2870,16 +2921,9 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
         {
         it.advance(vlen);
         auto tdatav = reinterpret_cast<cmplx<vtype_t<T>> *>(storage.data());
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            {
-            tdatav[i].r[j] = tin[it.iofs(j,i)].r;
-            tdatav[i].i[j] = tin[it.iofs(j,i)].i;
-            }
+        copy_input(it, tin, tdatav);
         plan->exec (tdatav, fct, forward);
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            out[it.oofs(j,i)].Set(tdatav[i].r[j],tdatav[i].i[j]);
+        copy_output(it, tdatav, out);
         }
 #endif
     while (it.remaining()>0)
@@ -2888,19 +2932,46 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
       auto buf = it.stride_out() == sizeof(cmplx<T>) ?
         &out[it.oofs(0)] : reinterpret_cast<cmplx<T> *>(storage.data());
 
-      if (buf != &tin[it.iofs(0)])
-        for (size_t i=0; i<len; ++i)
-          buf[i] = tin[it.iofs(i)];
+      copy_input(it, tin, buf);
       plan->exec (buf, fct, forward);
-      if (buf != &out[it.oofs(0)])
-        for (size_t i=0; i<len; ++i)
-          out[it.oofs(i)] = buf[i];
+      copy_output(it, buf, out);
       }
 } // end of parallel region
     fct = T(1); // factor has been applied, use 1 for remaining axes
     }
   }
 
+template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
+  const vtype_t<T> *POCKETFFT_RESTRICT src, ndarr<T> &dst)
+  {
+  for (size_t j=0; j<vlen; ++j)
+    dst[it.oofs(j,0)] = src[0][j];
+  size_t i=1, i1=1, i2=it.length_out()-1;
+  for (i=1; i<it.length_out()-1; i+=2, ++i1, --i2)
+    for (size_t j=0; j<vlen; ++j)
+      {
+        dst[it.oofs(j,i1)] = src[i][j]+src[i+1][j];
+        dst[it.oofs(j,i2)] = src[i][j]-src[i+1][j];
+      }
+  if (i<it.length_out())
+    for (size_t j=0; j<vlen; ++j)
+      dst[it.oofs(j,i1)] = src[i][j];
+  }
+
+template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
+  const T *POCKETFFT_RESTRICT src, ndarr<T> &dst)
+  {
+  dst[it.oofs(0)] = src[0];
+  size_t i=1, i1=1, i2=it.length_out()-1;
+  for (i=1; i<it.length_out()-1; i+=2, ++i1, --i2)
+    {
+      dst[it.oofs(i1)] = src[i]+src[i+1];
+      dst[it.oofs(i2)] = src[i]-src[i+1];
+    }
+  if (i<it.length_out())
+    dst[it.oofs(i1)] = src[i];
+  }
+
 template<typename T> POCKETFFT_NOINLINE void general_hartley(
   const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, T fct,
   size_t POCKETFFT_NTHREADS)
@@ -2927,41 +2998,18 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
         {
         it.advance(vlen);
         auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            tdatav[i][j] = tin[it.iofs(j,i)];
+        copy_input(it, tin, tdatav);
         plan->exec(tdatav, fct, true);
-        for (size_t j=0; j<vlen; ++j)
-          out[it.oofs(j,0)] = tdatav[0][j];
-        size_t i=1, i1=1, i2=len-1;
-        for (i=1; i<len-1; i+=2, ++i1, --i2)
-          for (size_t j=0; j<vlen; ++j)
-            {
-            out[it.oofs(j,i1)] = tdatav[i][j]+tdatav[i+1][j];
-            out[it.oofs(j,i2)] = tdatav[i][j]-tdatav[i+1][j];
-            }
-        if (i<len)
-          for (size_t j=0; j<vlen; ++j)
-            out[it.oofs(j,i1)] = tdatav[i][j];
+        copy_hartley(it, tdatav, out);
         }
 #endif
     while (it.remaining()>0)
       {
       it.advance(1);
       auto tdata = reinterpret_cast<T *>(storage.data());
-      for (size_t i=0; i<len; ++i)
-        tdata[i] = tin[it.iofs(i)];
+      copy_input(it, tin, tdata);
       plan->exec(tdata, fct, true);
-      // Hartley order
-      out[it.oofs(0)] = tdata[0];
-      size_t i=1, i1=1, i2=len-1;
-      for (i=1; i<len-1; i+=2, ++i1, --i2)
-        {
-        out[it.oofs(i1)] = tdata[i]+tdata[i+1];
-        out[it.oofs(i2)] = tdata[i]-tdata[i+1];
-        }
-      if (i<len)
-        out[it.oofs(i1)] = tdata[i];
+      copy_hartley(it, tdata, out);
       }
 } // end of parallel region
     fct = T(1); // factor has been applied, use 1 for remaining axes
@@ -2994,13 +3042,9 @@ template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
         {
         it.advance(vlen);
         auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            tdatav[i][j] = tin[it.iofs(j,i)];
+        copy_input(it, tin, tdatav);
         plan->exec(tdatav, fct, ortho, type, cosine);
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            out[it.oofs(j,i)] = tdatav[i][j];
+        copy_output(it, tdatav, out);
         }
 #endif
     while (it.remaining()>0)
@@ -3009,13 +3053,9 @@ template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
       auto buf = it.stride_out() == sizeof(T) ? &out[it.oofs(0)]
         : reinterpret_cast<T *>(storage.data());
 
-      if (buf != &tin[it.iofs(0)])
-        for (size_t i=0; i<len; ++i)
-          buf[i] = tin[it.iofs(i)];
+      copy_input(it, tin, buf);
       plan->exec(buf, fct, ortho, type, cosine);
-      if (buf != &out[it.oofs(0)])
-        for (size_t i=0; i<len; ++i)
-          out[it.oofs(i)] = buf[i];
+      copy_output(it, buf, out);
       }
 } // end of parallel region
     fct = T(1); // factor has been applied, use 1 for remaining axes
@@ -3041,9 +3081,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
       {
       it.advance(vlen);
       auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-      for (size_t i=0; i<len; ++i)
-        for (size_t j=0; j<vlen; ++j)
-          tdatav[i][j] = in[it.iofs(j,i)];
+      copy_input(it, in, tdatav);
       plan->exec(tdatav, fct, true);
       for (size_t j=0; j<vlen; ++j)
         out[it.oofs(j,0)].Set(tdatav[0][j]);
@@ -3065,8 +3103,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
     {
     it.advance(1);
     auto tdata = reinterpret_cast<T *>(storage.data());
-    for (size_t i=0; i<len; ++i)
-      tdata[i] = in[it.iofs(i)];
+    copy_input(it, in, tdata);
     plan->exec(tdata, fct, true);
     out[it.oofs(0)].Set(tdata[0]);
     size_t i=1, ii=1;
@@ -3117,9 +3154,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
           tdatav[i][j] = in[it.iofs(j,ii)].r;
       }
       plan->exec(tdatav, fct, false);
-      for (size_t i=0; i<len; ++i)
-        for (size_t j=0; j<vlen; ++j)
-          out[it.oofs(j,i)] = tdatav[i][j];
+      copy_output(it, tdatav, out);
       }
 #endif
   while (it.remaining()>0)
@@ -3139,8 +3174,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
       tdata[i] = in[it.iofs(ii)].r;
     }
     plan->exec(tdata, fct, false);
-    for (size_t i=0; i<len; ++i)
-      out[it.oofs(i)] = tdata[i];
+    copy_output(it, tdata, out);
     }
 } // end of parallel region
   }
@@ -3171,9 +3205,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
         {
         it.advance(vlen);
         auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            tdatav[i][j] = tin[it.iofs(j,i)];
+        copy_input(it, tin, tdatav);
         if ((!r2c) && forward)
           for (size_t i=2; i<len; i+=2)
             tdatav[i] = -tdatav[i];
@@ -3181,9 +3213,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
         if (r2c && (!forward))
           for (size_t i=2; i<len; i+=2)
             tdatav[i] = -tdatav[i];
-        for (size_t i=0; i<len; ++i)
-          for (size_t j=0; j<vlen; ++j)
-            out[it.oofs(j,i)] = tdatav[i][j];
+        copy_output(it, tdatav, out);
         }
 #endif
     while (it.remaining()>0)
@@ -3192,9 +3222,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
       auto buf = it.stride_out() == sizeof(T) ?
         &out[it.oofs(0)] : reinterpret_cast<T *>(storage.data());
 
-      if (buf != &tin[it.iofs(0)])
-        for (size_t i=0; i<len; ++i)
-          buf[i] = tin[it.iofs(i)];
+      copy_input(it, tin, buf);
       if ((!r2c) && forward)
         for (size_t i=2; i<len; i+=2)
           buf[i] = -buf[i];
@@ -3202,9 +3230,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r(
       if (r2c && (!forward))
         for (size_t i=2; i<len; i+=2)
           buf[i] = -buf[i];
-      if (buf != &out[it.oofs(0)])
-        for (size_t i=0; i<len; ++i)
-          out[it.oofs(i)] = buf[i];
+      copy_output(it, buf, out);
       }
 } // end of parallel region
     fct = T(1); // factor has been applied, use 1 for remaining axes
-- 
GitLab


From 4ef6e3bfb03adc4e0736ebce522bfbd8ca3c315f Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sat, 3 Aug 2019 20:29:47 +0100
Subject: [PATCH 3/7] Refactor n-dimensional general_x functions

---
 pocketfft_hdronly.h | 239 +++++++++++++++++---------------------------
 1 file changed, 91 insertions(+), 148 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 305fa4d..895ee0d 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2895,24 +2895,29 @@ template<typename T, size_t vlen> void copy_output(const multi_iter<vlen> &it,
     dst[it.oofs(i)] = src[i];
   }
 
-template<typename T> POCKETFFT_NOINLINE void general_c(
-  const cndarr<cmplx<T>> &in, ndarr<cmplx<T>> &out,
-  const shape_t &axes, bool forward, T fct, size_t POCKETFFT_NTHREADS)
+template <typename T> struct add_vec { using type = vtype_t<T>; };
+template <typename T> struct add_vec<cmplx<T>>
+  { using type = cmplx<vtype_t<T>>; };
+template <typename T> using add_vec_t = typename add_vec<T>::type;
+
+template<typename Tplan, typename T, typename T0, typename Exec>
+POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
+  const shape_t &axes, T0 fct, size_t POCKETFFT_NTHREADS, const Exec & exec)
   {
-  shared_ptr<pocketfft_c<T>> plan;
+  shared_ptr<Tplan> plan;
 
   for (size_t iax=0; iax<axes.size(); ++iax)
     {
-    constexpr auto vlen = VLEN<T>::val;
+    constexpr auto vlen = VLEN<T0>::val;
     size_t len=in.shape(axes[iax]);
     if ((!plan) || (len!=plan->length()))
-      plan = get_plan<pocketfft_c<T>>(len);
+      plan = get_plan<Tplan>(len);
 
 #ifdef POCKETFFT_OPENMP
 #pragma omp parallel num_threads(util::thread_count(nthreads, in.shape(), axes[iax]))
 #endif
 {
-    auto storage = alloc_tmp<T>(in.shape(), len, sizeof(cmplx<T>));
+    auto storage = alloc_tmp<T0>(in.shape(), len, sizeof(T));
     const auto &tin(iax==0? in : out);
     multi_iter<vlen> it(tin, out, axes[iax]);
 #ifndef POCKETFFT_NO_VECTORS
@@ -2920,27 +2925,43 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
       while (it.remaining()>=vlen)
         {
         it.advance(vlen);
-        auto tdatav = reinterpret_cast<cmplx<vtype_t<T>> *>(storage.data());
-        copy_input(it, tin, tdatav);
-        plan->exec (tdatav, fct, forward);
-        copy_output(it, tdatav, out);
+        auto tdatav = reinterpret_cast<add_vec_t<T> *>(storage.data());
+        exec(it, in, out, tdatav, *plan, fct);
         }
 #endif
     while (it.remaining()>0)
       {
       it.advance(1);
-      auto buf = it.stride_out() == sizeof(cmplx<T>) ?
-        &out[it.oofs(0)] : reinterpret_cast<cmplx<T> *>(storage.data());
-
-      copy_input(it, tin, buf);
-      plan->exec (buf, fct, forward);
-      copy_output(it, buf, out);
+      auto buf = it.stride_out() == sizeof(T) ?
+        &out[it.oofs(0)] : reinterpret_cast<T *>(storage.data());
+      exec(it, in, out, buf, *plan, fct);
       }
 } // end of parallel region
-    fct = T(1); // factor has been applied, use 1 for remaining axes
+    fct = T0(1); // factor has been applied, use 1 for remaining axes
     }
   }
 
+struct ExecC2C
+  {
+  bool forward;
+
+  template <typename T0, typename T, size_t vlen> void operator () (
+    const multi_iter<vlen> &it, const cndarr<cmplx<T0>> &in,
+    ndarr<cmplx<T0>> &out, T * buf, const pocketfft_c<T0> &plan, T0 fct) const
+    {
+    copy_input(it, in, buf);
+    plan.exec(buf, fct, forward);
+    copy_output(it, buf, out);
+    }
+  };
+
+template<typename T> void general_c(
+  const cndarr<cmplx<T>> &in, ndarr<cmplx<T>> &out,
+  const shape_t &axes, bool forward, T fct, size_t nthreads)
+  {
+  general_nd<pocketfft_c<T>>(in, out, axes, fct, nthreads, ExecC2C{forward});
+  }
+
 template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
   const vtype_t<T> *POCKETFFT_RESTRICT src, ndarr<T> &dst)
   {
@@ -2972,94 +2993,46 @@ template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
     dst[it.oofs(i1)] = src[i];
   }
 
+struct ExecHartley
+  {
+  template <typename T0, typename T, size_t vlen> void operator () (
+    const multi_iter<vlen> &it, const cndarr<T0> &in, ndarr<T0> &out,
+    T * buf, const pocketfft_r<T0> &plan, T0 fct) const
+    {
+    copy_input(it, in, buf);
+    plan.exec(buf, fct, true);
+    copy_hartley(it, buf, out);
+    }
+  };
+
 template<typename T> POCKETFFT_NOINLINE void general_hartley(
   const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, T fct,
-  size_t POCKETFFT_NTHREADS)
+  size_t nthreads)
   {
-  shared_ptr<pocketfft_r<T>> plan;
+  general_nd<pocketfft_r<T>>(in, out, axes, fct, nthreads, ExecHartley{});
+  }
 
-  for (size_t iax=0; iax<axes.size(); ++iax)
-    {
-    constexpr auto vlen = VLEN<T>::val;
-    size_t len=in.shape(axes[iax]);
-    if ((!plan) || (len!=plan->length()))
-      plan = get_plan<pocketfft_r<T>>(len);
+struct ExecDcst
+  {
+  bool ortho;
+  int type;
+  bool cosine;
 
-#ifdef POCKETFFT_OPENMP
-#pragma omp parallel num_threads(util::thread_count(nthreads, in.shape(), axes[iax]))
-#endif
-{
-    auto storage = alloc_tmp<T>(in.shape(), len, sizeof(T));
-    const auto &tin(iax==0 ? in : out);
-    multi_iter<vlen> it(tin, out, axes[iax]);
-#ifndef POCKETFFT_NO_VECTORS
-    if (vlen>1)
-      while (it.remaining()>=vlen)
-        {
-        it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        copy_input(it, tin, tdatav);
-        plan->exec(tdatav, fct, true);
-        copy_hartley(it, tdatav, out);
-        }
-#endif
-    while (it.remaining()>0)
-      {
-      it.advance(1);
-      auto tdata = reinterpret_cast<T *>(storage.data());
-      copy_input(it, tin, tdata);
-      plan->exec(tdata, fct, true);
-      copy_hartley(it, tdata, out);
-      }
-} // end of parallel region
-    fct = T(1); // factor has been applied, use 1 for remaining axes
+  template <typename T0, typename T, typename Tplan, size_t vlen>
+  void operator () (const multi_iter<vlen> &it, const cndarr<T0> &in,
+    ndarr<T0> &out, T * buf, const Tplan &plan, T0 fct) const
+    {
+    copy_input(it, in, buf);
+    plan.exec(buf, fct, ortho, type, cosine);
+    copy_output(it, buf, out);
     }
-  }
+  };
 
 template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
   const cndarr<T> &in, ndarr<T> &out, const shape_t &axes,
-  T fct, bool ortho, int type, bool cosine, size_t POCKETFFT_NTHREADS)
+  T fct, bool ortho, int type, bool cosine, size_t nthreads)
   {
-  shared_ptr<Trafo> plan;
-
-  for (size_t iax=0; iax<axes.size(); ++iax)
-    {
-    constexpr auto vlen = VLEN<T>::val;
-    size_t len=in.shape(axes[iax]);
-    if ((!plan) || (len!=plan->length()))
-      plan = get_plan<Trafo>(len);
-
-#ifdef POCKETFFT_OPENMP
-#pragma omp parallel num_threads(util::thread_count(nthreads, in.shape(), axes[iax]))
-#endif
-{
-    auto storage = alloc_tmp<T>(in.shape(), len, sizeof(T));
-    const auto &tin(iax==0 ? in : out);
-    multi_iter<vlen> it(tin, out, axes[iax]);
-#ifndef POCKETFFT_NO_VECTORS
-    if (vlen>1)
-      while (it.remaining()>=vlen)
-        {
-        it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        copy_input(it, tin, tdatav);
-        plan->exec(tdatav, fct, ortho, type, cosine);
-        copy_output(it, tdatav, out);
-        }
-#endif
-    while (it.remaining()>0)
-      {
-      it.advance(1);
-      auto buf = it.stride_out() == sizeof(T) ? &out[it.oofs(0)]
-        : reinterpret_cast<T *>(storage.data());
-
-      copy_input(it, tin, buf);
-      plan->exec(buf, fct, ortho, type, cosine);
-      copy_output(it, buf, out);
-      }
-} // end of parallel region
-    fct = T(1); // factor has been applied, use 1 for remaining axes
-    }
+  general_nd<Trafo>(in, out, axes, fct, nthreads, ExecDcst{ortho, type, cosine});
   }
 
 template<typename T> POCKETFFT_NOINLINE void general_r2c(
@@ -3179,62 +3152,32 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
 } // end of parallel region
   }
 
-template<typename T> POCKETFFT_NOINLINE void general_r(
-  const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, bool r2c,
-  bool forward, T fct, size_t POCKETFFT_NTHREADS)
+struct ExecR2R
   {
-  shared_ptr<pocketfft_r<T>> plan;
+  bool r2c, forward;
 
-  for (size_t iax=0; iax<axes.size(); ++iax)
+  template <typename T0, typename T, size_t vlen> void operator () (
+    const multi_iter<vlen> &it, const cndarr<T0> &in, ndarr<T0> &out, T * buf,
+    const pocketfft_r<T0> &plan, T0 fct) const
     {
-    constexpr auto vlen = VLEN<T>::val;
-    size_t len=in.shape(axes[iax]);
-    if ((!plan) || (len!=plan->length()))
-      plan = get_plan<pocketfft_r<T>>(len);
-
-#ifdef POCKETFFT_OPENMP
-#pragma omp parallel num_threads(util::thread_count(nthreads, in.shape(), axes[iax]))
-#endif
-{
-    auto storage = alloc_tmp<T>(in.shape(), len, sizeof(T));
-    const auto &tin(iax==0 ? in : out);
-    multi_iter<vlen> it(tin, out, axes[iax]);
-#ifndef POCKETFFT_NO_VECTORS
-    if (vlen>1)
-      while (it.remaining()>=vlen)
-        {
-        it.advance(vlen);
-        auto tdatav = reinterpret_cast<vtype_t<T> *>(storage.data());
-        copy_input(it, tin, tdatav);
-        if ((!r2c) && forward)
-          for (size_t i=2; i<len; i+=2)
-            tdatav[i] = -tdatav[i];
-        plan->exec (tdatav, fct, forward);
-        if (r2c && (!forward))
-          for (size_t i=2; i<len; i+=2)
-            tdatav[i] = -tdatav[i];
-        copy_output(it, tdatav, out);
-        }
-#endif
-    while (it.remaining()>0)
-      {
-      it.advance(1);
-      auto buf = it.stride_out() == sizeof(T) ?
-        &out[it.oofs(0)] : reinterpret_cast<T *>(storage.data());
-
-      copy_input(it, tin, buf);
-      if ((!r2c) && forward)
-        for (size_t i=2; i<len; i+=2)
-          buf[i] = -buf[i];
-      plan->exec(buf, fct, forward);
-      if (r2c && (!forward))
-        for (size_t i=2; i<len; i+=2)
-          buf[i] = -buf[i];
-      copy_output(it, buf, out);
-      }
-} // end of parallel region
-    fct = T(1); // factor has been applied, use 1 for remaining axes
+    copy_input(it, in, buf);
+    if ((!r2c) && forward)
+      for (size_t i=2; i<it.length_out(); i+=2)
+        buf[i] = -buf[i];
+    plan.exec(buf, fct, forward);
+    if (r2c && (!forward))
+      for (size_t i=2; i<it.length_out(); i+=2)
+        buf[i] = -buf[i];
+    copy_output(it, buf, out);
     }
+  };
+
+template<typename T> POCKETFFT_NOINLINE void general_r(
+  const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, bool r2c,
+  bool forward, T fct, size_t nthreads)
+  {
+  general_nd<pocketfft_r<T>>(in, out, axes, fct, nthreads,
+    ExecR2R{r2c, forward});
   }
 
 #undef POCKETFFT_NTHREADS
-- 
GitLab


From 1971f449e91330b4e92565c527caae832ad8815d Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sun, 4 Aug 2019 01:10:08 +0100
Subject: [PATCH 4/7] Remove specific general_ functions {c, r, hartley & dcst}

---
 pocketfft_hdronly.h | 50 +++++++++++----------------------------------
 1 file changed, 12 insertions(+), 38 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 895ee0d..7afc983 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2955,13 +2955,6 @@ struct ExecC2C
     }
   };
 
-template<typename T> void general_c(
-  const cndarr<cmplx<T>> &in, ndarr<cmplx<T>> &out,
-  const shape_t &axes, bool forward, T fct, size_t nthreads)
-  {
-  general_nd<pocketfft_c<T>>(in, out, axes, fct, nthreads, ExecC2C{forward});
-  }
-
 template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
   const vtype_t<T> *POCKETFFT_RESTRICT src, ndarr<T> &dst)
   {
@@ -3005,13 +2998,6 @@ struct ExecHartley
     }
   };
 
-template<typename T> POCKETFFT_NOINLINE void general_hartley(
-  const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, T fct,
-  size_t nthreads)
-  {
-  general_nd<pocketfft_r<T>>(in, out, axes, fct, nthreads, ExecHartley{});
-  }
-
 struct ExecDcst
   {
   bool ortho;
@@ -3028,13 +3014,6 @@ struct ExecDcst
     }
   };
 
-template<typename Trafo, typename T> POCKETFFT_NOINLINE void general_dcst(
-  const cndarr<T> &in, ndarr<T> &out, const shape_t &axes,
-  T fct, bool ortho, int type, bool cosine, size_t nthreads)
-  {
-  general_nd<Trafo>(in, out, axes, fct, nthreads, ExecDcst{ortho, type, cosine});
-  }
-
 template<typename T> POCKETFFT_NOINLINE void general_r2c(
   const cndarr<T> &in, ndarr<cmplx<T>> &out, size_t axis, bool forward, T fct,
   size_t POCKETFFT_NTHREADS)
@@ -3172,14 +3151,6 @@ struct ExecR2R
     }
   };
 
-template<typename T> POCKETFFT_NOINLINE void general_r(
-  const cndarr<T> &in, ndarr<T> &out, const shape_t &axes, bool r2c,
-  bool forward, T fct, size_t nthreads)
-  {
-  general_nd<pocketfft_r<T>>(in, out, axes, fct, nthreads,
-    ExecR2R{r2c, forward});
-  }
-
 #undef POCKETFFT_NTHREADS
 
 template<typename T> void c2c(const shape_t &shape, const stride_t &stride_in,
@@ -3191,7 +3162,7 @@ template<typename T> void c2c(const shape_t &shape, const stride_t &stride_in,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<cmplx<T>> ain(data_in, shape, stride_in);
   ndarr<cmplx<T>> aout(data_out, shape, stride_out);
-  general_c(ain, aout, axes, forward, fct, nthreads);
+  general_nd<pocketfft_c<T>>(ain, aout, axes, fct, nthreads, ExecC2C{forward});
   }
 
 template<typename T> void dct(const shape_t &shape,
@@ -3203,12 +3174,13 @@ template<typename T> void dct(const shape_t &shape,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
+  const ExecDcst exec{ortho, type, true};
   if (type==1)
-    general_dcst<T_dct1<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
+    general_nd<T_dct1<T>>(ain, aout, axes, fct, nthreads, exec);
   else if (type==4)
-    general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
+    general_nd<T_dcst4<T>>(ain, aout, axes, fct, nthreads, exec);
   else
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, true, nthreads);
+    general_nd<T_dcst23<T>>(ain, aout, axes, fct, nthreads, exec);
   }
 
 template<typename T> void dst(const shape_t &shape,
@@ -3220,12 +3192,13 @@ template<typename T> void dst(const shape_t &shape,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
+  const ExecDcst exec{ortho, type, false};
   if (type==1)
-    general_dcst<T_dst1<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
+    general_nd<T_dst1<T>>(ain, aout, axes, fct, nthreads, exec);
   else if (type==4)
-    general_dcst<T_dcst4<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
+    general_nd<T_dcst4<T>>(ain, aout, axes, fct, nthreads, exec);
   else
-    general_dcst<T_dcst23<T>>(ain, aout, axes, fct, ortho, type, false, nthreads);
+    general_nd<T_dcst23<T>>(ain, aout, axes, fct, nthreads, exec);
   }
 
 template<typename T> void r2c(const shape_t &shape_in,
@@ -3309,7 +3282,8 @@ template<typename T> void r2r_fftpack(const shape_t &shape,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
-  general_r(ain, aout, axes, real2hermitian, forward, fct, nthreads);
+  general_nd<pocketfft_r<T>>(ain, aout, axes, fct, nthreads,
+    ExecR2R{real2hermitian, forward});
   }
 
 template<typename T> void r2r_separable_hartley(const shape_t &shape,
@@ -3320,7 +3294,7 @@ template<typename T> void r2r_separable_hartley(const shape_t &shape,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
-  general_hartley(ain, aout, axes, fct, nthreads);
+  general_nd<pocketfft_r<T>>(ain, aout, axes, fct, nthreads, ExecHartley{});
   }
 
 } // namespace detail
-- 
GitLab


From ac891a3056b935185205626a6da9998e51dae8fb Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sun, 4 Aug 2019 02:34:43 +0100
Subject: [PATCH 5/7] Fix hartley transforms by disallowing inplace FFT output

---
 pocketfft_hdronly.h | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 7afc983..0b955ba 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2902,7 +2902,8 @@ template <typename T> using add_vec_t = typename add_vec<T>::type;
 
 template<typename Tplan, typename T, typename T0, typename Exec>
 POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
-  const shape_t &axes, T0 fct, size_t POCKETFFT_NTHREADS, const Exec & exec)
+  const shape_t &axes, T0 fct, size_t POCKETFFT_NTHREADS, const Exec & exec,
+  const bool allow_inplace=true)
   {
   shared_ptr<Tplan> plan;
 
@@ -2932,7 +2933,7 @@ POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
     while (it.remaining()>0)
       {
       it.advance(1);
-      auto buf = it.stride_out() == sizeof(T) ?
+      auto buf = allow_inplace && it.stride_out() == sizeof(T) ?
         &out[it.oofs(0)] : reinterpret_cast<T *>(storage.data());
       exec(it, in, out, buf, *plan, fct);
       }
@@ -2979,8 +2980,8 @@ template <typename T, size_t vlen> void copy_hartley(const multi_iter<vlen> &it,
   size_t i=1, i1=1, i2=it.length_out()-1;
   for (i=1; i<it.length_out()-1; i+=2, ++i1, --i2)
     {
-      dst[it.oofs(i1)] = src[i]+src[i+1];
-      dst[it.oofs(i2)] = src[i]-src[i+1];
+    dst[it.oofs(i1)] = src[i]+src[i+1];
+    dst[it.oofs(i2)] = src[i]-src[i+1];
     }
   if (i<it.length_out())
     dst[it.oofs(i1)] = src[i];
@@ -3294,7 +3295,8 @@ template<typename T> void r2r_separable_hartley(const shape_t &shape,
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   cndarr<T> ain(data_in, shape, stride_in);
   ndarr<T> aout(data_out, shape, stride_out);
-  general_nd<pocketfft_r<T>>(ain, aout, axes, fct, nthreads, ExecHartley{});
+  general_nd<pocketfft_r<T>>(ain, aout, axes, fct, nthreads, ExecHartley{},
+    false);
   }
 
 } // namespace detail
-- 
GitLab


From 15278280a70f11f1a0ecb0f17851f5dd13d3614c Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Sun, 4 Aug 2019 15:29:35 +0200
Subject: [PATCH 6/7] fix

---
 pocketfft_hdronly.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 0b955ba..d1eec5d 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -2927,7 +2927,7 @@ POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
         {
         it.advance(vlen);
         auto tdatav = reinterpret_cast<add_vec_t<T> *>(storage.data());
-        exec(it, in, out, tdatav, *plan, fct);
+        exec(it, tin, out, tdatav, *plan, fct);
         }
 #endif
     while (it.remaining()>0)
@@ -2935,7 +2935,7 @@ POCKETFFT_NOINLINE void general_nd(const cndarr<T> &in, ndarr<T> &out,
       it.advance(1);
       auto buf = allow_inplace && it.stride_out() == sizeof(T) ?
         &out[it.oofs(0)] : reinterpret_cast<T *>(storage.data());
-      exec(it, in, out, buf, *plan, fct);
+      exec(it, tin, out, buf, *plan, fct);
       }
 } // end of parallel region
     fct = T0(1); // factor has been applied, use 1 for remaining axes
-- 
GitLab


From 8aa2d5cd467187ec895a45498164aa27608dde41 Mon Sep 17 00:00:00 2001
From: Peter Bell <peterbell10@live.co.uk>
Date: Sun, 4 Aug 2019 16:57:20 +0100
Subject: [PATCH 7/7] xfail long double tests on WSL and platforms without
 extended precision

---
 test.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/test.py b/test.py
index c718c86..e33837b 100644
--- a/test.py
+++ b/test.py
@@ -58,9 +58,18 @@ def irfft_scipy(a, axis, inorm=0, out=None, nthreads=1):
 tol = {np.float32: 6e-7, np.float64: 1.5e-15, np.longfloat: 1e-18}
 ctype = {np.float32: np.complex64, np.float64: np.complex128, np.longfloat: np.longcomplex}
 
+import platform
+on_wsl = "microsoft" in platform.uname()[3].lower()
+true_long_double = (np.longfloat != np.float64 and not on_wsl)
+dtypes = [np.float32, np.float64,
+          pytest.param(np.longfloat, marks=pytest.mark.xfail(
+              not true_long_double,
+              reason="Long double doesn't offer more precision"))]
+
+
 @pmp("len", len1D)
 @pmp("inorm", [0, 1, 2])
-@pmp("dtype", [np.float32, np.float64, np.longfloat])
+@pmp("dtype", dtypes)
 def test1D(len, inorm, dtype):
     a = np.random.rand(len)-0.5 + 1j*np.random.rand(len)-0.5j
     a = a.astype(ctype[dtype])
@@ -199,7 +208,7 @@ def test_genuine_hartley_2D(shp, axes):
 @pmp("len", len1D)
 @pmp("inorm", [0, 1])  # inorm==2 not needed, tested via inverse
 @pmp("type", [1, 2, 3, 4])
-@pmp("dtype", [np.float32, np.float64, np.longfloat])
+@pmp("dtype", dtypes)
 def testdcst1D(len, inorm, type, dtype):
     a = (np.random.rand(len)-0.5).astype(dtype)
     eps = tol[dtype]
-- 
GitLab