Problem casting STL complex<double> to fftw_complex

Question

The FFTW manual says that its fftw_complex type is bit compatible to std::complex<double> class in STL. But that doesn't work for me:

#include <complex>
#include <fftw3.h>
int main()
{
   std::complex<double> x(1,0);
   fftw_complex fx;
   fx = reinterpret_cast<fftw_complex>(x);
}

This gives me an error:

error: invalid cast from type ‘std::complex<double>’ to type ‘double [2]’

What am I doing wrong?

Krantz

Answer 1

Re-write your code as follows:

#include <complex>
#include <fftw3.h>
int main()
{
   std::complex<double> x(1,0);
   fftw_complex fx;
   memcpy( &fx, &x, sizeof( fftw_complex ) );
}

Every compiler I've used will optimise out the memcpy because it is copying a fixed, ie at compile time, amount of data.

This avoids pointer aliasing issues.

Edit: You can also avoid strict aliasing issues using a union as follows:

#include <complex>
#include <fftw3.h>
int main()
{
   union stdfftw
   {
       std::complex< double > stdc;
       fftw_complex           fftw;
   };
   std::complex<double> x(1,0);
   stdfftw u;
   u.stdc = x;
   fftw_complex fx = u.fftw;
}

Though strictly this C99 rules (Not sure about C++) are broken as reading from a different member of a union to the one written too is undefined. It works on most compilers though. Personally I prefer my original method.

Swaraj

Ok, that works for me! 3 years not coding in pure C :) Thanks a lot. – Krantz

@Krantz Is the memcpy necessary? Why couldn't you do this: std::complex<double> x(1,0); fftw_complex* fx = &x – Sempach

what about with an array of complex numbers? – Gilding

Answer 2

The idea behind bit-compatibility of fftw_complex and C99 and C++ complex types is not that they can be easily created from one another, but that all functions in FFTW that take pointers to fftw_complex can also take pointers to c++ std::complex. Therefore the best approach is probably to use std::complex<> throughout your program and only convert pointers to these values when calling FFTW functions:

std::vector<std::complex<double> > a1, a2;
....
....
fftw_plan_dft(N, reinterpret_cast<fftw_complex*>(&a1[0]),
                 reinterpret_cast<fftw_complex*>(&a2[0]),
                 FFTW_FORWARD, FFTW_ESTIMATE);
....

Gunmaker

By far the best way to go. – Audry

Answer 3

reinterpret_cast only works for pointers and references. So you'd have to do this:

#include <complex>
#include <fftw3.h>
int main()
{
   std::complex<double> x(1,0);
   fftw_complex fx(*reinterpret_cast<fftw_complex*>(&x));
}

This assumes that fftw_complex has a copy constructor. To avoid problems with strict aliasing, Goz's solution should be preferred.

Similitude

What you do above breaks strict aliasing rules and should be avoided. GCC will warn you for doing it and will potentially break if you have strict aliasing turned on. – Swaraj

@Goz: You are right. I added a warning and referred to your answer. – Peppergrass

Works for references too: fftw_complex fx(reinterpret_cast<fftw_complex&>(x)); – Deledda