source: Sophya/trunk/AddOn/TAcq/tmtfft.cc@ 3916

//--------------------------------------------------------------------------------
// Programme d'exemple multi-threads de FFT par paquet
//  Test BAO-Radio/Acq
//   R. Ansari - C. Magneville  ,   Decembre 2007
// Exemples de commandes pour test (dans l'ordre ci-dessous :
// csh> tmtfft VSZ NTH NbFFT [PaqSz=4096] 
// csh> tmtfft 1310720 1 10000 16384
// csh> tmtfft 1310720 1 10000 8192
// csh> tmtfft 1310720 2 10000 8192
//--------------------------------------------------------------------------------

#include "sopnamsp.h"
#include "zthread.h"
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"

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

#include "tvector.h"
#include "fioarr.h"
#include "matharr.h"
#include "tarrinit.h"

#include "randfmt.h"
#include "fftpserver.h"
#include "fftwserver.h"

#include "FFTW/fftw3.h"


//---------------------------------------------------------------------------------
// Definition d'une classe heritant de ZThread, pour remplir un vecteur avec 
// des aleatoires

class MTDoFFT : public ZThread {
public:
  MTDoFFT(TVector<r_4>& v, sa_size_t nfft=10, sa_size_t paqsz=4096);
  virtual void run(); 
  inline TVector<r_4>& getSpectre() { return spectre_; }
protected:
  string nom_;
  TVector<r_4> vv_;
  sa_size_t nbfft_;  // Nb de FFT a faire
  sa_size_t paqsz_;  // taille des paquets pour FFT
  TVector<r_4> spectre_;   // spectre moyenne 
  TVector< complex<r_4> > cfour_;  // composant TF
  FFTPackServer ffts_;   // serveur FFT
};

static int mtrandId = 0;  // Pour donner un identificateur a chaque thread
ZMutex* pmtx_=NULL;  // pour synchroniser creation de plan_fftw;

MTDoFFT::MTDoFFT(TVector<r_4>& v, sa_size_t nfft, sa_size_t paqsz)
  : vv_(v,true),    // Partage de reference du vecteur 
    nbfft_(nfft), paqsz_(paqsz), ffts_(true) // preserve input=true
{
  if (pmtx_==NULL)  pmtx_ = new ZMutex;
  char buff[32];
  sprintf(buff, "MTDoFFT-Id=%d", mtrandId);
  mtrandId++;
  nom_ = buff;
  // Initialisation taille vecteurs 
  TVector<r_4> vx(paqsz_);
  vx = RegularSequence();
  ffts_.FFTForward(vx, cfour_);
  spectre_.ReSize(cfour_.Size());
  cout << " Thread MTDoFFT(" << nom_ << ", PaqSz=" << paqsz_ << " ) Created ... " << endl;
}

inline r_4 Zmod2(complex<r_4> z) 
{ return (z.real()*z.real()+z.imag()*z.imag()); }

// Le travail fait dans chaque thread (appele par start())
void MTDoFFT::run()
{
  Timer tm(nom_.c_str());
  cout << "MTDoFFT::run() - Nom= " << nom_ << " vv.Size()= " << vv_.Size() 
       << " PaqSz=" << paqsz_ << " NbFFT= " << nbfft_ << endl;
  FMTRandGen rgen;
  TVector<r_4> vx(paqsz_);
  
  pmtx_->lock();
  fftwf_plan plan = fftwf_plan_dft_r2c_1d(paqsz_, vx.Data(), 
                                          (fftwf_complex *)cfour_.Data(), FFTW_ESTIMATE);  
  pmtx_->unlock();

  sa_size_t prm = nbfft_/10;
  for(sa_size_t i=0; i<nbfft_; i++) {
    sa_size_t first = rgen.Flat01()*(vv_.Size()-paqsz_-10);
    sa_size_t last = first+paqsz_-1;
    TVector<r_4> vx = vv_(Range(first,last));
    fftwf_execute(plan); 

    //    ffts_.FFTForward(vx, cfour_);
//    for(sa_size_t j=0; j<spectre_.Size(); j++) 
//      spectre_(j) += Zmod2(cfour_(j)); 
    if (i%prm == 0)  
      cout << "MTDoFFT::run() - Nom= " << nom_ << " Done i= " << i << endl;
  }
}
// ----- Fin de la definition de la classe MTDoFFT ----
//---------------------------------------------------------------------------------

//----------------- Les fonction de test -----------------

//------- f3_tmtfft()
int f3_tmtfft(sa_size_t VSZ, int NTH, sa_size_t nfft=10, sa_size_t paqsz=4096) 
{
  cout << "[1] f3_tmtfft/starting, VSZ= " << VSZ << " NTH= " << NTH << endl;
  cout << "  ... NbFFT=" << nfft << " PaqSz=" << paqsz << endl;
  ResourceUsage res(ResourceUsage::RU_All);
  Timer tm("f3_tmtfft");
  vector<MTDoFFT *> vth;
  FMTRandGen rg;
  TVector<r_4> DATA(VSZ);
  for(sa_size_t i=0; i<VSZ; i++) 
    DATA(i) = sin(i*M_PI*0.08654)+2*cos(i*M_PI*0.27591)+rg.Flatpm1();
  tm.Split("DATA-OK");
  cout << "[2] f3_tmtfft/creating threads " << endl;
  for(int kt=0; kt<NTH; kt++) {
    vth.push_back(new MTDoFFT(DATA, nfft, paqsz) );
  }
  tm.Split("Thread-Created");
  cout << "[3] f3_tmtfft/starting threads " << endl;
  for(int kt=0; kt<NTH; kt++)  vth[kt]->start();
  
  cout << "[4] f3_tmtfft/waiting for all  threads to finish " << endl;
  sleep(1);
  for(int kt=0; kt<NTH; kt++)  vth[kt]->join();
  tm.Split("Threads-Finished");

  cout << "[5] f3_tmtfft/saving spectra to file mtspectra.ppf " << endl;
  POutPersist po("mtspectra.ppf");
  for(int kt=0; kt<NTH; kt++)  
    po << vth[kt]->getSpectre();
  
  cout << "[6] f3_tmtfft/deleting thread objects " << endl;
  for(int kt=0; kt<NTH; kt++)  delete vth[kt];
  
  cout << res;
  return 0;
}
//----------------- FIN des fonctions de test -----------------

//-------------------------------------------------------------
//---------------------- MAIN MAIN MAIN -----------------------
//-------------------------------------------------------------

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

{
  InitTim();
  SophyaInit();

  int rc = 0;
  if (narg < 4) {
    cout << " tmtfft/Error args - Usage: tmtfft VSZ NTH NbFFT [PaqSz=4096]" << endl;
    return 1;
  }
  sa_size_t VSZ = atoi(arg[1]);
  int NTH = atoi(arg[2]);
  sa_size_t nfft = atoi(arg[3]);
  sa_size_t paqsz = 4096;
  if (narg > 4) paqsz = atoi(arg[4]);
  try {
    f3_tmtfft(VSZ, NTH, nfft, paqsz); 
  }
  catch (PThrowable exc) {
    cerr << "tmtfft: catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (...) {
    cerr << " catched unknown (...) exception (tmtfft.cc) " << endl; 
    rc = 78; 
  } 
  cout << "----------- tmtfft/END ------- " << endl;
  PrtTim("---END tmtfft---");
  return(rc);
  
}