/********************************************************************                                            
 * Research Computing Center
 * The University of Chicago
 * Edelstone Center
 * 6030 S. Ellis Avenue, Room 126
 * Chicago, IL 60637
 * 773-795-2667
 * rcc@uchicago.edu
 *
 * Author: Douglas Rudd
 * CNET: drudd
 *
 * Created: Tue Jul 10 
 ********************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <fftw3.h>

int main( int argc, char *argv[] ) {
	int fft_size;
	int nthreads = 1;
	fftw_plan p;
	fftw_complex *in, *out;

	if ( argc < 2 || argc > 3 ) {
		fprintf(stderr,"Usage: ./a.out N (size of fft) [nthreads]\n");
		exit(1);
	}

	/* initialize threading support */
	fftw_init_threads();

	/* allocate space and initialize real array to transform */
	fft_size = atoi(argv[1]);
	in = (fftw_complex *)fftw_malloc( sizeof(fftw_complex)*fft_size );
	out = (fftw_complex*)fftw_malloc( sizeof(fftw_complex)*fft_size );

	if ( in == NULL || out == NULL ) {
		fprintf(stderr,"Unable to allocate space for real arrays!\n");
		exit(1);
	}
	memset( in, 0, fft_size*sizeof(fftw_complex) );

	/* check if number of threads was specified */
	if ( argc == 3 ) {
        nthreads = atoi(argv[2]);
    }
	printf("Using %u threads.\n", nthreads);

	/* create fftw plan, compute fft using threads, and destroy plan */
	fftw_plan_with_nthreads(nthreads);
	p = fftw_plan_dft_1d( fft_size, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
	fftw_execute(p);
	fftw_destroy_plan(p);
	fftw_cleanup_threads();

	fftw_free(in);
	fftw_free(out);

	return 0;
}