From fb991029fd64149ee5af64e659a042764a40b6eb Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 2 Aug 2019 16:19:05 +0200
Subject: [PATCH 1/9] code cosmetics

---
 pocketfft_hdronly.h | 200 +++++++++++++++++++++-----------------------
 1 file changed, 94 insertions(+), 106 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 403d8c9..480087e 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -235,14 +235,19 @@ template<typename T> struct cmplx {
                : Tres(r*other.r-i*other.i, r*other.i+i*other.r);
     }
 };
+template<typename T> inline void PM(T &a, T &b, T c, T d)
+  { a=c+d; b=c-d; }
 template<typename T> inline void PMINPLACE(T &a, T &b)
   { T t = a; a+=b; b=t-b; }
 template<typename T> inline void MPINPLACE(T &a, T &b)
   { T t = a; a-=b; b=t+b; }
-template<typename T> void PMC(T &a, T &b, const T &c, const T &d)
-  { a = c+d; b = c-d; }
 template<typename T> cmplx<T> conj(const cmplx<T> &a)
   { return {a.r, -a.i}; }
+template<bool fwd, typename T, typename T2> void special_mul (const cmplx<T> &v1, const cmplx<T2> &v2, cmplx<T> &res)
+  {
+  res = fwd ? cmplx<T>(v1.r*v2.r+v1.i*v2.i, v1.i*v2.r-v1.r*v2.i)
+            : cmplx<T>(v1.r*v2.r-v1.i*v2.i, v1.r*v2.i+v1.i*v2.r);
+  }
 
 template<typename T> void ROT90(cmplx<T> &a)
   { auto tmp_=a.r; a.r=-a.i; a.i=tmp_; }
@@ -615,30 +620,27 @@ template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
       for (size_t i=1; i<ido; ++i)
         {
         CH(i,k,0) = CC(i,0,k)+CC(i,1,k);
-        CH(i,k,1) = (CC(i,0,k)-CC(i,1,k)).template special_mul<fwd>(WA(0,i));
+        special_mul<fwd>(CC(i,0,k)-CC(i,1,k),WA(0,i),CH(i,k,1));
         }
       }
   }
 
 #define POCKETFFT_PREP3(idx) \
         T t0 = CC(idx,0,k), t1, t2; \
-        PMC (t1,t2,CC(idx,1,k),CC(idx,2,k)); \
+        PM (t1,t2,CC(idx,1,k),CC(idx,2,k)); \
         CH(idx,k,0)=t0+t1;
 #define POCKETFFT_PARTSTEP3a(u1,u2,twr,twi) \
         { \
-        T ca,cb; \
-        ca=t0+t1*twr; \
-        cb=t2*twi; ROT90(cb); \
-        PMC(CH(0,k,u1),CH(0,k,u2),ca,cb) ;\
+        T ca=t0+t1*twr; \
+        T cb{-t2.i*twi, t2.r*twi}; \
+        PM(CH(0,k,u1),CH(0,k,u2),ca,cb) ;\
         }
 #define POCKETFFT_PARTSTEP3b(u1,u2,twr,twi) \
         { \
-        T ca,cb,da,db; \
-        ca=t0+t1*twr; \
-        cb=t2*twi; ROT90(cb); \
-        PMC(da,db,ca,cb); \
-        CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
-        CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
+        T ca=t0+t1*twr; \
+        T cb{-t2.i*twi, t2.r*twi}; \
+        special_mul<fwd>(ca+cb,WA(u1-1,i),CH(i,k,u1)); \
+        special_mul<fwd>(ca-cb,WA(u2-1,i),CH(i,k,u2)); \
         }
 template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
@@ -697,43 +699,42 @@ template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
     for (size_t k=0; k<l1; ++k)
       {
       T t1, t2, t3, t4;
-      PMC(t2,t1,CC(0,0,k),CC(0,2,k));
-      PMC(t3,t4,CC(0,1,k),CC(0,3,k));
+      PM(t2,t1,CC(0,0,k),CC(0,2,k));
+      PM(t3,t4,CC(0,1,k),CC(0,3,k));
       ROTX90<fwd>(t4);
-      PMC(CH(0,k,0),CH(0,k,2),t2,t3);
-      PMC(CH(0,k,1),CH(0,k,3),t1,t4);
+      PM(CH(0,k,0),CH(0,k,2),t2,t3);
+      PM(CH(0,k,1),CH(0,k,3),t1,t4);
       }
   else
     for (size_t k=0; k<l1; ++k)
       {
       {
       T t1, t2, t3, t4;
-      PMC(t2,t1,CC(0,0,k),CC(0,2,k));
-      PMC(t3,t4,CC(0,1,k),CC(0,3,k));
+      PM(t2,t1,CC(0,0,k),CC(0,2,k));
+      PM(t3,t4,CC(0,1,k),CC(0,3,k));
       ROTX90<fwd>(t4);
-      PMC(CH(0,k,0),CH(0,k,2),t2,t3);
-      PMC(CH(0,k,1),CH(0,k,3),t1,t4);
+      PM(CH(0,k,0),CH(0,k,2),t2,t3);
+      PM(CH(0,k,1),CH(0,k,3),t1,t4);
       }
       for (size_t i=1; i<ido; ++i)
         {
-        T c2, c3, c4, t1, t2, t3, t4;
+        T t1, t2, t3, t4;
         T cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k);
-        PMC(t2,t1,cc0,cc2);
-        PMC(t3,t4,cc1,cc3);
+        PM(t2,t1,cc0,cc2);
+        PM(t3,t4,cc1,cc3);
         ROTX90<fwd>(t4);
-        PMC(CH(i,k,0),c3,t2,t3);
-        PMC(c2,c4,t1,t4);
-        CH(i,k,1) = c2.template special_mul<fwd>(WA(0,i));
-        CH(i,k,2) = c3.template special_mul<fwd>(WA(1,i));
-        CH(i,k,3) = c4.template special_mul<fwd>(WA(2,i));
+        CH(i,k,0) = t2+t3;
+        special_mul<fwd>(t1+t4,WA(0,i),CH(i,k,1));
+        special_mul<fwd>(t2-t3,WA(1,i),CH(i,k,2));
+        special_mul<fwd>(t1-t4,WA(2,i),CH(i,k,3));
         }
       }
   }
 
 #define POCKETFFT_PREP5(idx) \
         T t0 = CC(idx,0,k), t1, t2, t3, t4; \
-        PMC (t1,t4,CC(idx,1,k),CC(idx,4,k)); \
-        PMC (t2,t3,CC(idx,2,k),CC(idx,3,k)); \
+        PM (t1,t4,CC(idx,1,k),CC(idx,4,k)); \
+        PM (t2,t3,CC(idx,2,k),CC(idx,3,k)); \
         CH(idx,k,0).r=t0.r+t1.r+t2.r; \
         CH(idx,k,0).i=t0.i+t1.i+t2.i;
 
@@ -744,7 +745,7 @@ template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
         ca.i=t0.i+twar*t1.i+twbr*t2.i; \
         cb.i=twai*t4.r twbi*t3.r; \
         cb.r=-(twai*t4.i twbi*t3.i); \
-        PMC(CH(0,k,u1),CH(0,k,u2),ca,cb); \
+        PM(CH(0,k,u1),CH(0,k,u2),ca,cb); \
         }
 
 #define POCKETFFT_PARTSTEP5b(u1,u2,twar,twbr,twai,twbi) \
@@ -754,9 +755,8 @@ template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
         ca.i=t0.i+twar*t1.i+twbr*t2.i; \
         cb.i=twai*t4.r twbi*t3.r; \
         cb.r=-(twai*t4.i twbi*t3.i); \
-        PMC(da,db,ca,cb); \
-        CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
-        CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
+        special_mul<fwd>(ca+cb,WA(u1-1,i),CH(i,k,u1)); \
+        special_mul<fwd>(ca-cb,WA(u2-1,i),CH(i,k,u2)); \
         }
 template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
@@ -805,9 +805,9 @@ template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
 
 #define POCKETFFT_PREP7(idx) \
         T t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7; \
-        PMC (t2,t7,CC(idx,1,k),CC(idx,6,k)); \
-        PMC (t3,t6,CC(idx,2,k),CC(idx,5,k)); \
-        PMC (t4,t5,CC(idx,3,k),CC(idx,4,k)); \
+        PM (t2,t7,CC(idx,1,k),CC(idx,6,k)); \
+        PM (t3,t6,CC(idx,2,k),CC(idx,5,k)); \
+        PM (t4,t5,CC(idx,3,k),CC(idx,4,k)); \
         CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r; \
         CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i;
 
@@ -818,7 +818,7 @@ template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
         ca.i=t1.i+x1*t2.i+x2*t3.i+x3*t4.i; \
         cb.i=y1*t7.r y2*t6.r y3*t5.r; \
         cb.r=-(y1*t7.i y2*t6.i y3*t5.i); \
-        PMC(out1,out2,ca,cb); \
+        PM(out1,out2,ca,cb); \
         }
 #define POCKETFFT_PARTSTEP7a(u1,u2,x1,x2,x3,y1,y2,y3) \
         POCKETFFT_PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,CH(0,k,u1),CH(0,k,u2))
@@ -826,8 +826,8 @@ template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
         { \
         T da,db; \
         POCKETFFT_PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,da,db) \
-        CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
-        CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
+        special_mul<fwd>(da,WA(u1-1,i),CH(i,k,u1)); \
+        special_mul<fwd>(db,WA(u2-1,i),CH(i,k,u2)); \
         }
 
 template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
@@ -915,8 +915,8 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
     for (size_t k=0; k<l1; ++k)
       {
       T a0, a1, a2, a3, a4, a5, a6, a7;
-      PMC(a1,a5,CC(0,1,k),CC(0,5,k));
-      PMC(a3,a7,CC(0,3,k),CC(0,7,k));
+      PM(a1,a5,CC(0,1,k),CC(0,5,k));
+      PM(a3,a7,CC(0,3,k),CC(0,7,k));
       PMINPLACE(a1,a3);
       ROTX90<fwd>(a3);
 
@@ -925,21 +925,21 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
       ROTX45<fwd>(a5);
       ROTX135<fwd>(a7);
 
-      PMC(a0,a4,CC(0,0,k),CC(0,4,k));
-      PMC(a2,a6,CC(0,2,k),CC(0,6,k));
-      PMC(CH(0,k,0),CH(0,k,4),a0+a2,a1);
-      PMC(CH(0,k,2),CH(0,k,6),a0-a2,a3);
+      PM(a0,a4,CC(0,0,k),CC(0,4,k));
+      PM(a2,a6,CC(0,2,k),CC(0,6,k));
+      PM(CH(0,k,0),CH(0,k,4),a0+a2,a1);
+      PM(CH(0,k,2),CH(0,k,6),a0-a2,a3);
       ROTX90<fwd>(a6);
-      PMC(CH(0,k,1),CH(0,k,5),a4+a6,a5);
-      PMC(CH(0,k,3),CH(0,k,7),a4-a6,a7);
+      PM(CH(0,k,1),CH(0,k,5),a4+a6,a5);
+      PM(CH(0,k,3),CH(0,k,7),a4-a6,a7);
       }
   else
     for (size_t k=0; k<l1; ++k)
       {
       {
       T a0, a1, a2, a3, a4, a5, a6, a7;
-      PMC(a1,a5,CC(0,1,k),CC(0,5,k));
-      PMC(a3,a7,CC(0,3,k),CC(0,7,k));
+      PM(a1,a5,CC(0,1,k),CC(0,5,k));
+      PM(a3,a7,CC(0,3,k),CC(0,7,k));
       PMINPLACE(a1,a3);
       ROTX90<fwd>(a3);
 
@@ -948,38 +948,38 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
       ROTX45<fwd>(a5);
       ROTX135<fwd>(a7);
 
-      PMC(a0,a4,CC(0,0,k),CC(0,4,k));
-      PMC(a2,a6,CC(0,2,k),CC(0,6,k));
-      PMC(CH(0,k,0),CH(0,k,4),a0+a2,a1);
-      PMC(CH(0,k,2),CH(0,k,6),a0-a2,a3);
+      PM(a0,a4,CC(0,0,k),CC(0,4,k));
+      PM(a2,a6,CC(0,2,k),CC(0,6,k));
+      PM(CH(0,k,0),CH(0,k,4),a0+a2,a1);
+      PM(CH(0,k,2),CH(0,k,6),a0-a2,a3);
       ROTX90<fwd>(a6);
-      PMC(CH(0,k,1),CH(0,k,5),a4+a6,a5);
-      PMC(CH(0,k,3),CH(0,k,7),a4-a6,a7);
+      PM(CH(0,k,1),CH(0,k,5),a4+a6,a5);
+      PM(CH(0,k,3),CH(0,k,7),a4-a6,a7);
       }
       for (size_t i=1; i<ido; ++i)
         {
         T a0, a1, a2, a3, a4, a5, a6, a7;
-        PMC(a1,a5,CC(i,1,k),CC(i,5,k));
-        PMC(a3,a7,CC(i,3,k),CC(i,7,k));
+        PM(a1,a5,CC(i,1,k),CC(i,5,k));
+        PM(a3,a7,CC(i,3,k),CC(i,7,k));
         ROTX90<fwd>(a7);
         PMINPLACE(a1,a3);
         ROTX90<fwd>(a3);
         PMINPLACE(a5,a7);
         ROTX45<fwd>(a5);
         ROTX135<fwd>(a7);
-        PMC(a0,a4,CC(i,0,k),CC(i,4,k));
-        PMC(a2,a6,CC(i,2,k),CC(i,6,k));
+        PM(a0,a4,CC(i,0,k),CC(i,4,k));
+        PM(a2,a6,CC(i,2,k),CC(i,6,k));
         PMINPLACE(a0,a2);
         CH(i,k,0) = a0+a1;
-        CH(i,k,4) = (a0-a1).template special_mul<fwd>(WA(3,i));
-        CH(i,k,2) = (a2+a3).template special_mul<fwd>(WA(1,i));
-        CH(i,k,6) = (a2-a3).template special_mul<fwd>(WA(5,i));
+        special_mul<fwd>(a0-a1,WA(3,i),CH(i,k,4));
+        special_mul<fwd>(a2+a3,WA(1,i),CH(i,k,2));
+        special_mul<fwd>(a2-a3,WA(5,i),CH(i,k,6));
         ROTX90<fwd>(a6);
         PMINPLACE(a4,a6);
-        CH(i,k,1) = (a4+a5).template special_mul<fwd>(WA(0,i));
-        CH(i,k,5) = (a4-a5).template special_mul<fwd>(WA(4,i));
-        CH(i,k,3) = (a6+a7).template special_mul<fwd>(WA(2,i));
-        CH(i,k,7) = (a6-a7).template special_mul<fwd>(WA(6,i));
+        special_mul<fwd>(a4+a5,WA(0,i),CH(i,k,1));
+        special_mul<fwd>(a4-a5,WA(4,i),CH(i,k,5));
+        special_mul<fwd>(a6+a7,WA(2,i),CH(i,k,3));
+        special_mul<fwd>(a6-a7,WA(6,i),CH(i,k,7));
         }
       }
    }
@@ -987,11 +987,11 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
 
 #define POCKETFFT_PREP11(idx) \
         T t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7, t8, t9, t10, t11; \
-        PMC (t2,t11,CC(idx,1,k),CC(idx,10,k)); \
-        PMC (t3,t10,CC(idx,2,k),CC(idx, 9,k)); \
-        PMC (t4,t9 ,CC(idx,3,k),CC(idx, 8,k)); \
-        PMC (t5,t8 ,CC(idx,4,k),CC(idx, 7,k)); \
-        PMC (t6,t7 ,CC(idx,5,k),CC(idx, 6,k)); \
+        PM (t2,t11,CC(idx,1,k),CC(idx,10,k)); \
+        PM (t3,t10,CC(idx,2,k),CC(idx, 9,k)); \
+        PM (t4,t9 ,CC(idx,3,k),CC(idx, 8,k)); \
+        PM (t5,t8 ,CC(idx,4,k),CC(idx, 7,k)); \
+        PM (t6,t7 ,CC(idx,5,k),CC(idx, 6,k)); \
         CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r+t5.r+t6.r; \
         CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i+t5.i+t6.i;
 
@@ -1001,7 +1001,7 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
           cb; \
         cb.i=y1*t11.r y2*t10.r y3*t9.r y4*t8.r y5*t7.r; \
         cb.r=-(y1*t11.i y2*t10.i y3*t9.i y4*t8.i y5*t7.i ); \
-        PMC(out1,out2,ca,cb); \
+        PM(out1,out2,ca,cb); \
         }
 #define POCKETFFT_PARTSTEP11a(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \
         POCKETFFT_PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,CH(0,k,u1),CH(0,k,u2))
@@ -1009,8 +1009,8 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
         { \
         T da,db; \
         POCKETFFT_PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,da,db) \
-        CH(i,k,u1) = da.template special_mul<fwd>(WA(u1-1,i)); \
-        CH(i,k,u2) = db.template special_mul<fwd>(WA(u2-1,i)); \
+        special_mul<fwd>(da,WA(u1-1,i),CH(i,k,u1)); \
+        special_mul<fwd>(db,WA(u2-1,i),CH(i,k,u2)); \
         }
 
 template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
@@ -1105,7 +1105,7 @@ template<bool fwd, typename T> void passg (size_t ido, size_t ip,
   for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc)
     for (size_t k=0; k<l1; ++k)
       for (size_t i=0; i<ido; ++i)
-        PMC(CH(i,k,j),CH(i,k,jc),CC(i,j,k),CC(i,jc,k));
+        PM(CH(i,k,j),CH(i,k,jc),CC(i,j,k),CC(i,jc,k));
   for (size_t k=0; k<l1; ++k)
     for (size_t i=0; i<ido; ++i)
       {
@@ -1161,7 +1161,7 @@ template<bool fwd, typename T> void passg (size_t ido, size_t ip,
       for (size_t ik=0; ik<idl1; ++ik)
         {
         T t1=CX2(ik,j), t2=CX2(ik,jc);
-        PMC(CX2(ik,j),CX2(ik,jc),t1,t2);
+        PM(CX2(ik,j),CX2(ik,jc),t1,t2);
         }
   else
     {
@@ -1169,15 +1169,15 @@ template<bool fwd, typename T> void passg (size_t ido, size_t ip,
       for (size_t k=0; k<l1; ++k)
         {
         T t1=CX(0,k,j), t2=CX(0,k,jc);
-        PMC(CX(0,k,j),CX(0,k,jc),t1,t2);
+        PM(CX(0,k,j),CX(0,k,jc),t1,t2);
         for (size_t i=1; i<ido; ++i)
           {
           T x1, x2;
-          PMC(x1,x2,CX(i,k,j),CX(i,k,jc));
+          PM(x1,x2,CX(i,k,j),CX(i,k,jc));
           size_t idij=(j-1)*(ido-1)+i-1;
-          CX(i,k,j) = x1.template special_mul<fwd>(wa[idij]);
+          special_mul<fwd>(x1,wa[idij],CX(i,k,j));
           idij=(jc-1)*(ido-1)+i-1;
-          CX(i,k,jc) = x2.template special_mul<fwd>(wa[idij]);
+          special_mul<fwd>(x2,wa[idij],CX(i,k,jc));
           }
         }
     }
@@ -1333,9 +1333,6 @@ template<typename T0> class rfftp
     void add_factor(size_t factor)
       { fact.push_back({factor, nullptr, nullptr}); }
 
-template<typename T> inline void PM(T &a, T &b, T c, T d)
-  { a=c+d; b=c-d; }
-
 /* (a+ib) = conj(c+id) * (e+if) */
 template<typename T1, typename T2, typename T3> inline void MULPM
   (T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f)
@@ -1557,15 +1554,13 @@ template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
         for (size_t i=1; i<=ido-2; i+=2)                      // 112
           {
           T t1=C1(i,k,j ), t2=C1(i+1,k,j ),
-                 t3=C1(i,k,jc), t4=C1(i+1,k,jc);
+            t3=C1(i,k,jc), t4=C1(i+1,k,jc);
           T x1=wa[idij]*t1 + wa[idij+1]*t2,
-                 x2=wa[idij]*t2 - wa[idij+1]*t1,
-                 x3=wa[idij2]*t3 + wa[idij2+1]*t4,
-                 x4=wa[idij2]*t4 - wa[idij2+1]*t3;
-          C1(i  ,k,j ) = x1+x3;
-          C1(i  ,k,jc) = x2-x4;
-          C1(i+1,k,j ) = x2+x4;
-          C1(i+1,k,jc) = x3-x1;
+            x2=wa[idij]*t2 - wa[idij+1]*t1,
+            x3=wa[idij2]*t3 + wa[idij2+1]*t4,
+            x4=wa[idij2]*t4 - wa[idij2+1]*t3;
+          PM(C1(i,k,j),C1(i+1,k,jc),x3,x1);
+          PM(C1(i+1,k,j),C1(i,k,jc),x2,x4);
           idij+=2;
           idij2+=2;
           }
@@ -1575,11 +1570,7 @@ template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
 
   for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)                // 123
     for (size_t k=0; k<l1; ++k)                              // 122
-      {
-      T t1=C1(0,k,j), t2=C1(0,k,jc);
-      C1(0,k,j ) = t1+t2;
-      C1(0,k,jc) = t2-t1;
-      }
+      MPINPLACE(C1(0,k,jc), C1(0,k,j));
 
 //everything in C
 //memset(ch,0,ip*l1*ido*sizeof(double));
@@ -1961,10 +1952,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
       CH2(ik,0) += CH2(ik,j);
   for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc)   // 124
     for (size_t k=0; k<l1; ++k)
-      {
-      CH(0,k,j ) = C1(0,k,j)-C1(0,k,jc);
-      CH(0,k,jc) = C1(0,k,j)+C1(0,k,jc);
-      }
+      PM(CH(0,k,jc),CH(0,k,j),C1(0,k,j),C1(0,k,jc));
 
   if (ido==1) return;
 
@@ -2171,7 +2159,7 @@ template<typename T0> class fftblue
 
       /* initialize a_k and FFT it */
       for (size_t m=0; m<n; ++m)
-        akf[m] = c[m].template special_mul<fwd>(bk[m]);
+        special_mul<fwd>(c[m],bk[m],akf[m]);
       auto zero = akf[0]*T0(0);
       for (size_t m=n; m<n2; ++m)
         akf[m]=zero;
@@ -2366,7 +2354,7 @@ template<typename T0> class T_dct1
       for (size_t i=1; i<n; ++i)
         c[i] = tmp[2*i-1];
       if (ortho)
-        { c[0]/=sqrt2; c[n-1]/=sqrt2; }
+        { c[0]*=sqrt2*T0(0.5); c[n-1]*=sqrt2*T0(0.5); }
       }
 
     size_t length() const { return fftplan.length()/2+1; }
@@ -2439,7 +2427,7 @@ template<typename T0> class T_dcst23
         if (!cosine)
           for (size_t k=0, kc=N-1; k<kc; ++k, --kc)
             swap(c[k], c[kc]);
-        if (ortho) c[0]/=sqrt2;
+        if (ortho) c[0]*=sqrt2*T0(0.5);
         }
       else
         {
-- 
GitLab


From 8a1d38e43a6befe830dc0fbeafddbc121833e7fa Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 2 Aug 2019 18:04:01 +0200
Subject: [PATCH 2/9] make code a bit more compact

---
 pocketfft_hdronly.h | 253 ++++++++++++++++++++------------------------
 1 file changed, 113 insertions(+), 140 deletions(-)

diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h
index 480087e..e6bfd73 100644
--- a/pocketfft_hdronly.h
+++ b/pocketfft_hdronly.h
@@ -595,7 +595,7 @@ template<typename T0> class cfftp
 
 template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=2;
 
@@ -644,7 +644,7 @@ template<bool fwd, typename T> void pass2 (size_t ido, size_t l1,
         }
 template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=3;
   constexpr T0 tw1r=-0.5,
@@ -684,7 +684,7 @@ template<bool fwd, typename T> void pass3 (size_t ido, size_t l1,
 
 template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=4;
 
@@ -760,7 +760,7 @@ template<bool fwd, typename T> void pass4 (size_t ido, size_t l1,
         }
 template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=5;
   constexpr T0 tw1r= T0(0.3090169943749474241022934171828191L),
@@ -832,7 +832,7 @@ template<bool fwd, typename T> void pass5 (size_t ido, size_t l1,
 
 template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=7;
   constexpr T0 tw1r= T0(0.6234898018587335305250048840042398L),
@@ -881,7 +881,7 @@ template<bool fwd, typename T> void pass7(size_t ido, size_t l1,
 #undef POCKETFFT_PARTSTEP7a
 #undef POCKETFFT_PREP7
 
-template <bool fwd, typename T> void ROTX45(T &a)
+template <bool fwd, typename T> void ROTX45(T &a) const
   {
   constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
   if (fwd)
@@ -889,7 +889,7 @@ template <bool fwd, typename T> void ROTX45(T &a)
   else
     { auto tmp_=a.r; a.r=hsqt2*(a.r-a.i); a.i=hsqt2*(a.i+tmp_); }
   }
-template <bool fwd, typename T> void ROTX135(T &a)
+template <bool fwd, typename T> void ROTX135(T &a) const
   {
   constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
   if (fwd)
@@ -900,7 +900,7 @@ template <bool fwd, typename T> void ROTX135(T &a)
 
 template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=8;
 
@@ -1015,7 +1015,7 @@ template<bool fwd, typename T> void pass8 (size_t ido, size_t l1,
 
 template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const cmplx<T0> * POCKETFFT_RESTRICT wa)
+  const cmplx<T0> * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=11;
   constexpr T0 tw1r= T0(0.8412535328311811688618116489193677L),
@@ -1077,7 +1077,7 @@ template<bool fwd, typename T> void pass11 (size_t ido, size_t l1,
 template<bool fwd, typename T> void passg (size_t ido, size_t ip,
   size_t l1, T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
   const cmplx<T0> * POCKETFFT_RESTRICT wa,
-  const cmplx<T0> * POCKETFFT_RESTRICT csarr)
+  const cmplx<T0> * POCKETFFT_RESTRICT csarr) const
   {
   const size_t cdim=ip;
   size_t ipph = (ip+1)/2;
@@ -1183,7 +1183,7 @@ template<bool fwd, typename T> void passg (size_t ido, size_t ip,
     }
   }
 
-template<bool fwd, typename T> void pass_all(T c[], T0 fct)
+template<bool fwd, typename T> void pass_all(T c[], T0 fct) const
   {
   if (length==1) { c[0]*=fct; return; }
   size_t l1=1;
@@ -1232,11 +1232,8 @@ template<bool fwd, typename T> void pass_all(T c[], T0 fct)
   }
 
   public:
-    template<typename T> void forward(T c[], T0 fct)
-      { pass_all<true>(c, fct); }
-
-    template<typename T> void backward(T c[], T0 fct)
-      { pass_all<false>(c, fct); }
+    template<typename T> void exec(T c[], T0 fct, bool fwd) const
+      { fwd ? pass_all<true>(c, fct) : pass_all<false>(c, fct); }
 
   private:
     POCKETFFT_NOINLINE void factorize()
@@ -1335,12 +1332,12 @@ template<typename T0> class rfftp
 
 /* (a+ib) = conj(c+id) * (e+if) */
 template<typename T1, typename T2, typename T3> inline void MULPM
-  (T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f)
+  (T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f) const
   {  a=c*e+d*f; b=c*f-d*e; }
 
 template<typename T> void radf2 (size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=2;
 
@@ -1379,7 +1376,7 @@ template<typename T> void radf2 (size_t ido, size_t l1,
 
 template<typename T> void radf3(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=3;
   constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L);
@@ -1419,7 +1416,7 @@ template<typename T> void radf3(size_t ido, size_t l1,
 
 template<typename T> void radf4(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=4;
   constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L);
@@ -1467,7 +1464,7 @@ template<typename T> void radf4(size_t ido, size_t l1,
 
 template<typename T> void radf5(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=5;
   constexpr T0 tr11= T0(0.3090169943749474241022934171828191L),
@@ -1524,7 +1521,7 @@ template<typename T> void radf5(size_t ido, size_t l1,
 
 template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
   T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr)
+  const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const
   {
   const size_t cdim=ip;
   size_t ipph=(ip+1)/2;
@@ -1666,7 +1663,7 @@ template<typename T> void radfg(size_t ido, size_t ip, size_t l1,
 
 template<typename T> void radb2(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=2;
 
@@ -1698,7 +1695,7 @@ template<typename T> void radb2(size_t ido, size_t l1,
 
 template<typename T> void radb3(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=3;
   constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L);
@@ -1739,7 +1736,7 @@ template<typename T> void radb3(size_t ido, size_t l1,
 
 template<typename T> void radb4(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=4;
   constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L);
@@ -1792,7 +1789,7 @@ template<typename T> void radb4(size_t ido, size_t l1,
 
 template<typename T> void radb5(size_t ido, size_t l1,
   const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa)
+  const T0 * POCKETFFT_RESTRICT wa) const
   {
   constexpr size_t cdim=5;
   constexpr T0 tr11= T0(0.3090169943749474241022934171828191L),
@@ -1852,7 +1849,7 @@ template<typename T> void radb5(size_t ido, size_t l1,
 
 template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
   T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch,
-  const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr)
+  const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const
   {
   const size_t cdim=ip;
   size_t ipph=(ip+1)/ 2;
@@ -1985,7 +1982,7 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
     }
   }
 
-    template<typename T> void copy_and_norm(T *c, T *p1, size_t n, T0 fct)
+    template<typename T> void copy_and_norm(T *c, T *p1, size_t n, T0 fct) const
       {
       if (p1!=c)
         {
@@ -2002,60 +1999,51 @@ template<typename T> void radbg(size_t ido, size_t ip, size_t l1,
       }
 
   public:
-    template<typename T> void forward(T c[], T0 fct)
+    template<typename T> void exec(T c[], T0 fct, bool r2hc) const
       {
       if (length==1) { c[0]*=fct; return; }
-      size_t n=length;
-      size_t l1=n, nf=fact.size();
+      size_t n=length, nf=fact.size();
       arr<T> ch(n);
       T *p1=c, *p2=ch.data();
 
-      for(size_t k1=0; k1<nf;++k1)
-        {
-        size_t k=nf-k1-1;
-        size_t ip=fact[k].fct;
-        size_t ido=n / l1;
-        l1 /= ip;
-        if(ip==4)
-          radf4(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==2)
-          radf2(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==3)
-          radf3(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==5)
-          radf5(ido, l1, p1, p2, fact[k].tw);
-        else
-          { radfg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws); swap (p1,p2); }
-        swap (p1,p2);
-        }
-      copy_and_norm(c,p1,n,fct);
-      }
-
-    template<typename T> void backward(T c[], T0 fct)
-      {
-      if (length==1) { c[0]*=fct; return; }
-      size_t n=length;
-      size_t l1=1, nf=fact.size();
-      arr<T> ch(n);
-      T *p1=c, *p2=ch.data();
+      if (r2hc)
+        for(size_t k1=0, l1=n; k1<nf;++k1)
+          {
+          size_t k=nf-k1-1;
+          size_t ip=fact[k].fct;
+          size_t ido=n / l1;
+          l1 /= ip;
+          if(ip==4)
+            radf4(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==2)
+            radf2(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==3)
+            radf3(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==5)
+            radf5(ido, l1, p1, p2, fact[k].tw);
+          else
+            { radfg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws); swap (p1,p2); }
+          swap (p1,p2);
+          }
+      else
+        for(size_t k=0, l1=1; k<nf; k++)
+          {
+          size_t ip = fact[k].fct,
+                 ido= n/(ip*l1);
+          if(ip==4)
+            radb4(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==2)
+            radb2(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==3)
+            radb3(ido, l1, p1, p2, fact[k].tw);
+          else if(ip==5)
+            radb5(ido, l1, p1, p2, fact[k].tw);
+          else
+            radbg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws);
+          swap (p1,p2);
+          l1*=ip;
+          }
 
-      for(size_t k=0; k<nf; k++)
-        {
-        size_t ip = fact[k].fct,
-               ido= n/(ip*l1);
-        if(ip==4)
-          radb4(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==2)
-          radb2(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==3)
-          radb3(ido, l1, p1, p2, fact[k].tw);
-        else if(ip==5)
-          radb5(ido, l1, p1, p2, fact[k].tw);
-        else
-          radbg(ido, ip, l1, p1, p2, fact[k].tw, fact[k].tws);
-        swap (p1,p2);
-        l1*=ip;
-        }
       copy_and_norm(c,p1,n,fct);
       }
 
@@ -2153,7 +2141,7 @@ template<typename T0> class fftblue
     arr<cmplx<T0>> mem;
     cmplx<T0> *bk, *bkf;
 
-    template<bool fwd, typename T> void fft(cmplx<T> c[], T0 fct)
+    template<bool fwd, typename T> void fft(cmplx<T> c[], T0 fct) const
       {
       arr<cmplx<T>> akf(n2);
 
@@ -2164,14 +2152,14 @@ template<typename T0> class fftblue
       for (size_t m=n; m<n2; ++m)
         akf[m]=zero;
 
-      plan.forward (akf.data(),1.);
+      plan.exec (akf.data(),1.,true);
 
       /* do the convolution */
       for (size_t m=0; m<n2; ++m)
         akf[m] = akf[m].template special_mul<!fwd>(bkf[m]);
 
       /* inverse FFT */
-      plan.backward (akf.data(),1.);
+      plan.exec (akf.data(),1.,false);
 
       /* multiply by b_k */
       for (size_t m=0; m<n; ++m)
@@ -2203,40 +2191,38 @@ template<typename T0> class fftblue
         bkf[m] = bkf[n2-m] = bk[m]*xn2;
       for (size_t m=n;m<=(n2-n);++m)
         bkf[m].Set(0.,0.);
-      plan.forward(bkf,1.);
+      plan.exec(bkf,1.,true);
       }
 
-    template<typename T> void backward(cmplx<T> c[], T0 fct)
-      { fft<false>(c,fct); }
+    template<typename T> void exec(cmplx<T> c[], T0 fct, bool fwd) const
+      { fwd ? fft<true>(c,fct) : fft<false>(c,fct); }
 
-    template<typename T> void forward(cmplx<T> c[], T0 fct)
-      { fft<true>(c,fct); }
-
-    template<typename T> void backward_r(T c[], T0 fct)
-      {
-      arr<cmplx<T>> tmp(n);
-      tmp[0].Set(c[0],c[0]*0);
-      memcpy (reinterpret_cast<void *>(tmp.data()+1),
-              reinterpret_cast<void *>(c+1), (n-1)*sizeof(T));
-      if ((n&1)==0) tmp[n/2].i=T0(0)*c[0];
-      for (size_t m=1; 2*m<n; ++m)
-        tmp[n-m].Set(tmp[m].r, -tmp[m].i);
-      fft<false>(tmp.data(),fct);
-      for (size_t m=0; m<n; ++m)
-        c[m] = tmp[m].r;
-      }
-
-    template<typename T> void forward_r(T c[], T0 fct)
+    template<typename T> void exec_r(T c[], T0 fct, bool fwd)
       {
       arr<cmplx<T>> tmp(n);
-      auto zero = T0(0)*c[0];
-      for (size_t m=0; m<n; ++m)
-        tmp[m].Set(c[m], zero);
-      fft<true>(tmp.data(),fct);
-      c[0] = tmp[0].r;
-      memcpy (c+1, tmp.data()+1, (n-1)*sizeof(T));
+      if (fwd)
+        {
+        auto zero = T0(0)*c[0];
+        for (size_t m=0; m<n; ++m)
+          tmp[m].Set(c[m], zero);
+        fft<true>(tmp.data(),fct);
+        c[0] = tmp[0].r;
+        memcpy (c+1, tmp.data()+1, (n-1)*sizeof(T));
+        }
+      else
+        {
+        tmp[0].Set(c[0],c[0]*0);
+        memcpy (reinterpret_cast<void *>(tmp.data()+1),
+                reinterpret_cast<void *>(c+1), (n-1)*sizeof(T));
+        if ((n&1)==0) tmp[n/2].i=T0(0)*c[0];
+        for (size_t m=1; 2*m<n; ++m)
+          tmp[n-m].Set(tmp[m].r, -tmp[m].i);
+        fft<false>(tmp.data(),fct);
+        for (size_t m=0; m<n; ++m)
+          c[m] = tmp[m].r;
+        }
       }
-    };
+  };
 
 //
 // flexible (FFTPACK/Bluestein) complex 1D transform
@@ -2269,11 +2255,8 @@ template<typename T0> class pocketfft_c
         packplan=unique_ptr<cfftp<T0>>(new cfftp<T0>(length));
       }
 
-    template<typename T> POCKETFFT_NOINLINE void backward(cmplx<T> c[], T0 fct) const
-      { packplan ? packplan->backward(c,fct) : blueplan->backward(c,fct); }
-
-    template<typename T> POCKETFFT_NOINLINE void forward(cmplx<T> c[], T0 fct) const
-      { packplan ? packplan->forward(c,fct) : blueplan->forward(c,fct); }
+    template<typename T> POCKETFFT_NOINLINE void exec(cmplx<T> c[], T0 fct, bool fwd) const
+      { packplan ? packplan->exec(c,fct,fwd) : blueplan->exec(c,fct,fwd); }
 
     size_t length() const { return len; }
   };
@@ -2309,17 +2292,8 @@ template<typename T0> class pocketfft_r
         packplan=unique_ptr<rfftp<T0>>(new rfftp<T0>(length));
       }
 
-    template<typename T> POCKETFFT_NOINLINE void backward(T c[], T0 fct) const
-      {
-      packplan ? packplan->backward(c,fct)
-               : blueplan->backward_r(c,fct);
-      }
-
-    template<typename T> POCKETFFT_NOINLINE void forward(T c[], T0 fct) const
-      {
-      packplan ? packplan->forward(c,fct)
-               : blueplan->forward_r(c,fct);
-      }
+    template<typename T> POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool fwd) const
+      { packplan ? packplan->exec(c,fct,fwd) : blueplan->exec_r(c,fct,fwd); }
 
     size_t length() const { return len; }
   };
@@ -2349,7 +2323,7 @@ template<typename T0> class T_dct1
       tmp[0] = c[0];
       for (size_t i=1; i<n; ++i)
         tmp[i] = tmp[N-i] = c[i];
-      fftplan.forward(tmp.data(), fct);
+      fftplan.exec(tmp.data(), fct, true);
       c[0] = tmp[0];
       for (size_t i=1; i<n; ++i)
         c[i] = tmp[2*i-1];
@@ -2377,7 +2351,7 @@ template<typename T0> class T_dst1
       tmp[0] = tmp[n+1] = c[0]*0;
       for (size_t i=0; i<n; ++i)
         { tmp[i+1]=c[i]; tmp[N-1-i]=-c[i]; }
-      fftplan.forward(tmp.data(), fct);
+      fftplan.exec(tmp.data(), fct, true);
       for (size_t i=0; i<n; ++i)
         c[i] = -tmp[2*i+2];
       }
@@ -2415,7 +2389,7 @@ template<typename T0> class T_dcst23
         if ((N&1)==0) c[N-1]*=2;
         for (size_t k=1; k<N-1; k+=2)
           MPINPLACE(c[k+1], c[k]);
-        fftplan.backward(c, fct);
+        fftplan.exec(c, fct, false);
         for (size_t k=1, kc=N-1; k<NS2; ++k, --kc)
           {
           T t1 = twiddle[k-1]*c[kc]+twiddle[kc-1]*c[k];
@@ -2443,7 +2417,7 @@ template<typename T0> class T_dcst23
           }
         if ((N&1)==0)
           c[NS2] *= 2*twiddle[NS2-1];
-        fftplan.forward(c, fct);
+        fftplan.exec(c, fct, true);
         for (size_t k=1; k<N-1; k+=2)
           MPINPLACE(c[k], c[k+1]);
         if (!cosine)
@@ -2511,7 +2485,7 @@ template<typename T0> class T_dcst4
         for (; i<N; ++i, m+=4)
           y[i] = c[m];
         }
-        rfft->forward(y.data(), fct);
+        rfft->exec(y.data(), fct, true);
         {
         size_t i=0;
         for (; i+i+1<n2; ++i)
@@ -2544,7 +2518,7 @@ template<typename T0> class T_dcst4
           y[i].Set(c[2*i],c[N-1-2*i]);
           y[i] *= C2[i];
           }
-        fft->forward(y.data(), fct);
+        fft->exec(y.data(), fct, true);
         for(size_t i=0, ic=n2-1; i<n2; ++i, --ic)
           {
           c[2*i  ] =  2*(y[i ].r*C2[i ].r-y[i ].i*C2[i ].i);
@@ -2901,7 +2875,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
             tdatav[i].r[j] = tin[it.iofs(j,i)].r;
             tdatav[i].i[j] = tin[it.iofs(j,i)].i;
             }
-        forward ? plan->forward (tdatav, fct) : plan->backward(tdatav, fct);
+        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]);
@@ -2916,7 +2890,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c(
       if (buf != &tin[it.iofs(0)])
         for (size_t i=0; i<len; ++i)
           buf[i] = tin[it.iofs(i)];
-      forward ? plan->forward (buf, fct) : plan->backward(buf, fct);
+      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];
@@ -2956,7 +2930,7 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
         for (size_t i=0; i<len; ++i)
           for (size_t j=0; j<vlen; ++j)
             tdatav[i][j] = tin[it.iofs(j,i)];
-        plan->forward(tdatav, fct);
+        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;
@@ -2977,7 +2951,7 @@ template<typename T> POCKETFFT_NOINLINE void general_hartley(
       auto tdata = reinterpret_cast<T *>(storage.data());
       for (size_t i=0; i<len; ++i)
         tdata[i] = tin[it.iofs(i)];
-      plan->forward(tdata, fct);
+      plan->exec(tdata, fct, true);
       // Hartley order
       out[it.oofs(0)] = tdata[0];
       size_t i=1, i1=1, i2=len-1;
@@ -3072,7 +3046,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
       for (size_t i=0; i<len; ++i)
         for (size_t j=0; j<vlen; ++j)
           tdatav[i][j] = in[it.iofs(j,i)];
-      plan->forward(tdatav, fct);
+      plan->exec(tdatav, fct, true);
       for (size_t j=0; j<vlen; ++j)
         out[it.oofs(j,0)].Set(tdatav[0][j]);
       size_t i=1, ii=1;
@@ -3095,7 +3069,7 @@ template<typename T> POCKETFFT_NOINLINE void general_r2c(
     auto tdata = reinterpret_cast<T *>(storage.data());
     for (size_t i=0; i<len; ++i)
       tdata[i] = in[it.iofs(i)];
-    plan->forward(tdata, fct);
+    plan->exec(tdata, fct, true);
     out[it.oofs(0)].Set(tdata[0]);
     size_t i=1, ii=1;
     if (forward)
@@ -3145,7 +3119,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
         for (size_t j=0; j<vlen; ++j)
           tdatav[i][j] = in[it.iofs(j,ii)].r;
       }
-      plan->backward(tdatav, fct);
+      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];
@@ -3167,7 +3141,7 @@ template<typename T> POCKETFFT_NOINLINE void general_c2r(
     if (i<len)
       tdata[i] = in[it.iofs(ii)].r;
     }
-    plan->backward(tdata, fct);
+    plan->exec(tdata, fct, false);
     for (size_t i=0; i<len; ++i)
       out[it.oofs(i)] = tdata[i];
     }
@@ -3207,8 +3181,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];
-        forward ? plan->forward (tdatav, fct)
-                : plan->backward(tdatav, fct);
+        plan->exec (tdatav, fct, forward);
         if (r2c && (!forward))
           for (size_t i=2; i<len; i+=2)
             tdatav[i] = -tdatav[i];
@@ -3229,7 +3202,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];
-      forward ? plan->forward(buf, fct) : plan->backward(buf, fct);
+      plan->exec(buf, fct, forward);
       if (r2c && (!forward))
         for (size_t i=2; i<len; i+=2)
           buf[i] = -buf[i];
-- 
GitLab


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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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