Using FFTW in OOP with multithreading

Question

Consider a container for two-dimensional complex-valued arrays

#include <vector>
#include <array>

struct Array2D {
  typedef std::array<double, 2> complex;

  int X, Y;
  std::vector<complex> values;
};

and a function that uses FFTW to compute the Fourier transform:

#include <fftw3.h>

void FourierTransform(const Array2D& source, Array2D *dest) {
  fftw_plan p = fftw_plan_dft_2d(source.X, source.Y, (fftw_complex*)source.values.data(), (fftw_complex*)dest->values.data(), FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(p);
  fftw_destroy_plan(p);
}

This is not an optimal implementation because it creates and destroys a plan every time the Fourier transform is computed, but it is convenient and easy to use. However, creating and destroying FFTW plans is not thread safe (see here), so this function should not be called simultaneously from different threads.

The problem:

I have many different arrays in my application of which I need to calculate Fourier transforms and
there are several function calls in between the parallelization of the main loop with #pragma omp parallel for and the places in the code where I need to calculate the Fourier transform.

Initializing all of these arrays along with the fftw plans before #pragma omp parallel for and then passing them as arguments to the subsequent functions would make the code a lot less legible. Is there another way that would allow me to encapsulate and hide the FFTW function calls in a thread-safe way?

Aside: you have undefined behaviour in your casts. You are reinterpret_casting std::array<double, 2> * to fftw_complex*, which is invalid. Why aren't you using std::complex<double> or _Complex? — Caleth

Jim Cownie Jim Cownie · Accepted Answer · 2021-03-15T10:22:15

The link into FFTW documentation suggests an answer, which is to ensure that the unsafe calls into FFTW are serialised. It suggests doing that using semaphores, but since you are in OpenMP you can do it with #pragma omp critical, something like this

    void FourierTransform(const Array2D& source, Array2D *dest) {
      fftw_plan p;
      #pragma omp critical (FFTW)
      {
        p = fftw_plan_dft_2d(source.X, source.Y, 
                             (fftw_complex*)source.values.data(),
                             (fftw_complex*)dest->values.data(), 
                             FFTW_FORWARD, FFTW_ESTIMATE);
      }
      fftw_execute(p);
      #pragma omp critical (FFTW)
      {
         fftw_destroy_plan(p);
      }
    }

I have used a named critical section there (with the name FFTW) so that these critical sections do not reduce parallelism against any other critical regions in the code, but are protected against each other.

More generally, it might be nice to provide wrapped, thread-safe, versions of the relevant functions and call those, rather than the unwrapped ones. That would be something like

static fftw_plan_p TS_fftw_plan_dft_2d(... arguments ... ) {
    fftw_plan p;
    #pragma omp critical (FFTW)
    {
        p = fftw_plan_dft_2d(... arguments ...);
    }
    return p;
}

You'd then simply call TS_fftw_foo instead of fftw_foo in the rest of your code.

If you want to get extreme, you'd put all of these wrappers into a header (tagging them as inline, too), and then also have macros to expand the relevant FFTW functions into the TS_fftw ones. Then including your header would do the job. That is rather more obscure and extreme, though. Just wrapping the functions you need is possibly clearer about what is going on.

Obviously all of this code is untested and uncompiled, but I hope you get the idea.

Using FFTW in OOP with multithreading

1 Answers