#include "fourier.h"

namespace mimas {

template<>
void dft_< float >::operator()( int rank, const int *n,
                                const std::complex< float > *in,
                                std::complex< float > *out, int sign )
{
  fftwf_plan plan =
    fftwf_plan_dft( rank, n, (fftwf_complex *)in, (fftwf_complex *)out,
                    sign, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT );

  assert( plan != NULL );
  
  fftwf_execute( plan );
  fftwf_destroy_plan( plan );
}

template<>
void dft_< double >::operator()( int rank, const int *n,
                                 const std::complex< double > *in,
                                 std::complex< double > *out, int sign )
{
  fftw_plan plan =
    fftw_plan_dft( rank, n, (fftw_complex *)in, (fftw_complex *)out,
                   sign, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT );

  assert( plan != NULL );
  
  fftw_execute( plan );
  fftw_destroy_plan( plan );
}

template<>
void dft_r2c_< float >::operator()( int rank, const int *n,
                                    const float *in,
                                    std::complex< float > *out )
{
  fftwf_plan plan =
    fftwf_plan_dft_r2c( rank, n, (float *)in, (fftwf_complex *)out,
                        FFTW_ESTIMATE | FFTW_PRESERVE_INPUT );

  assert( plan != NULL );
  
  fftwf_execute( plan );
  fftwf_destroy_plan( plan );
}

template<>
void dft_r2c_< double >::operator()( int rank, const int *n,
                                     const double *in,
                                     std::complex< double > *out )
{
  assert( sizeof( std::complex< double > ) == sizeof( fftw_complex ) );

  fftw_plan plan =
    fftw_plan_dft_r2c( rank, n, (double *)in, (fftw_complex *)out,
                       FFTW_ESTIMATE | FFTW_PRESERVE_INPUT );

  assert( plan != NULL );
  
  fftw_execute( plan );
  fftw_destroy_plan( plan );
}

template<>
void dft_c2r_< float >::operator()( int rank, const int *n,
                                    const std::complex< float > *in,
                                    float *out )
{
  assert( sizeof( std::complex< float > ) == sizeof( fftwf_complex ) );
  fftwf_plan plan =
    fftwf_plan_dft_c2r( rank, n, (fftwf_complex *)in, (float *)out,
                        FFTW_ESTIMATE  );// FFTW_PRESERVE_INPUT

  assert( plan != NULL );
  
  fftwf_execute( plan );
  fftwf_destroy_plan( plan );
}

template<>
void dft_c2r_< double >::operator()( int rank, const int *n,
                                     const std::complex< double > *in,
                                     double *out )
{
  assert( sizeof( std::complex< double > ) == sizeof( fftw_complex ) );
  fftw_plan plan =
    fftw_plan_dft_c2r( rank, n, (fftw_complex *)in, (double *)out,
                       FFTW_ESTIMATE );// FFTW_PRESERVE_INPUT

  assert( plan != NULL );
  
  fftw_execute( plan );
  fftw_destroy_plan( plan );
}

};