#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#ifndef SINGLE
#include <dfftw.h>
#else
#include <sfftw.h>
#endif

/*
FFTW_MEASURE: find the optimal plan by actually computing several FFTs 
FFTW_ESTIMATE: do not run any FFT and provide a "reasonable" plan
FFTW_OUT_OF_PLACE: a plan assumes that the in and out are distinct 
FFTW_IN_PLACE: a plan assumes that the in and out are same
*/

int main(int argc, char *argv[]) {

  int 		N=0,nthreads=1,i,j;
  fftw_complex 	*in,*out;
  fftw_plan	fwd_plan, bck_plan;
  double	err=0.0;

  if(argc>2) {
   nthreads = atoi(argv[1]);
   N = atoi(argv[2]);
  } else {
   printf("usage: %s nthreads N\n\t nthreads = number of threads\n\t N = length of vector\n",argv[0]);
   exit;
  }

  if(nthreads<1) {
  }

  if(N<1) {
   printf("\tN = %d ... exiting\n",N);
   exit(1);
  }

  in = (fftw_complex *)malloc(N*sizeof(fftw_complex));
  out = (fftw_complex *)malloc(N*sizeof(fftw_complex));

  for (i=0; i<N; ++i) { in[i].re = 1.0; in[i].im = 0.0; }

  fwd_plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
  bck_plan = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);

/*
  for (i=0; i<N; ++i) { printf("in  %.3e %.3e \n",in[i].re, in[i].im); }
*/
  fftw_threads_one(nthreads,fwd_plan, in, out); 
/*
  for (i=0; i<N; ++i) { printf("out %.3e %.3e \n",out[i].re, out[i].im); }
*/
  fftw_threads_one(nthreads,bck_plan, out, in); 
/*
  for (i=0; i<N; ++i) { printf("inv %.3e %.3e \n",in[i].re, in[i].im); }
*/

  for (i=0; i<N; ++i) {
   err = pow((N-in[i].re),2.0) + pow(in[i].im,2.0);
   if(err > 0.0) printf("ERROR:: index=%d err=%.3e \n", i,err); 
  }



  fftw_destroy_plan(fwd_plan);
  fftw_destroy_plan(bck_plan);

  return 0;
}