#include "ctime.h"
#include "math.h"
#include <iostream>
#include <fstream>

#include "signal/fft.h"
#include "signal/dct.h"
#include "fftw3.h"

#define FFTW_estimate FFTW_ESTIMATE 
#define FFTW_measure  FFTW_MEASURE  

const double default_measure_time=1;
const double default_max_error=0.000;


static DenseVector<int>::self vector_pow2_size(18, "2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536 131072 262144");
static DenseVector<int>::self vector_non_pow2_size(12, "6 9 12 15 18 24 36 80 108 1000 10368 27000");

static DenseVector<MatrixSize<int> >::self matrix_pow2_size(15, "(4x4) (8x8) (8x16) (16x16) (16x32) (32x32) (32x64) (64x64) (64x128) (128x128) (128x256) (256x256)  (256x512) (512x512) (1024x1024)");
static DenseVector<MatrixSize<int> >::self matrix_non_pow2_size(19, "(5x5) (6x6) (9x9) (10x10) (12x12) (15x15) (25x24) (48x48) (60x60) (75x75) (80x80) (96x96) (100x100) (120x120) (144x144) (180x180) (240x240) (360x360) (1000x1000)");


double log2(double length)
{
  return log10(length)/log10(2.0);
}

template<class T> T random_range( T lowest_number, T highest_number)
{
	if(lowest_number > highest_number){
        swap(lowest_number, highest_number);
    }
	return lowest_number+ (highest_number-lowest_number)*rand()/(RAND_MAX+1);
}

template<class T> 
struct vector_random_function : public unary_function<int,T>
{
  typedef unary_function<int,T> base;
  typedef typename base::result_type result_type;
  typedef typename base::argument_type argument_type;
  
  result_type operator()(const argument_type &) const
  {
    return result_type(random_range(0,10));
  }
};

template<class T> 
struct vector_random_function<complex<T> > : public unary_function<int,complex<T> >
{
  typedef unary_function<int,complex<T> > base;
  typedef typename base::result_type result_type;
  typedef typename base::argument_type argument_type;

  result_type operator()(const argument_type &) const
  {
    return result_type(random_range(0,10),random_range(0,10));
  }
};

template<class T> 
struct matrix_random_function : public binary_function<int,int,T>
{
  typedef binary_function<int,int,T> base;
  typedef typename base::result_type result_type;
  typedef typename base::first_argument_type first_argument_type;
  typedef typename base::second_argument_type second_argument_type;

  result_type operator()(const first_argument_type &,const second_argument_type &) const
  {
    return result_type(random_range(0,10));
  }
};

template<class T> 
struct matrix_random_function<complex<T> > : public binary_function<int,int,complex<T> >
{
  typedef binary_function<int,int,complex<T> > base;
  typedef typename base::result_type result_type;
  typedef typename base::first_argument_type first_argument_type;
  typedef typename base::second_argument_type second_argument_type;

  result_type operator()(const first_argument_type &,const second_argument_type &) const
  {
    return result_type(random_range(0,10),random_range(0,10));
  }
};

template<class E,class Tr,class G1,class G2>
void excel_output(basic_ostream<E,Tr> &os, const Vector<G1> &X, const Matrix<G2> &M)
{
	for (int i=0; i<X.size(); ++i)
	{
		os<< X[i] << "\t";
		for(int n=0;n<M.ncols();++n)
			os<< round(M(i,n)) << "\t";
		os<< endl;
	}
}


template<class G1,class G2>
void half_to_full(const Vector<G1> &X, Vector<G2> &Y)
{
  int n=Y.size();
  Y[0]=X[0];
  for (int i=1, imax=(n+1)/2; i<imax; ++i)
  {
    Y[i]=X[i];
    Y[n-i]=conj(X[i]);
  }
  if (is_even(n)) Y[n/2]=X[n/2];
}

template<class G1,class G2>
void half_to_full(const Matrix<G1> &X, Matrix<G2> &Y)
{
  typedef typename Matrix<G1>::value_type value_type;

  int m=Y.nrows(); int m2=m/2;
  int n=Y.ncols(); int n2=n/2;

  for (int i=0; i<m; ++i)
  {
    Y(i,0)=X(i,0);
  }

  for (int j=1; j<=n2; ++j)
  {
    value_type x;
    x=X(0 ,j); Y(0,j)=x; Y(0,n-j)=conj(x);
    for (int i=1; i<=m2; ++i)
    {
      x=X(i  ,j); Y(i  ,j)=x; Y(m-i,n-j)=conj(x);
      x=X(m-i,j); Y(m-i,j)=x; Y(i  ,n-j)=conj(x);
    }    
  }
}




template<class E,class Tr,class G>
void complex_vector_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "complex vector fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
    int n = S[i];
    typedef complex<float> value_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		DenseVector<value_type> ::self Y(X.size(), 0);
		DenseVector<value_type> ::self Z(X.size(), 0);
		double t0;
		double sum;

		fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fft(X,Y);
		M(i,0)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_1d(X.size(),(fftwf_complex *)&X[0],(fftwf_complex *)&Z[0], FFTW_FORWARD, FFTW_estimate);     
		fftwf_execute(p1);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
		fftwf_destroy_plan(p1);
		M(i,1)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
    double d = dist(Y,Z);
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
    
		fftwf_plan p2 = fftwf_plan_dft_1d(X.size(),(fftwf_complex *)&X[0],(fftwf_complex *)&Z[0], FFTW_FORWARD, FFTW_measure);
		fftwf_execute(p2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
		fftwf_destroy_plan(p2);
		M(i,2)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}	
	
	for(int i=0;i<M.nrows();i++)
	{
	  int n=S[i];
    typedef complex<double> value_type;
		DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		DenseVector<value_type> ::self Y(X.size(), 0);
		DenseVector<value_type> ::self Z(X.size(), 0);
		double t0,sum;
  
		fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fft(X,Y);
		M(i,3)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		fftw_plan p1 = fftw_plan_dft_1d(X.size(),(fftw_complex *)&X[0],(fftw_complex *)&Z[0], FFTW_FORWARD, FFTW_estimate);     
		fftw_execute(p1);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
		fftw_destroy_plan(p1);
		M(i,4)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
    double d =dist(Y,Z);
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
  
		fftw_plan p2 = fftw_plan_dft_1d(X.size(),(fftw_complex *)&X[0],(fftw_complex *)&Z[0], FFTW_FORWARD, FFTW_measure);
		fftw_execute(p2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
		fftw_destroy_plan(p2);
		M(i,5)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}	
	
	os << "fft: complex-data, 1D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}


template<class E,class Tr,class G>
void complex_vector_ifft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "complex vector ifft" << endl;

  for(int i=0;i<M.nrows();i++)
  {
    int n = S[i];
    typedef complex<float> value_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0;
    double sum;

    ifft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) ifft(X,Y);
    M(i,0)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftwf_plan p1 = fftwf_plan_dft_1d(X.size(),(fftwf_complex *)&X[0],(fftwf_complex *)&Z[0], FFTW_BACKWARD, FFTW_estimate);     
    fftwf_execute(p1);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
    fftwf_destroy_plan(p1);
    M(i,1)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
    double d = dist(Y,Z/float(n));

    fftwf_plan p2 = fftwf_plan_dft_1d(X.size(),(fftwf_complex *)&X[0],(fftwf_complex *)&Z[0], FFTW_BACKWARD, FFTW_measure);
    fftwf_execute(p2);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);

    fftwf_destroy_plan(p2);
    M(i,2)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
  }	
	
  for(int i=0;i<M.nrows();i++)
  {
    int n=S[i];
    typedef complex<double> value_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    ifft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) ifft(X,Y);
    M(i,3)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_dft_1d(X.size(),(fftw_complex *)&X[0],(fftw_complex *)&Z[0], FFTW_BACKWARD, FFTW_estimate);     
    fftw_execute(p1);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
    double d =dist(Y,Z/double(n));

    fftw_plan p2 = fftw_plan_dft_1d(X.size(),(fftw_complex *)&X[0],(fftw_complex *)&Z[0], FFTW_BACKWARD, FFTW_measure);
    fftw_execute(p2);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
  }	

	os << "IFFT: complex-data, 1D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}


template<class E,class Tr,class G>
void real_vector_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
  DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "real vector fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
    int n = S[i];
    typedef complex<float> value_type;
    DenseVector<float> ::self X=dense_vector_generate(n, vector_random_function<float>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    half_real_fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
    M(i,0)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftwf_plan p1 = fftwf_plan_dft_r2c_1d(X.size(),&X[0],(fftwf_complex *)&Z[0], FFTW_estimate);     
    fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
    fftwf_destroy_plan(p1);
    M(i,1)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);


    double d = dist(sub(Y,0,n/2+1),sub(Z,0,n/2+1));
    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftwf_plan p2 = fftwf_plan_dft_r2c_1d(X.size(),&X[0],(fftwf_complex *)&Z[0], FFTW_measure);

    fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
    fftwf_destroy_plan(p2);
    M(i,2)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}	
	
  for(int i=0;i<M.nrows();i++)
  {
    int n=S[i];
    typedef complex<double> value_type;
    DenseVector<double> ::self X=dense_vector_generate(n, vector_random_function<double>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    half_real_fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
    M(i,3)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_dft_r2c_1d(X.size(),&X[0],(fftw_complex *)&Z[0], FFTW_estimate);     
    fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    double d=dist(sub(Y,0,n/2+1),sub(Z,0,n/2+1));
    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftw_plan p2 = fftw_plan_dft_r2c_1d(X.size(),&X[0],(fftw_complex *)&Z[0], FFTW_measure);
    fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
  }			
	os << "fft: real-data, 1D transforms" <<endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}

template<class E,class Tr,class G>
void full_real_vector_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
  DenseMatrix<double>::self M(S.size(),4,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "full real vector fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int n = S[i];
    typedef complex<float> value_type;
		DenseVector<float> ::self X=dense_vector_generate(n, vector_random_function<float>());
		DenseVector<value_type>::self Y(n    , 0); DenseVector<value_type>::self Y2(n, 0);
		DenseVector<value_type>::self Z(n/2+1, 0); DenseVector<value_type>::self Z2(n, 0);
		double t0,sum;

		half_real_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
		M(i,0)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		real_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_fft(X,Y2);
		M(i,1)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		fftwf_plan p1 = fftwf_plan_dft_r2c_1d(X.size(),&X[0],(fftwf_complex *)&Z[0], FFTW_estimate);     
		{ fftwf_execute(p1); half_to_full(Z,Z2); }
		fftwf_destroy_plan(p1);

    double d = dist(Y2,Z2);
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
	}
	
  for(int i=0;i<M.nrows();i++)
  {
    int n=S[i];
    typedef complex<double> value_type;
    DenseVector<double> ::self X=dense_vector_generate(n, vector_random_function<double>());
    DenseVector<value_type> ::self Y(n    , 0); DenseVector<value_type> ::self Y2(n, 0);
    DenseVector<value_type> ::self Z(n/2+1, 0); DenseVector<value_type> ::self Z2(n, 0);
    double t0,sum;

    half_real_fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
    M(i,2)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fft(X,Y2);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fft(X,Y2);
    M(i,3)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_dft_r2c_1d(X.size(),&X[0],(fftw_complex *)&Z[0], FFTW_estimate);     
    { fftw_execute(p1); half_to_full(Z,Z2); }
    fftw_destroy_plan(p1);

    double d=dist(Y2,Z2);
    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
  }
	os << "full fft: real-data, 1D transforms" <<endl;
	os << "\t" << "genial float" << "\t" << "genial+adaptor float" << "\t" << "genial double" << "\t" << "genial+adaptor double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}

template<class E,class Tr,class G>
void real_vector_ifft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
  DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "real vector ifft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int n = S[i];
    typedef complex<float> value_type;
		DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		//imag(X[0])=0; if (is_even(n)) imag(X[n/2])=0;
		X=0;
		DenseVector<float> ::self Y(X.size(), 0);
		DenseVector<float> ::self Z(X.size(), 0);
		double t0,sum;

		real_ifft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_ifft(X,Y);
		M(i,0)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_c2r_1d(X.size(),(fftwf_complex *)&X[0],&Z[0], FFTW_estimate);     
		//fftwf_plan p1 = fftwf_plan_dft_c2r_1d(X.size(),(fftwf_complex *)&X[0],&Z[0], FFTW_estimate|FFTW_PRESERVE_INPUT);     
		fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
		fftwf_destroy_plan(p1);
		M(i,1)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		double d = dist(Y,Z/float(n));
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

		fftwf_plan p2 = fftwf_plan_dft_c2r_1d(X.size(),(fftwf_complex *)&X[0],&Z[0], FFTW_measure);
		fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
		fftwf_destroy_plan(p2);
		M(i,2)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}	
	
	for(int i=0;i<M.nrows();i++)
	{
		int n=S[i];
    typedef complex<double> value_type;
		DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		//imag(X[0])=0; if (is_even(n)) imag(X[n/2])=0;		
		X=0;
		DenseVector<double> ::self Y(X.size(), 0);
		DenseVector<double> ::self Z(X.size(), 0);
		double t0,sum;

		real_ifft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_ifft(X,Y);
		M(i,3)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftw_plan p1 = fftw_plan_dft_c2r_1d(X.size(),(fftw_complex *)&X[0],&Z[0], FFTW_estimate);     
		//fftw_plan p1 = fftw_plan_dft_c2r_1d(X.size(),(fftw_complex *)&X[0],&Z[0], FFTW_estimate|FFTW_PRESERVE_INPUT);     
		fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
		fftw_destroy_plan(p1);
		M(i,4)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		double d=dist(Y,Z/double(n));
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

		fftw_plan p2 = fftw_plan_dft_c2r_1d(X.size(),(fftw_complex *)&X[0],&Z[0], FFTW_measure);
		fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
		fftw_destroy_plan(p2);
		M(i,5)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}			
	os << "ifft: real-data, 1D transforms" <<endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}


template<class E,class Tr,class G>
void complex_matrix_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "complex matrix fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<float> value_type;
		DenseMatrix<value_type> ::self X=dense_matrix_generate(m, n, matrix_random_function<value_type>());
		DenseMatrix<value_type> ::self Y(X.size(), 0);
		DenseMatrix<value_type> ::self Z(X.size(), 0);
		double t0,sum;

		fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fft(X,Y);
		M(i,0)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_2d(X.nrows(),X.ncols(),(fftwf_complex *)&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_FORWARD, FFTW_estimate);     
		fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
		fftwf_destroy_plan(p1);
		M(i,1)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
		
    double d = dist(Y,Z);
		if (d/(n*m)>=default_max_error) cout << S[i] << "\t" << d << endl;

		fftwf_plan p2 = fftwf_plan_dft_2d(X.nrows(),X.ncols(),(fftwf_complex *)&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_FORWARD, FFTW_measure);
		fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
		fftwf_destroy_plan(p2);
		M(i,2)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
	}	

  for(int i=0;i<M.nrows();i++)
  {
    int m=S[i].nrows();
    int n=S[i].ncols();
    typedef complex<double> value_type;
    DenseMatrix<value_type> ::self X=dense_matrix_generate(m, n, matrix_random_function<value_type>());
    DenseMatrix<value_type> ::self Y(X.size(), 0);
    DenseMatrix<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fft(X,Y);
    M(i,3)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_dft_2d(X.nrows(),X.ncols(),(fftw_complex *)&X(0,0),(fftw_complex *)&Z(0,0), FFTW_FORWARD, FFTW_estimate);     
    fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    double d=dist(Y,Z);
    if (d/(n*m)>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftw_plan p2 = fftw_plan_dft_2d(X.nrows(),X.ncols(),(fftw_complex *)&X(0,0),(fftw_complex *)&Z(0,0), FFTW_FORWARD, FFTW_measure);
    fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
  }		

	os << "fft: complex-data, 2D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl;
}

template<class E,class Tr,class G>
void complex_matrix_ifft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "complex matrix ifft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<float> value_type;
		DenseMatrix<value_type> ::self X=dense_matrix_generate(m, n, matrix_random_function<value_type>());
		DenseMatrix<value_type> ::self Y(X.size(), 0);
		DenseMatrix<value_type> ::self Z(X.size(), 0);
		double t0,sum;

		ifft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) ifft(X,Y);
		M(i,0)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_2d(X.nrows(),X.ncols(),(fftwf_complex *)&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_BACKWARD, FFTW_estimate);     
		fftwf_execute(p1); 
    double d = dist(Y,Z/float(m*n));
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
		fftwf_destroy_plan(p1);
		M(i,1)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
		
    
		fftwf_plan p2 = fftwf_plan_dft_2d(X.nrows(),X.ncols(),(fftwf_complex *)&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_BACKWARD, FFTW_measure);
		fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
		fftwf_destroy_plan(p2);
		M(i,2)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
		
		if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	

	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<double> value_type;
		DenseMatrix<value_type> ::self X=dense_matrix_generate(m, n, matrix_random_function<value_type>());
		DenseMatrix<value_type> ::self Y(X.size(), 0);
		DenseMatrix<value_type> ::self Z(X.size(), 0);
		double t0,sum;

		ifft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) ifft(X,Y);
		M(i,3)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

		fftw_plan p1 = fftw_plan_dft_2d(X.nrows(),X.ncols(),(fftw_complex *)&X(0,0),(fftw_complex *)&Z(0,0), FFTW_BACKWARD, FFTW_estimate);     
		fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
		fftw_destroy_plan(p1);
		M(i,4)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    double d=dist(Y,Z/double(m*n));
		if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;

		fftw_plan p2 = fftw_plan_dft_2d(X.nrows(),X.ncols(),(fftw_complex *)&X(0,0),(fftw_complex *)&Z(0,0), FFTW_BACKWARD, FFTW_measure);
		fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
		fftw_destroy_plan(p2);
		M(i,5)=5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
	}		
	
	os << "ifft: complex-data, 2D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl;
}


template<class E,class Tr,class G>
void real_matrix_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "real matrix fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
    int m=S[i].nrows();
    int n=S[i].ncols();
    typedef complex<float> value_type;
    DenseMatrix<float> ::self X=dense_matrix_generate(m, n, matrix_random_function<float>());
    //DenseMatrix<value_type> ::self Y(X.size(), 0);
    DenseMatrix<value_type> ::self Y(m,((n/2)/2+1)*2, 0);
    DenseMatrix<value_type> ::self Z(m,n/2+1, 0);
    typedef DenseMatrix<value_type>::self::index_type index_type;
    typedef DenseMatrix<value_type>::self::size_type  size_type;
    double t0,sum;

    half_real_fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
    M(i,0)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    fftwf_plan p1 = fftwf_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_estimate);     
    fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
    fftwf_destroy_plan(p1);
    M(i,1)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    double d=dist(sub(Y,index_type(0,0),size_type(m,n/2+1)),sub(Z,index_type(0,0),size_type(m,n/2+1)));
    if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftwf_plan p2 = fftwf_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_measure);
    fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
    fftwf_destroy_plan(p2);
    M(i,2)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
	}	
	
  for(int i=0;i<M.nrows();i++)
  {
    int m=S[i].nrows();
    int n=S[i].ncols();
    typedef complex<double> value_type;
    DenseMatrix<double> ::self X=dense_matrix_generate(m,n, matrix_random_function<double>());
    //DenseMatrix<value_type> ::self Y(m,n, 0);
    DenseMatrix<value_type> ::self Y(m,n/2+1, 0);
    DenseMatrix<value_type> ::self Z(m,n/2+1, 0);
    typedef DenseMatrix<value_type>::self::index_type index_type;
    typedef DenseMatrix<value_type>::self::size_type  size_type;
    double t0,sum;

    half_real_fft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
    M(i,3)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftw_complex *)&Z(0,0), FFTW_estimate);   
    fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);

    double d=dist(sub(Y,index_type(0,0),size_type(m,n/2+1)),sub(Z,index_type(0,0),size_type(m,n/2+1)));
    if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftw_plan p2 = fftw_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftw_complex *)&Z(0,0), FFTW_measure);
    fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
  }		
	
  os << "fft: real-data, 2D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl;
}


template<class E,class Tr,class G>
void real_matrix_ifft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "real matrix ifft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<float> value_type;
    typedef float real_type;
		DenseMatrix<value_type> ::self X1(m,((n/2)/2+1)*2, 0);
		DenseMatrix<value_type> ::self X2(m,n/2+1, 0);
		DenseMatrix<real_type> ::self Y(m,n, 0);
		DenseMatrix<real_type> ::self Z=dense_matrix_generate(m,n, matrix_random_function<real_type>());
    fftwf_plan p = fftwf_plan_dft_r2c_2d(m,n,&Z(0,0),(fftwf_complex *)&X2(0,0), FFTW_estimate);     
    fftwf_execute(p); fftwf_destroy_plan(p);
    sub(X1,0,0,X2.nrows(),X2.ncols())=X2;
		double t0,sum;

    real_ifft(X1,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_ifft(X1,Y);
    M(i,0)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

    double d=dist(Y,Z);
    if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;    

    fftwf_plan p1 = fftwf_plan_dft_c2r_2d(m,n,(fftwf_complex *)&X2(0,0),&Z(0,0), FFTW_estimate);     
    fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
    fftwf_destroy_plan(p1);
    M(i,1)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);

    fftwf_plan p2 = fftwf_plan_dft_c2r_2d(m,n,(fftwf_complex *)&X2(0,0),&Z(0,0), FFTW_measure);
    fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
    fftwf_destroy_plan(p2);
    M(i,2)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);
	}	

  for(int i=0;i<M.nrows();i++)
  {
    int m=S[i].nrows();
    int n=S[i].ncols();
    typedef complex<double> value_type;
    typedef double real_type;
    DenseMatrix<value_type> ::self X(m,n/2+1, 0);
    DenseMatrix<real_type> ::self Y(m,n, 0);
    //DenseMatrix<real_type> ::self Z(m,n, 0);
    DenseMatrix<real_type> ::self Z=dense_matrix_generate(m,n, matrix_random_function<real_type>());
    fftw_plan p = fftw_plan_dft_r2c_2d(m,n,&Z(0,0),(fftw_complex *)&X(0,0), FFTW_estimate);   
    fftw_execute(p); fftw_destroy_plan(p);
    double t0,sum;

    real_ifft(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_ifft(X,Y);
    M(i,3)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);

    double d=dist(Y,Z);
    if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftw_plan p1 = fftw_plan_dft_c2r_2d(m,n,(fftw_complex *)&X(0,0),&Z(0,0), FFTW_estimate);   
    fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);
         		
    fftw_plan p2 = fftw_plan_dft_c2r_2d(m,n,(fftw_complex *)&X(0,0),&Z(0,0), FFTW_measure);
    fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);
  }		
	
  os << "ifft: real-data, 2D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl;
}


template<class E,class Tr,class G>
void full_real_matrix_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),4,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "full real matrix fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<float> value_type;
		DenseMatrix<float> ::self X=dense_matrix_generate(m, n, matrix_random_function<float>());
		DenseMatrix<value_type> ::self Y(m,n    , 0); DenseMatrix<value_type> ::self Y2(m,n, 0);
		DenseMatrix<value_type> ::self Z(m,n/2+1, 0); DenseMatrix<value_type> ::self Z2(m,n, 0);
		typedef DenseMatrix<value_type>::self::index_type index_type;
		typedef DenseMatrix<value_type>::self::size_type  size_type;
		double t0,sum;

    half_real_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
		M(i,0)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);
		
    real_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_fft(X,Y2);
		M(i,1)=2.5*m*n*log2( m*n )/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftwf_complex *)&Z(0,0), FFTW_estimate);     
		{ fftwf_execute(p1); half_to_full(Z,Z2); }
		fftwf_destroy_plan(p1);
		
		double d=dist(Y2,Z2);
		if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	
	
	for(int i=0;i<M.nrows();i++)
	{
	  int m=S[i].nrows();
	  int n=S[i].ncols();
    typedef complex<double> value_type;
		DenseMatrix<double> ::self X=dense_matrix_generate(m,n, matrix_random_function<double>());
		DenseMatrix<value_type> ::self Y(m,n    , 0); DenseMatrix<value_type> ::self Y2(m,n, 0);
		DenseMatrix<value_type> ::self Z(m,n/2+1, 0); DenseMatrix<value_type> ::self Z2(m,n, 0);
		typedef DenseMatrix<value_type>::self::index_type index_type;
		typedef DenseMatrix<value_type>::self::size_type  size_type;
		double t0,sum;

    half_real_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
		M(i,2)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);

    real_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) real_fft(X,Y2);
		M(i,3)=2.5*m*n*log2(m*n)/(1000000*(get_clock()-t0)/sum);
		
		fftw_plan p1 = fftw_plan_dft_r2c_2d(X.nrows(),X.ncols(),&X(0,0),(fftw_complex *)&Z(0,0), FFTW_estimate);   
		{ fftw_execute(p1); half_to_full(Z,Z2); }
		fftw_destroy_plan(p1);
		
		double d=dist(Y2,Z2);
		if (d/(m*n)>=default_max_error) cout << S[i] << "\t" << d << endl;
	}		
	
  os << "full fft: real-data, 2D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "genial+adaptor float" << "\t" << "genial double" << "\t" << "genial+adaptor double" << endl;
	excel_output(os,S,M);
	os << endl << endl;
}




template<class E,class Tr,class G>
void real_vector_dct_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
  DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

  cout << endl << "real vector dct" << endl;

	for(int i=0;i<M.nrows();i++)
	{
    int n = S[i];
    typedef float value_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    dct(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) dct(X,Y);
    M(i,0)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftwf_plan p1 = fftwf_plan_r2r_1d(X.size(),&X[0],&Z[0],FFTW_REDFT10,FFTW_estimate);     
    fftwf_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p1);
    fftwf_destroy_plan(p1);
    M(i,1)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
    Z/=value_type(sqrt(value_type(2*n))); Z[0]/=value_type(sqrt(value_type(2))); 

    double d = dist(Y,Z);
    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftwf_plan p2 = fftwf_plan_r2r_1d(X.size(),&X[0],&Z[0],FFTW_REDFT10,FFTW_measure);
    fftwf_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftwf_execute(p2);
    fftwf_destroy_plan(p2);
    M(i,2)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}	
	
  for(int i=0;i<M.nrows();i++)
	{
    int n = S[i];
    typedef double value_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
    DenseVector<value_type> ::self Y(X.size(), 0);
    DenseVector<value_type> ::self Z(X.size(), 0);
    double t0,sum;

    dct(X,Y);
    for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) dct(X,Y);
    M(i,3)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

    fftw_plan p1 = fftw_plan_r2r_1d(X.size(),&X[0],&Z[0],FFTW_REDFT10,FFTW_estimate);     
    fftw_execute(p1); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p1);
    fftw_destroy_plan(p1);
    M(i,4)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
    Z/=value_type(sqrt(value_type(2*n))); Z[0]/=value_type(sqrt(value_type(2))); 

    double d = dist(Y,Z);
    if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;

    fftw_plan p2 = fftw_plan_r2r_1d(X.size(),&X[0],&Z[0],FFTW_REDFT10,FFTW_measure);
    fftw_execute(p2); for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) fftw_execute(p2);
    fftw_destroy_plan(p2);
    M(i,5)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
	}			
	os << "dct: real-data, 1D transforms" <<endl;
	os << "\t" << "genial float" << "\t" << "fftw estimate float" << "\t" << "fftw measure float" << "\t" << "genial double" << "\t" << "fftw estimate double" << "\t" << "fftw measure double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}

template<class E,class Tr,class G>
void complex_vector_norm_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

	cout << endl << "complex vector norm fft" << endl;

  for(int i=0;i<M.nrows();i++)
	{
	  int n = S[i];
    typedef complex<float> value_type;
    typedef float real_type;
    DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		DenseVector<value_type> ::self Y(X.size(), 0); DenseVector<real_type> ::self Y2(Y.size(), 0);
		DenseVector<value_type> ::self Z(X.size(), 0); DenseVector<real_type> ::self Z2(Z.size(), 0);
		double t0;
		double sum;

		complex_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) complex_fft(X,Y);
		M(i,0)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		norm_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) norm_fft(X,Y2);
		M(i,1)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_1d(X.size(),(fftwf_complex *)&X[0],(fftwf_complex *)&Z[0], FFTW_FORWARD, FFTW_estimate);     
		{ fftwf_execute(p1); Z2=norm(Z); }
		fftwf_destroy_plan(p1);
		
    double d = dist(Y2,Z2);
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	
	
	for(int i=0;i<M.nrows();i++)
	{
	  int n=S[i];
    typedef complex<double> value_type;
    typedef double real_type;
		DenseVector<value_type> ::self X=dense_vector_generate(n, vector_random_function<value_type>());
		DenseVector<value_type> ::self Y(X.size(), 0); DenseVector<real_type> ::self Y2(Y.size(), 0);
		DenseVector<value_type> ::self Z(X.size(), 0); DenseVector<real_type> ::self Z2(Z.size(), 0);
		double t0,sum;
		
		complex_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) complex_fft(X,Y);
		M(i,2)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		norm_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) norm_fft(X,Y2);
		M(i,3)=5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		fftw_plan p1 = fftw_plan_dft_1d(X.size(),(fftw_complex *)&X[0],(fftw_complex *)&Z[0], FFTW_FORWARD, FFTW_estimate);     
		{ fftw_execute(p1); Z2=norm(Z); }
		fftw_destroy_plan(p1);

    double d =dist(Y2,Z2);
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	
	
	os << "norm fft: complex-data, 1D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "genial+adaptor float" << "\t" << "genial double" << "\t" << "genial+adaptor double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}




template<class E,class Tr,class G>
void real_vector_norm_fft_benchmark(basic_ostream<E,Tr> &os, const Vector<G> &S, double dt=default_measure_time)
{
	DenseMatrix<double>::self M(S.size(),6,-1);
	srand(static_cast<unsigned>(time(0)));

	cout << endl << "real vector norm fft" << endl;

	for(int i=0;i<M.nrows();i++)
	{
	  int n = S[i];
    typedef float real_type;
    typedef complex<float> complex_type;
    DenseVector<real_type> ::self X=dense_vector_generate(n, vector_random_function<real_type>());
		DenseVector<complex_type> ::self Y(n    , 0); DenseVector<real_type> ::self Y2(Y.size(), 0);
		DenseVector<complex_type> ::self Z(n/2+1, 0); DenseVector<real_type> ::self Z2(Z.size(), 0);
		double t0;
		double sum;

		half_real_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
		M(i,0)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);
		
		half_real_norm_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_norm_fft(X,Y2);
		M(i,1)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftwf_plan p1 = fftwf_plan_dft_r2c_1d(X.size(),&X[0],(fftwf_complex *)&Z[0], FFTW_estimate);     
		{ fftwf_execute(p1); Z2=norm(Z); }
		fftwf_destroy_plan(p1);
		
    double d = dist(sub(Y2,0,n/2+1),sub(Z2,0,n/2+1));
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	
	
	for(int i=0;i<M.nrows();i++)
	{
	  int n=S[i];
    typedef double real_type;
    typedef complex<double> complex_type;
		DenseVector<real_type> ::self X=dense_vector_generate(n, vector_random_function<real_type>());
		DenseVector<complex_type> ::self Y(X.size()    , 0); DenseVector<real_type> ::self Y2(Y.size(), 0);
		DenseVector<complex_type> ::self Z(X.size()/2+1, 0); DenseVector<real_type> ::self Z2(Z.size(), 0);
		double t0,sum;
		
		half_real_fft(X,Y);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_fft(X,Y);
		M(i,2)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		half_real_norm_fft(X,Y2);
		for (sum=0,t0=get_clock(); (get_clock()-t0)<dt; ++sum) half_real_norm_fft(X,Y2);
		M(i,3)=2.5*n*log2(n)/(1000000*(get_clock()-t0)/sum);

		fftw_plan p1 = fftw_plan_dft_r2c_1d(X.size(),&X[0],(fftw_complex *)&Z[0], FFTW_estimate);     
		{ fftw_execute(p1); Z2=norm(Z); }
		fftw_destroy_plan(p1);
		
    double d = dist(sub(Y2,0,n/2+1),sub(Z2,0,n/2+1));
		if (d/n>=default_max_error) cout << S[i] << "\t" << d << endl;
	}	
	
	os << "norm_fft: real-data, 1D transforms" << endl;
	os << "\t" << "genial float" << "\t" << "genial+adaptor float" << "\t" << "genial double" << "\t" << "genial+adaptor double" << endl;
	excel_output(os,S,M);
	os << endl << endl << endl << endl << endl << endl;
}


int main()
{
  string s="Benchmark FFT ";
  s+=to_string(FFT_LEVEL);
#if defined(SSE3)
  s+=" SSE3";
#elif defined(SSE2)
  s+=" SSE2";
#endif
#ifdef FFT_THREADING
  #ifdef OPENMP
    s+=" OpenMP"+to_string(max_threads());
  #else
    s+=" MT"+to_string(max_threads());
  #endif
#endif
  s+=".txt";
  
	ofstream fout(s.c_str());	
  cout << s << endl;
		 
  complex_vector_fft_benchmark    (fout,vector_pow2_size);
  real_vector_fft_benchmark       (fout,vector_pow2_size);
  //full_real_vector_fft_benchmark  (fout,vector_pow2_size);
  //complex_vector_ifft_benchmark   (fout,vector_pow2_size);
  //real_vector_ifft_benchmark      (fout,vector_pow2_size);
  //real_vector_dct_benchmark       (fout,vector_pow2_size);

  complex_matrix_fft_benchmark    (fout,matrix_pow2_size);
  real_matrix_fft_benchmark       (fout,matrix_pow2_size);
  //full_real_matrix_fft_benchmark  (fout,matrix_pow2_size);
  //complex_matrix_ifft_benchmark   (fout,matrix_pow2_size);
  //real_matrix_ifft_benchmark      (fout,matrix_pow2_size);

  //complex_vector_norm_fft_benchmark(fout,vector_pow2_size);	
  //real_vector_norm_fft_benchmark   (fout,vector_pow2_size);  



  //complex_vector_fft_benchmark    (fout,vector_non_pow2_size);
  //real_vector_fft_benchmark       (fout,vector_non_pow2_size);
  //full_real_vector_fft_benchmark  (fout,vector_non_pow2_size);
  //complex_vector_ifft_benchmark   (fout,vector_non_pow2_size);
  //real_vector_dct_benchmark       (fout,vector_non_pow2_size);

  //complex_matrix_fft_benchmark    (fout,matrix_non_pow2_size);
  //real_matrix_fft_benchmark       (fout,matrix_non_pow2_size);
  //full_real_matrix_fft_benchmark  (fout,matrix_non_pow2_size);
  //complex_matrix_ifft_benchmark   (fout,matrix_non_pow2_size);


  //complex_vector_norm_fft_benchmark(fout,vector_non_pow2_size);	
  //real_vector_norm_fft_benchmark   (fout,vector_non_pow2_size);  
}