BenchmarkFFT.C

Go to the documentation of this file.

#ifndef _HAVE_FFTW3_
#error Program requires FFTW3 headers to run. Check your environment settings and that DIANA_HAVE_FFTW3 is set
#endif
 
#include <iostream>
#include "QVector.hh"
#include "QVectorC.hh"
#include "QRealComplexFFT.hh"
#include "QRealComplexFFTW3.hh"
#include "QRealComplexFFTGSL.hh"
#include "QFFTWPlan.hh"
#include <TRandom.h>
#include <TTree.h>
#include <TFile.h>
#include <vector>
#include <getopt.h>
 
const char *progname;
 
 
struct BenchmarkData {
  int fNGSLTests;
  int fNFFTWTests;
  int fNFFTWPlanTests;
  int fWindowSize;
  double fGSLPrepTime;
  double fFFTWPrepTime;
  double fGSLForwardTimePerPulse;
  double fFFTWForwardTimePerPulse;
  double fGSLBackwardTimePerPulse;
  double fFFTWBackwardTimePerPulse;
  double fFFTWForwardPlanTimePerPulse;
  double fFFTWBackwardPlanTimePerPulse;
  
};
 
void usage() {
  std::cout << progname << std::endl;
  std::cout << "Benchmark FFT algorithms for different size windows.\n" << std::endl;
  std::cout << "Options: " << std::endl;
  std::cout << "  -n  :  Number of runs per test." << std::endl;
  std::cout << "  -w  :  Add a window size to test. (Can specify multiple times)." << std::endl;
  std::cout << "  -o  :  Output to a root file." << std::endl;
  std::cout << "  -h  :  This helpful message." << std::endl;
}
 
 
double GetTime_ms(time_t t0,clock_t c0,time_t t1,clock_t c1){
  clock_t cDiff=c1-c0;
  double tDiff=difftime(t1,t0);
  int wrappedTime=((int)(tDiff/4294))*4294;
  return (1.E3*cDiff/CLOCKS_PER_SEC)+1.E3*wrappedTime;
}
 
 
void BenchmarkTransform(int NRuns,int WinSize,BenchmarkData &results) {
 
  std::cout << "Benchmarking window size " << WinSize << " ... " << std::endl;
  std::cout << "Running " << NRuns << " transforms per test ..."  << std::endl;
 
  results.fNGSLTests = NRuns;
  results.fNFFTWTests = NRuns;
  results.fWindowSize = WinSize;
 
  bool compress = false;
  Diana::QVector X1(WinSize),X2(WinSize),X3(WinSize);
  Diana::QVectorC Y1(compress ? WinSize/2+1 : WinSize),
    Y2(compress ? WinSize/2+1 : WinSize),
    Y3(WinSize/2+1);
  for(int j=0;j<WinSize;j++)
    X1[j]=X2[j]=X3[j]=gRandom->Gaus();
 
  clock_t c1 = ::clock();
  time_t t1 = ::time(NULL);
  Diana::QRealComplexFFTW3 *FFTW = new Diana::QRealComplexFFTW3(WinSize);
  FFTW->TransformToFreq(X1,Y1);
  FFTW->TransformFromFreq(Y1,X1);
  clock_t c2 = ::clock();
  time_t t2 = ::time(NULL);
  results.fFFTWPrepTime=1000.*GetTime_ms(t1,c1,t2,c2);
  
  c1 = ::clock();
  t1 = ::time(NULL);
  Diana::QRealComplexFFTGSL *GSL = new Diana::QRealComplexFFTGSL(WinSize);
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fGSLPrepTime=GetTime_ms(t1,c1,t2,c2);
 
  Diana::QFFTWRealComplexPlan FFTWForwardPlan = Diana::QFFTWRealComplexPlan::CreateForwardPlan(X3,Y3);
  Diana::QFFTWRealComplexPlan FFTWBackwardPlan = Diana::QFFTWRealComplexPlan::CreateBackwardPlan(Y3,X3);
  
  std::cout << "Running " << NRuns << " Forward Tests With GSL ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    GSL->TransformToFreq(X2,Y2,compress);
  c2 = ::clock();
  t1 = ::time(NULL);
  results.fGSLForwardTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
  
  std::cout << "Running " << NRuns << " Forward Tests With FFTW ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    FFTW->TransformToFreq(X1,Y1,compress);
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fFFTWForwardTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
 
  std::cout << "Running " << NRuns << " Forward Tests With FFTW Plan ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    FFTWForwardPlan.Execute();
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fFFTWForwardPlanTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
 
  std::cout << "Running " << NRuns << " Backward Tests With GSL ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    GSL->TransformFromFreq(Y2,X2,compress);
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fGSLBackwardTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
 
  std::cout << "RUnning " << NRuns << " Backward Tests With FFTW ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    FFTW->TransformFromFreq(Y1,X1,compress);
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fFFTWBackwardTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
 
  std::cout << "Running " << NRuns << " Backward Tests With FFTW Plan ..." << std::endl;
  c1 = ::clock();
  t1 = ::time(NULL);
  for(int j=0;j<NRuns;j++)
    FFTWBackwardPlan.Execute();
  c2 = ::clock();
  t2 = ::time(NULL);
  results.fFFTWBackwardPlanTimePerPulse = 1000.*GetTime_ms(t1,c1,t2,c2)/NRuns;
 
  
  std::cout << "GSL Prep:\t" << results.fGSLPrepTime/1000 << "ms" << std::endl;
  std::cout << "FFTW Prep:\t" << results.fFFTWPrepTime/1000 << "ms" << std::endl;
  
  std::cout << Form("GSL Forward:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fGSLForwardTimePerPulse/1000,results.fGSLForwardTimePerPulse)
            << std::endl;
  std::cout << Form("FFTW Forward:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fFFTWForwardTimePerPulse/1000,results.fFFTWForwardTimePerPulse)
            << std::endl;
  std::cout << Form("FFTW Forward Plan:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fFFTWForwardPlanTimePerPulse/1000,results.fFFTWForwardPlanTimePerPulse)
            << std::endl;
  Diana::QVectorC Diff=Y1-Y2;
  std::cout << "MeanSqDif: " 
            << Diff.InterleavedVector().GetRMS(2*Diff.Size()) << std::endl;
  
  std::cout << Form("GSL Backward:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fGSLBackwardTimePerPulse/1000,results.fGSLBackwardTimePerPulse)
            << std::endl;
  std::cout << Form("FFTW Backward:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fFFTWBackwardTimePerPulse/1000,results.fFFTWBackwardTimePerPulse)
            << std::endl;
  std::cout << Form("FFTW Backward Plan:\t%.2f ms (%4.2f us/pulse)",
                    NRuns*results.fFFTWBackwardPlanTimePerPulse/1000,results.fFFTWBackwardPlanTimePerPulse)
            << std::endl;
  
  Diana::QVector Diff2=X1-X2;
  std::cout << "MeanSqDif: " << Diff2.GetRMS(Diff2.Size()) << std::endl;
  
  std::cout<<"GSL Total:\t"<< results.fGSLPrepTime/1000 + NRuns * (results.fGSLForwardTimePerPulse + results.fGSLBackwardTimePerPulse)/1000
           << "ms" <<std::endl;
  std::cout<<"FFTW Total:\t"<< results.fFFTWPrepTime/1000 + NRuns * (results.fFFTWForwardTimePerPulse + results.fFFTWBackwardTimePerPulse)/1000
           << "ms" <<std::endl;
  std::cout << std::endl;
 
 
 
  
 
  
  
  delete GSL;delete FFTW;
 
}
  
int main(int argc,char *argv[]){
 
  progname=(progname=strrchr(argv[0],'/'))==NULL ? argv[0] : progname+1;
 
  int NRuns = 10000;
  std::vector<int> windowSizes;
  TFile *outFile = NULL;
  TTree *tree = NULL;
 
  int c;
  while((c = getopt(argc,argv,":n:w:o:h"))!=-1) {
    switch(c) {
    case 'n':
      NRuns = atoi(optarg);
      break;
    case 'w':
      windowSizes.push_back(atoi(optarg));
      break;
    case 'o':
      if(!tree) {
        outFile = new TFile(optarg,"RECREATE");
        tree = new TTree("FFTBenchmark","FFT Benchmark Results");
      }
      break;
    case 'h':
      usage();
      return 0;
    case ':':
      std::cout << "Option -" << (char)optopt << " requires an argument." << std::endl;
      return 2;
    case '?':
      std::cout << "Unknown option: -" << (char)optopt << std::endl;
      return 2;
    }
  }
 
  std::cout << "Testing " << NRuns << " FFTs for window sizes: ";
  for(size_t i=0;i<windowSizes.size();i++)
    std::cout << windowSizes[i] << " ";
  std::cout << std::endl;
 
  BenchmarkData benchmarkData;
  if(tree){
    tree->Branch("NGSLTests",&NRuns,"NGSLTests/I");
    tree->Branch("NFFTWTests",&NRuns,"NFFTWTests/I");
    tree->Branch("WindowSize",&benchmarkData.fWindowSize,"WindowSize/I");
    tree->Branch("GSLForwardTimePerPulse_us",&benchmarkData.fGSLForwardTimePerPulse,"GSLForwardTimePerPulse_us/D");
    tree->Branch("GSLBackwardTimePerPulse_us",&benchmarkData.fGSLBackwardTimePerPulse,"GSLBackwardTimePerPulse_us/D");
    tree->Branch("FFTWForwardTimePerPulse_us",&benchmarkData.fFFTWForwardTimePerPulse,"FFTWForwardTimePerPulse_us/D");
    tree->Branch("FFTWBackwardTimePerPulse_us",&benchmarkData.fFFTWBackwardTimePerPulse,"FFTWBackwardTimePerPulse_us/D");
  }
 
  for(size_t i=0;i<windowSizes.size();i++){
    BenchmarkTransform(NRuns,windowSizes[i],benchmarkData);
    if(tree) {
      tree->Fill();
      tree->Write("",TTree::kOverwrite);
    }
  }
 
  if(tree) {
    tree->Write();
    outFile->Close();
  }
 
  return 0;
}