Added the possibility to shift to different fft backend. Now set to fftw. Seems to work properly

CMakeLists.txt

@@ -23,12 +23,27 @@ add_definitions(-DLASP_MAX_NFFT=33554432) # 2**25
 
 add_definitions(-D_REENTRANT)
 
+# ############### Choose an fft backend here
+#set(LASP_FFT_BACKEND fftpack)
+set(LASP_FFT_BACKEND "fftw")
+
+if(LASP_FFT_BACKEND STREQUAL "fftw")
+    find_library(FFTW_LIBRARY NAMES fftw3 fftw)
+    set(FFTW_LIBRARIES "${FFTW_LIBRARY}")
+    set(LASP_FFT_LIBRARY "${FFTW_LIBRARIES}")
+    add_definitions(-DLASP_FFT_BACKEND_FFTW)
+elseif(LASP_FFT_BACKEND STREQUAL "fftpack")
+    add_definitions(-DLASP_FFT_BACKEND_FFTPACK)
+    set(LASP_FFT_LIBRARY "fftpack")
+endif()
+
 
 if(LASP_FLOAT STREQUAL "double")
     add_definitions(-DLASP_FLOAT=64)
     add_definitions(-DLASP_DOUBLE_PRECISION)
 else()
     add_definitions(-DLASP_FLOAT=32)
+    add_definitions(-DLASP_SINGLE_PRECISION)
 endif(LASP_FLOAT STREQUAL "double")
 
 if(NOT DEFINED LASP_DEBUG)
@@ -70,8 +85,7 @@ endif(LASP_DEBUG)
 # The last argument here takes care of calling SIGABRT when an integer overflow
 # occures.
 ############################## General compilation flags (independent of debug mode, windows or linux)
-set(CYTHON_EXTRA_C_FLAGS "-Wno-sign-compare -Wno-cpp -Wno-implicit-fallthrough\
- -Wno-incompatible-pointer-types -Wno-strict-aliasing")
+set(CYTHON_EXTRA_C_FLAGS "-Wno-sign-compare -Wno-cpp -Wno-implicit-fallthrough -Wno-incompatible-pointer-types -Wno-strict-aliasing")
 
 # General make flags
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fPIC -std=c11  -Wall -Wextra -Wno-type-limits \

lasp/c/CMakeLists.txt

@@ -25,4 +25,4 @@ add_library(lasp_lib
 
 
 
-target_link_libraries(lasp_lib fftpack openblas)
+target_link_libraries(lasp_lib ${LASP_FFT_LIBRARY} openblas)

lasp/c/lasp_fft.c

@@ -9,12 +9,24 @@
 #include "lasp_tracer.h"
 #include "lasp_fft.h"
 #include "lasp_types.h"
-#include "fftpack.h"
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
+#include "fftpack.h"
 typedef struct Fft_s {
     us nfft;
-    vd fft_work;
-} Fft;
+    vd fft_work; // Storage memory for fftpack
+};
+#elif defined LASP_FFT_BACKEND_FFTW
+#include <fftw3.h>
+typedef struct Fft_s {
+    us nfft;
+    fftw_plan forward_plan;
+    fftw_plan reverse_plan;
+    c* complex_storage;
+    d* real_storage;
+};
+
+#endif
 
 Fft* Fft_create(const us nfft) {
     fsTRACE(15);
@@ -24,18 +36,28 @@ Fft* Fft_create(const us nfft) {
     }
 
     Fft* fft = a_malloc(sizeof(Fft));
-    if(fft==NULL) {
-        WARN("Fft allocation failed");
-        return NULL;
-    }
 
     fft->nfft = nfft;
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
     /* Initialize foreign fft lib */
     fft->fft_work = vd_alloc(2*nfft+15);
-
     npy_rffti(nfft,getvdval(&fft->fft_work,0));
     check_overflow_vx(fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+    fft->complex_storage = fftw_malloc(sizeof(c) * (nfft/2 + 1));
+    fft->real_storage = fftw_malloc(sizeof(d) * nfft);
+
+    fft->forward_plan = fftw_plan_dft_r2c_1d(nfft,
+            fft->real_storage,
+            fft->complex_storage,
+            FFTW_MEASURE);
+    fft->reverse_plan = fftw_plan_dft_c2r_1d(nfft,
+            fft->complex_storage,
+            fft->real_storage,
+            FFTW_MEASURE);
+
+#endif
 
     /* print_vd(&fft->fft_work); */
 
@@ -45,11 +67,20 @@ Fft* Fft_create(const us nfft) {
 void Fft_free(Fft* fft) {
     fsTRACE(15);
     dbgassert(fft,NULLPTRDEREF);
+#ifdef LASP_FFT_BACKEND_FFTPACK
     vd_free(&fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+    fftw_free(fft->complex_storage);
+    fftw_free(fft->real_storage);
+    fftw_destroy_plan(fft->forward_plan);
+    fftw_destroy_plan(fft->reverse_plan);
+#endif
     a_free(fft);
     feTRACE(15);
 }
+
 us Fft_nfft(const Fft* fft) {return fft->nfft;}
+
 void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
     fsTRACE(15);
     dbgassert(fft && freqdata && result,NULLPTRDEREF);
@@ -62,9 +93,10 @@ void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
             " result array");
 
 
-    d* freqdata_ptr = (d*) getvcval(freqdata,0);
     d* result_ptr = getvdval(result,0);
 
+#ifdef LASP_FFT_BACKEND_FFTPACK
+    d* freqdata_ptr = (d*) getvcval(freqdata,0);
     /* Copy freqdata, to fft_result. */
     d_copy(&result_ptr[1],&freqdata_ptr[2],nfft-1,1,1);
     result_ptr[0] = freqdata_ptr[0];
@@ -74,11 +106,22 @@ void Fft_ifft_single(const Fft* fft,const vc* freqdata,vd* result) {
             result_ptr,
             getvdval(&fft->fft_work,0));
 
+
+#elif defined LASP_FFT_BACKEND_FFTW
+    c* freqdata_ptr = (c*) getvcval(freqdata,0);
+
+    c_copy(fft->complex_storage, freqdata_ptr,nfft/2+1);
+
+    fftw_execute(fft->reverse_plan);
+    
+    d_copy(result_ptr, fft->real_storage, nfft, 1, 1);
+    
+#endif
+    check_overflow_vx(*result);
+
     /* Scale by dividing by nfft. Checked with numpy implementation
      * that this indeed needs to be done for FFTpack. */
     d_scale(result_ptr,1/((d) nfft),nfft);
-    
-    check_overflow_vx(*result);
     feTRACE(15);
 }
 void Fft_ifft(const Fft* fft,const cmat* freqdata,dmat* timedata) {
@@ -122,6 +165,8 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
     dbgassert(result->n_rows == (nfft/2+1),"Invalid number of rows in"
             " result array");
 
+
+#ifdef LASP_FFT_BACKEND_FFTPACK
     d* result_ptr = (d*) getvcval(result,0);
 
     /* Fftpack stores the data a bit strange, the resulting array
@@ -150,9 +195,19 @@ void Fft_fft_single(const Fft* fft,const vd* timedata,vc* result) {
     if(likely(nfft%2 == 0)) {
         result_ptr[nfft+1] = 0;
     }
-
-    check_overflow_vx(*result);
     check_overflow_vx(fft->fft_work);
+#elif defined LASP_FFT_BACKEND_FFTW
+
+    d* timedata_ptr = getvdval(timedata,0);
+    c* result_ptr = getvcval(result,0);
+    d_copy(fft->real_storage,timedata_ptr, nfft, 1, 1);
+
+    fftw_execute(fft->forward_plan);
+    
+    c_copy(result_ptr, fft->complex_storage, nfft/2+1);
+    
+#endif
+    check_overflow_vx(*result);
     feTRACE(15);
 
 }

lasp/c/lasp_mat.h

@@ -21,6 +21,7 @@ typedef struct {
     us stride;
     d* _data;
 } dmat;
+
 /// Dense matrix of complex floating point values
 typedef struct {
     us n_rows;

lasp/c/lasp_math_raw.h

@@ -173,6 +173,8 @@ static inline d d_dot(const d a[],const d b[],const us size){
  * @param to : Array to write to
  * @param from : Array to read from
  * @param size : Size of arrays
+ * @param to_inc : Mostly equal to 1, the stride of the array to copy to
+ * @param from_inc : Mostly equal to 1, the stride of the array to copy from
  */
 static inline void d_copy(d to[],
                           const d from[],

lasp/c/lasp_python.h

@@ -20,10 +20,11 @@
 /** 
  * Create a numpy array from an existing dmat.
  *
- * @param mat 
- * @param transfer_ownership 
+ * @param mat dmat struccture containing array data and metadata.
+ * @param transfer_ownership If set to true, Numpy array will be responsible
+ * for freeing the data.
  *
- * @return 
+ * @return Numpy array
  */
 static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
     fsTRACE(15);
@@ -50,7 +51,7 @@ static inline PyObject* dmat_to_ndarray(dmat* mat,bool transfer_ownership) {
     }
 
     // Transpose the array
-    PyObject* arr = PyArray_Transpose(arr_t,NULL);
+    PyObject* arr = PyArray_Transpose((PyArrayObject*) arr_t,NULL);
     if(!arr) {
         WARN("Array transpose failure");
         feTRACE(15);

test/fft_test.py

@@ -6,7 +6,7 @@ Created on Mon Jan 15 19:45:33 2018
 @author: anne
 """
 import numpy as np
-from lasp import Fft
+from lasp.wrappers import Fft
 
 nfft=9
 print('nfft:',nfft)