switch to new FFT

Makefile.am

@@ -3,8 +3,6 @@ ACLOCAL_AMFLAGS = -I m4
 lib_LTLIBRARIES = libsharp2.la
 
 libsharp2_la_SOURCES = \
-  libsharp2/pocketfft.cc \
-  libsharp2/pocketfft.h \
   libsharp2/sharp_utils.cc \
   libsharp2/sharp_utils.h \
   libsharp2/sharp.cc \

libsharp2/pocketfft.h

@@ -1,50 +0,0 @@
-/*
- * This file is part of pocketfft.
- * Licensed under a 3-clause BSD style license - see LICENSE.md
- */
-
-/*! \file pocketfft.h
- *  Public interface of the pocketfft library
- *
- *  Copyright (C) 2008-2019 Max-Planck-Society
- *  \author Martin Reinecke
- */
-
-#ifndef POCKETFFT_H
-#define POCKETFFT_H
-
-#include <stdlib.h>
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-struct pocketfft_plan_c_i;
-typedef struct pocketfft_plan_c_i * pocketfft_plan_c;
-pocketfft_plan_c pocketfft_make_plan_c (size_t length);
-void pocketfft_delete_plan_c (pocketfft_plan_c plan);
-int pocketfft_backward_c(pocketfft_plan_c plan, double c[], double fct);
-int pocketfft_forward_c(pocketfft_plan_c plan, double c[], double fct);
-size_t pocketfft_length_c(pocketfft_plan_c plan);
-
-struct pocketfft_plan_r_i;
-typedef struct pocketfft_plan_r_i * pocketfft_plan_r;
-pocketfft_plan_r pocketfft_make_plan_r (size_t length);
-void pocketfft_delete_plan_r (pocketfft_plan_r plan);
-/*! Computes a real backward FFT on \a c, using \a plan
-    and assuming the FFTPACK storage scheme:
-    - on entry, \a c has the form <tt>r0, r1, i1, r2, i2, ...</tt>
-    - on exit, it has the form <tt>r0, r1, ..., r[length-1]</tt>. */
-int pocketfft_backward_r(pocketfft_plan_r plan, double c[], double fct);
-/*! Computes a real forward FFT on \a c, using \a plan
-    and assuming the FFTPACK storage scheme:
-    - on entry, \a c has the form <tt>r0, r1, ..., r[length-1]</tt>;
-    - on exit, it has the form <tt>r0, r1, i1, r2, i2, ...</tt> */
-int pocketfft_forward_r(pocketfft_plan_r plan, double c[], double fct);
-size_t pocketfft_length_r(pocketfft_plan_r plan);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif

libsharp2/sharp.cc

@@ -25,9 +25,10 @@
  *  \author Martin Reinecke \author Dag Sverre Seljebotn
  */
 
-#include <math.h>
-#include <string.h>
-#include "libsharp2/pocketfft.h"
+#include <cmath>
+#include <cstring>
+#include <atomic>
+#include "mr_util/fft.h"
 #include "libsharp2/sharp_ylmgen_c.h"
 #include "libsharp2/sharp_internal.h"
 #include "libsharp2/sharp_utils.h"
@@ -78,7 +79,7 @@ typedef struct
   double phi0_;
   dcmplx *shiftarr;
   int s_shift;
-  pocketfft_plan_r plan;
+  mr::detail_fft::rfftp<double> *plan;
   int length;
   int norot;
   } ringhelper;
@@ -91,7 +92,7 @@ static void ringhelper_init (ringhelper *self)
 
 static void ringhelper_destroy (ringhelper *self)
   {
-  if (self->plan) pocketfft_delete_plan_r(self->plan);
+  delete self->plan;
   DEALLOC(self->shiftarr);
   ringhelper_init(self);
   }
@@ -114,8 +115,8 @@ MRUTIL_NOINLINE static void ringhelper_update (ringhelper *self, int nph, int mm
 //  if (!self->plan) self->plan=pocketfft_make_plan_r(nph);
   if (nph!=(int)self->length)
     {
-    if (self->plan) pocketfft_delete_plan_r(self->plan);
-    self->plan=pocketfft_make_plan_r(nph);
+    if (self->plan) delete self->plan;
+    self->plan=new mr::detail_fft::rfftp<double>(nph);
     self->length=nph;
     }
   }
@@ -332,7 +333,7 @@ MRUTIL_NOINLINE static void ringhelper_phase2ring (ringhelper *self,
       }
     }
   data[1]=data[0];
-  pocketfft_backward_r (self->plan, &(data[1]), 1.);
+  self->plan->exec(&(data[1]), 1., false);
   }
 
 MRUTIL_NOINLINE static void ringhelper_ring2phase (ringhelper *self,
@@ -351,7 +352,7 @@ MRUTIL_NOINLINE static void ringhelper_ring2phase (ringhelper *self,
   if (flags&SHARP_REAL_HARMONICS)
     wgt *= sqrt_two;
 
-  pocketfft_forward_r (self->plan, &(data[1]), 1.);
+  self->plan->exec (&(data[1]), 1., true);
   data[0]=data[1];
   data[1]=data[nph+1]=0.;
 

libsharp2/sharp_geomhelpers.cc

@@ -29,7 +29,7 @@
 #include "libsharp2/sharp_geomhelpers.h"
 #include "libsharp2/sharp_legendre_roots.h"
 #include "libsharp2/sharp_utils.h"
-#include "libsharp2/pocketfft.h"
+#include "mr_util/fft.h"
 
 void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
   const int *rings, const double *weight, sharp_geom_info **geom_info)
@@ -155,9 +155,7 @@ void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
     weight[2*k  ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
     }
   if ((nrings&1)==0) weight[nrings-1]=0.;
-  pocketfft_plan_r plan = pocketfft_make_plan_r(nrings);
-  pocketfft_backward_r(plan,weight,1.);
-  pocketfft_delete_plan_r(plan);
+  mr::r2r_fftpack({size_t(nrings)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight, weight, 1.);
 
   for (int m=0; m<(nrings+1)/2; ++m)
     {
@@ -202,9 +200,7 @@ void sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
   for (int k=1; k<=(n/2-1); ++k)
     weight[2*k-1]=2./(1.-4.*k*k) + dw;
   weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
-  pocketfft_plan_r plan = pocketfft_make_plan_r(n);
-  pocketfft_backward_r(plan,weight,1.);
-  pocketfft_delete_plan_r(plan);
+  mr::r2r_fftpack({size_t(n)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight, weight, 1.);
   weight[n]=weight[0];
 
   for (int m=0; m<(nrings+1)/2; ++m)
@@ -250,9 +246,7 @@ void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
   for (int k=1; k<=(n/2-1); ++k)
     weight[2*k-1]=2./(1.-4.*k*k);
   weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1.;
-  pocketfft_plan_r plan = pocketfft_make_plan_r(n);
-  pocketfft_backward_r(plan,weight,1.);
-  pocketfft_delete_plan_r(plan);
+  mr::r2r_fftpack({size_t(n)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight, weight, 1.);
   for (int m=0; m<nrings; ++m)
     weight[m]=weight[m+1];