tweaks

bench.py

@@ -2,7 +2,6 @@ import numpy as np
 import pypocketfft
 from time import time
 import matplotlib.pyplot as plt
-import math
 
 
 def _l2error(a, b):
@@ -43,9 +42,10 @@ def measure_fftw_np_interface(a, nrepeat, nthr):
 def measure_pypocketfft(a, nrepeat, nthr):
     import pypocketfft as ppf
     tmin = 1e38
+    b = a.copy()
     for i in range(nrepeat):
         t0 = time()
-        b = ppf.c2c(a, forward=True, nthreads=nthr)
+        b = ppf.c2c(a, out=b, forward=True, nthreads=nthr)
         t1 = time()
         tmin = min(tmin, t1-t0)
     return tmin, b

bench_r2c.py

@@ -2,7 +2,13 @@ import numpy as np
 import pypocketfft
 from time import time
 import matplotlib.pyplot as plt
-import math
+
+
+def get_complex_array(real_array, allocfunc):
+    tocomplex = { np.float32: np.complex64, np.float64: np.complex128 }
+    shape = list(real_array.shape)
+    shape[-1] = shape[-1]//2 + 1
+    return allocfunc(shape, dtype=tocomplex[real_array.dtype.type])
 
 
 def _l2error(a, b):
@@ -12,11 +18,7 @@ def _l2error(a, b):
 def measure_fftw(a, nrepeat, nthr, flags=('FFTW_MEASURE',)):
     import pyfftw
     f1 = pyfftw.empty_aligned(a.shape, dtype=a.dtype)
-    tval = np.ones(1).astype(f1.dtype)
-    t2 = (tval+1j*tval).dtype
-    shape_out=list(a.shape)
-    shape_out[-1] = shape_out[-1]//2 + 1
-    f2 = pyfftw.empty_aligned(shape_out, dtype=t2)
+    f2 = get_complex_array(a, pyfftw.empty_aligned)
     fftw = pyfftw.FFTW(f1, f2, flags=flags, axes=range(a.ndim), threads=nthr)
     f1[()] = a
     tmin = 1e38
@@ -47,9 +49,10 @@ def measure_fftw_np_interface(a, nrepeat, nthr):
 def measure_pypocketfft(a, nrepeat, nthr):
     import pypocketfft as ppf
     tmin = 1e38
+    b = get_complex_array(a, np.empty)
     for i in range(nrepeat):
         t0 = time()
-        b = ppf.r2c(a, forward=True, nthreads=nthr)
+        b = ppf.r2c(a, forward=True, nthreads=nthr, out=b)
         t1 = time()
         tmin = min(tmin, t1-t0)
     return tmin, b

pocketfft_hdronly.h

@@ -3392,6 +3392,9 @@ template<typename T> void r2r_genuine_hartley(const shape_t &shape,
   const T *data_in, T *data_out, T fct, size_t nthreads=1)
   {
   if (util::prod(shape)==0) return;
+  if (axes.size()==1)
+    return r2r_separable_hartley(shape, stride_in, stride_out, axes, data_in,
+      data_out, fct, nthreads);
   util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
   shape_t tshp(shape);
   tshp[axes.back()] = tshp[axes.back()]/2+1;