template_real.c

// This is the template used by PP to generate the FFTW routines.
// FFTW3.pd includes this file

// This is passed into pp_def() in the 'Code' key. Before this file is passed to
// pp_def, the following strings are replaced:
//
// INVERSE      If this is a c2r transform rather than r2c
// RANK         The rank of this transform

#ifndef __TEMPLATE_ALREADY_INCLUDED__

/* the Linux kernel does something similar to assert at compile time */
#define static_assert_fftw(x) (void)( sizeof( int[ 1 - 2* !(x) ]) )

#define __TEMPLATE_ALREADY_INCLUDED__
#endif


{
  // make sure the PDL data type I'm using matches the FFTW data type
  static_assert_fftw( sizeof($GENERIC())*2 == sizeof($TFD(fftwf_,fftw_)complex) );

  $TFD(fftwf_,fftw_)plan plan = INT2PTR( $TFD(fftwf_,fftw_)plan, $COMP(plan));

  if( !INVERSE )
    $TFD(fftwf_,fftw_)execute_dft_r2c( plan,
                                       ($TFD(float,double)*)$P(real),
                                       ($TFD(fftwf_,fftw_)complex*)$P(complexv) );
  else
  {
    // FFTW inverse real transforms clobber their input. I thus make a new
    // buffer and transform from there
    unsigned long nelem = 1;
    for( int i=0; i<=RANK; i++ )
      nelem *= $PDL(complexv)->dims[i];
    unsigned long elem_scale = sizeof($GENERIC()) / sizeof( $TFD(float,double) ); /* native complex */
    /* explicit C types here means when changed to $TGC will still be right */
    $TFD(float,double)* input_copy = $TFD(fftwf_,fftw_)alloc_real( nelem * elem_scale );
    memcpy( input_copy, $P(complexv), sizeof($GENERIC()) * nelem );

    $TFD(fftwf_,fftw_)execute_dft_c2r( plan,
                                       ($TFD(fftwf_,fftw_)complex*)input_copy,
                                       ($TFD(float,double)*)$P(real) );

    $TFD(fftwf_,fftw_)free( input_copy );
  }
}