source: Sophya/trunk/SophyaExt/IFFTW/fftw3server.cc@ 3593

#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)
#ifdef ALSO_FFTW_FLOAT_EXTSOP 
  : FFTServerInterface("FFTServer using FFTW3 package (r_4 and r_8)") ,
#else
  : FFTServerInterface("FFTServer using FFTW3 package (r_8 only)") ,
#endif
  ckR8("FFTWServer - ", true, false),
  ckR4("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(in.NbDimensions(), 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(in.NbDimensions(), 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(in.NbDimensions(), 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(out.NbDimensions(), 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;
  //unused: 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.);

}


#ifdef ALSO_FFTW_FLOAT_EXTSOP 
 /* ---------------------------------------------------------------------------
   We define here single precision (float) version of FFTWServr methods 
   Needs the libfftw3f.a , in addition to libfftw3.a   
   --------------------------------------------------------------------------- */
/* --Methode-- */
void 
FFTWServer::FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
{
  int rank = ckR4.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(), 
                         (fftwf_complex *)out.Data(), 
                         FFTW_FORWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftwf_plan p = fftwf_plan_dft(in.NbDimensions(), sz,
                      (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(), 
                      FFTW_FORWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }  
  if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
  return;
}

/* --Methode-- */
void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
{
  int rank = ckR4.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(), 
                                   (fftwf_complex *)out.Data(), 
                                   FFTW_BACKWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftwf_plan p = fftwf_plan_dft(in.NbDimensions(), sz,
                                (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(), 
                                FFTW_BACKWARD, FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }  

  return;
}


void FFTWServer::FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out)
{  
  int rank = ckR4.CheckResize(in, out);
  if (rank == 1) { // One dimensional transform 
    fftwf_plan p = fftwf_plan_dft_r2c_1d(in.Size(), in.Data(), 
                                         (fftwf_complex *)out.Data(), FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);

    
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(in, sz);
    fftwf_plan p = fftwf_plan_dft_r2c(in.NbDimensions(), sz, in.Data(),
                                      (fftwf_complex *)out.Data(), FFTW_ESTIMATE);
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }
  if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
  return;

}



void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & 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 = ckR4.CheckResize(in, out, usoutsz);
  bool share = (_preserve_input) ? false : true;
  TArray< complex<r_4> > inp(in, share);
  if (rank == 1) { // One dimensional transform 
    fftwf_plan p = fftwf_plan_dft_c2r_1d(out.Size(), (fftwf_complex *)inp.Data(), 
                                         out.Data(), FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }
  else {   // Multi dimensional 
    if (in.NbDimensions() > MAXND_FFTW) 
      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) rank > MAXND_FFTW !"); 
    int sz[MAXND_FFTW];
    FillSizes4FFTW(out, sz);
    fftwf_plan p = fftwf_plan_dft_c2r(out.NbDimensions(), sz, (fftwf_complex *)inp.Data(), 
                          out.Data(), FFTW_ESTIMATE);  
    if (p == NULL) 
      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan"); 
    fftwf_execute(p); 
    fftwf_destroy_plan(p);
  }
  return;
}



/* --Methode-- */
void FFTWServer::ReShapetoReal(TArray< complex<r_4> > const & ina, TArray< r_4 > & outa, 
                               bool usz)
{
  TVector< complex<r_4> > in(ina);
  TVector< r_4> out(outa);
  int n = in.NElts();
  r_4 thr = FFTArrayChecker<r_4>::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_4 > const & ina, TArray< complex<r_4> > & outa)
{
  TVector< r_4> in(ina);
  TVector< complex<r_4> > out(outa);

  sa_size_t n = in.NElts();
  sa_size_t ncs = n/2+1;
  //unused: 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_4> (in(0),0.);
  for(int k=1; k<ncs; k++) 
    out(k) = complex<r_4>(in(k), in(n-k));
  if (n%2 == 0) out(ncs-1) = complex<r_4>(in(n/2), 0.);

}

#endif