#include "FFTW/fftw3.h"


static void FillSizes4FFTW(const BaseArray & in, int * sz)
{
  int k1 = 0;
  int k2 = 0;
  for(k1=in.NbDimensions()-1; k1>=0; k1--) {
    sz[k2] = in.Size(k1); k2++; 
  }
}

/* --Methode-- */
//! Constructor - If preserve_input==true, input arrays will not be overwritten.
FFTWServer::FFTWServer(bool preserve_input)
  : FFTServerInterface("FFTServer using FFTW3 package")
  , ckR8("FFTWServer - ", true, false),
  _preserve_input(preserve_input)  
{
}


/* --Methode-- */
FFTWServer::~FFTWServer()
{
}

/* --Methode-- */
FFTServerInterface * FFTWServer::Clone()
{
  return (new FFTWServer) ;
}

/* --Methode-- */
void 
FFTWServer::FFTForward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
{
  int rank = ckR8.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftw_plan p = fftw_plan_dft_1d(in.Size(), (fftw_complex *)in.Data(), 
			 (fftw_complex *)out.Data(), 
			 FFTW_FORWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftw_plan p = fftw_plan_dft(rank, sz,
		      (fftw_complex *)in.Data(), (fftw_complex *)out.Data(), 
		      FFTW_FORWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }  
  if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
  return;
}

/* --Methode-- */
void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
{
  int rank = ckR8.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftw_plan p = fftw_plan_dft_1d(in.Size(), (fftw_complex *)in.Data(), 
				   (fftw_complex *)out.Data(), 
				   FFTW_BACKWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftw_plan p = fftw_plan_dft(rank, sz,
				(fftw_complex *)in.Data(), (fftw_complex *)out.Data(), 
				FFTW_BACKWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }  

  return;
}


void FFTWServer::FFTForward(TArray< r_8 > & in, TArray< complex<r_8> > & out)
{  
  int rank = ckR8.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftw_plan p = fftw_plan_dft_r2c_1d(in.Size(), in.Data(), 
				       (fftw_complex *)out.Data(),
				       FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);

    
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftw_plan p = fftw_plan_dft_r2c(rank, sz, in.Data(),
				      (fftw_complex *)out.Data(),
				      FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }
  if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
  return;

}



void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< r_8 > & out,
			     bool usoutsz)
{
  // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
  // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
  int rank = ckR8.CheckResize(in, out, usoutsz);
  bool share = (_preserve_input) ? false : true;
  TArray< complex<r_8> > inp(in, share);
  if (rank == 1) { // One dimensional transform 
    fftw_plan p = fftw_plan_dft_c2r_1d(out.Size(), (fftw_complex *)inp.Data(), 
			     out.Data(),
			     FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(out, sz);
    fftw_plan p = fftw_plan_dft_c2r(rank, sz, (fftw_complex *)inp.Data(), 
			  out.Data(),
			  FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan"); 
    fftw_execute(p); 
    fftw_destroy_plan(p);
  }
  return;
}



/* --Methode-- */
void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa, 
			       bool usz)
{
  TVector< complex<r_8> > in(ina);
  TVector< r_8> out(outa);
  int n = in.NElts();
  r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
  sa_size_t ncs;
  if (usz) {
    if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
      throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
    ncs = out.NElts();
  }
  else {
    ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
                    ncs = 2*n-1 : ncs = 2*n-2;

    if (out.NElts() != ncs) {
      cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs 
	   << " out.NElts()= " << out.NElts() << endl;
      throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
    }
  }
  //  if (ncs == 2*n-2)  {
  out(0) = in(0).real();
  for(int k=1; k<n; k++) {
      out(k) = in(k).real();
      out(ncs-k) = in(k).imag();
  }
  if (ncs == 2*n-2)
    out(n-1) = in(n-1).real();
  
  //  for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
}


/* --Methode-- */
void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
{
  TVector< r_8> in(ina);
  TVector< complex<r_8> > out(outa);

  sa_size_t n = in.NElts();
  sa_size_t ncs = n/2+1;
  sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
  if (out.NElts() != ncs)
    throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");

  out(0) = complex<r_8> (in(0),0.);
  for(int k=1; k<ncs; k++) 
    out(k) = complex<r_8>(in(k), in(n-k));
  if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);

}