source: Sophya/trunk/SophyaProg/Tests/tfft.cc@ 1956

#include "machdefs.h"

#include <math.h>
#include <iostream.h>

#include "nbrandom.h"
#include "matharr.h"
#include "fftpserver.h"
#include "fftmserver.h"
#include "fftwserver.h"
#include "ntoolsinit.h"

#include "timing.h"


static bool inp_typ_random = false ;  // true -> random input

template <class T>
inline T module(complex<T> c) 
{
  return (sqrt(c.real()*c.real()+c.imag()*c.imag()));
}

// Max Matrix print elts
static int nprt = 2;
static int nprtfc = 8;

template <class T>
void TestFFTPack(T seuil, int num)
{
  int i;
  T fact = 1./num;

  TVector< complex<T> > inc(num), bkc(num), difc(num);
  TVector< T > in(num), ino(num), bk(num),dif(num);
  TVector< complex<T> > outc(num);

  cout << " DBG/1 outc " << outc.NElts()  << endl; 
  outc.ReSize(32);
  cout << " DBG/2 outc " << outc.NElts()  << endl; 
  outc.ReSize(10);
  cout << " DBG/3 outc " << outc.NElts()  << endl; 
  outc.ReSize(48);
  cout << " DBG/4 outc " << outc.NElts()  << endl; 

  if (inp_typ_random) 
    for (i=0; i<num ; i++){
      ino[i] = in[i] = GauRnd(0., 1.);
      inc[i] = complex<T> (in[i], 0.);
    }
  else for (i=0; i<num ; i++){
    ino[i] = in[i] = 0.5 + cos(2*M_PI*(double)i/(double)num) 
                         + 2*sin(4*M_PI*(double)i/(double)num);
    inc[i] = complex<T> (in[i], 0.);
  }
  

  cout << "Input / L = " << num << in; 
  cout << endl;

  cout << " >>>> Testing FFTPackServer "  << endl;
  FFTPackServer fftp;
  cout << " Testing FFTPackServer "  << endl;
  fftp.fftf(in.NElts(), in.Data());
  //  in /= (num/2.);
  cout << " fftp.fftf(in.NElts(), in.Data()) FORWARD: " << in << endl; 
  cout << endl;
  fftp.fftb(in.NElts(), in.Data());
  cout << " fftp.fftb(in.NElts(), in.Data()) BACKWARD: " << in <<endl; 
  cout << endl;
  dif = ino-in;
  cout << " dif , NElts= " << dif.NElts() << dif << endl;

  int ndiff = 0;
  T maxdif=0., vdif;
  for(i=0; i<num; i++) {
    vdif = fabs(dif(i));
    if (vdif > seuil)  ndiff++;
    if (vdif > maxdif) maxdif = vdif;
  }
  cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff 
       << " MaxDiff= " << maxdif << endl;

}

template <class T>
void TestFFTS(T seuil, FFTServerInterface & ffts, int num)
{

  cout <<" ===>  TestFFTS " << ffts.getInfo() << " ArrSz= " << num << endl;  
  int i;

  T fact = 1.;

  TVector< complex<T> > inc(num), bkc(num), difc(num);
  TVector< T > in(num), ino(num), bk(num),dif(num);
  TVector< complex<T> > outc(num);

  for (i=0; i<num ; i++){
    ino[i] = in[i] = 0.5 + cos(2*M_PI*(double)i/(double)num) 
                         + 2*sin(4*M_PI*(double)i/(double)num);
    inc[i] = complex<T> (in[i], 0.);
  }
  

  cout << " Testing FFTServer " << ffts.getInfo() << endl;

  cout << "Input / L = " << num << in << endl; 
  cout << endl;

  int ndiff = 0;

  cout << "\n ----  Testing FFT(T, complex<T>) ---- " << endl;
  ffts.FFTForward(in, outc);
  cout << " FourierCoefs , NElts= " << outc.NElts() << outc << endl;

  ffts.FFTBackward(outc, bk);
  cout << " Backward , NElts= " << bk.NElts() << bk << endl;
 
  dif = bk*fact - in;
  cout << " Difference , NElts= " << dif.NElts() << dif << endl;

  ndiff = 0;
  T maxdif=0., vdif;
  for(i=0; i<num; i++) {
    vdif = fabs(dif(i));
    if (vdif > seuil)  ndiff++;
    if (vdif > maxdif) maxdif = vdif;
  }
  cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff 
       << " MaxDiff= " << maxdif << endl;
 
  cout << "\n ----  Testing FFT(complex<T>, complex<T>) ---- " << endl;
  ffts.FFTForward(inc, outc);
  cout << " FourierCoef , NElts= " << outc.NElts() << outc << endl;

  ffts.FFTBackward(outc, bkc);
  cout << " Backward , NElts= " << bkc.NElts() << bkc << endl;
 
  difc = bkc*complex<T>(fact,0.) - inc;
  cout << " Difference , NElts= " << difc.NElts() << difc << endl;

  ndiff = 0;
  maxdif=0., vdif;
  for(i=0; i<num; i++) {
    vdif = fabs(module(difc(i)));
    if (vdif > seuil)  ndiff++;
    if (vdif > maxdif) maxdif = vdif;
  }
  cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff 
       << " MaxDiff= " << maxdif << endl;
}





int main(int narg, char* arg[])
{

  SophyaInit();
  InitTim();   // Initializing the CPU timer

  if (narg < 4) {
    cout << "tfft/ args error - \n  Usage tfft size px/Px/Mx/Wx f/d [NPrt=50 PrtLev=0] \n" 
         << " p=FFTPackTest  P=FFTPack, M=FFTMayer, W= FFTWServer" 
         << " x=0 -> Random input - f=float, d=double " << endl;
    exit(0);
  }

  FFTPackServer fftp;
  FFTMayerServer fftm;

  FFTWServer fftw;

  inp_typ_random = false;
  if (arg[2][1] == '0')  inp_typ_random = true;

  int sz = atoi(arg[1]);
  int nprt = 50;
  int prtlev = 0;
  if (narg > 4) nprt = atoi(arg[4]);
  if (narg > 5) prtlev = atoi(arg[5]);
  BaseArray::SetMaxPrint(nprt, prtlev);

  if (sz < 2) sz = 2;
  FFTServerInterface * ffts;
  if (*arg[2] == 'M')  ffts = fftm.Clone();
  else if (*arg[2] == 'W')  ffts = fftw.Clone();
  else ffts = fftp.Clone();

  float fs = 1.e-4;
  double ds = 1.e-6;

  try {
    if (*arg[3] == 'd') {
      cout << "\n ========================================== \n"
           << "    ------ Testing FFTServer for double ----- \n"
           << " ============================================ " << endl;
      if (*arg[2] == 'p') TestFFTPack(ds, sz);
      else TestFFTS(ds, *ffts, sz);
    }
    else {
      cout << "\n ========================================= \n"
           << "    ------ Testing FFTServer for float ----- \n"
           << " =========================================== " << endl;
      if (*arg[2] == 'p') TestFFTPack(fs, sz);
      else TestFFTS(fs, *ffts, sz);
    }
  }
  catch(PThrowable exc ) {
    cerr << "TestFFT-main() , Catched exception: \n" << exc.Msg() << endl;
  }
  catch(std::exception ex) {
    cerr << "TestFFT-main() , Catched exception ! " << (string)(ex.what()) << endl;
  }

  /*
  catch(...) {
    cerr << "TestFFT-main() , Catched ... exception ! " << endl;
  }
  */
  PrtTim("End of tfft ");
  delete ffts;
}