LMU/musrfit
5
0
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity

first work to add GPU support via DKS for Fourier.

Browse Source
This commit is contained in:
suter_a
2015-04-07 16:44:20 +02:00
parent e4f11aca8c
commit 75578f1977
7 changed files with 545 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/Makefile.am
View File

@@ -54,18 +54,18 @@ dump_header_SOURCES = dump_header.cpp
xmldir = $(bindir)
xml_DATA = musrfit_startup.xml
LIBADD = $(PMUSR_LIBS) $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS)
LIBADD = $(PMUSR_LIBS) $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS) $(DKS_LIB)
AM_CXXFLAGS = $(LOCAL_BIN_CXXFLAGS)
AM_LDFLAGS = $(LOCAL_BIN_LDFLAGS)
AM_CPPFLAGS = $(MUSR_ROOT_CFLAGS) $(LEM_CFLAGS) $(MUD_CFLAGS) $(PSIBIN_CFLAGS) $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(BOOST_CFLAGS) $(ROOT_CFLAGS) \
$(LIBXML2_CFLAGS)
$(LIBXML2_CFLAGS) $(DKS_CFLAGS)
if PNEXUS_ENABLED
AM_CPPFLAGS += $(HDF5_CFLAGS) $(NEXUS_CFLAGS) $(PNEXUS_CXXFLAGS)
endif
LIBS = $(PMUSR_LIBS) $(USERFCN_LIBS) $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS) \
$(FFTW3_LIBS) $(GSL_LIBS) $(ROOT_LIBS) $(LIBXML2_LIBS)
$(FFTW3_LIBS) $(GSL_LIBS) $(DKS_CFLAGS) $(ROOT_LIBS) $(LIBXML2_LIBS)
install-xmlDATA: $(xml_DATA)
test -z "$(xmldir)" || $(mkdir_p) "$(DESTDIR)$(xmldir)"

4
src/classes/Makefile.am
View File

@@ -83,7 +83,7 @@ dict_cpp_sources_userFcn = \
include_HEADERS = $(h_sources) $(h_sources_userFcn)
noinst_HEADERS = $(h_linkdef) $(dict_h_sources) $(h_linkdef_userFcn) $(dict_h_sources_userFcn)
AM_CPPFLAGS = -I$(top_srcdir)/src/include $(MUSR_ROOT_CFLAGS) $(PSIBIN_CFLAGS) $(MUD_CFLAGS) $(LEM_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(BOOST_CFLAGS) -I$(ROOTINCDIR) $(PNEXUS_CXXFLAGS) $(NEXUS_CFLAGS)
AM_CPPFLAGS = -I$(top_srcdir)/src/include $(MUSR_ROOT_CFLAGS) $(PSIBIN_CFLAGS) $(MUD_CFLAGS) $(LEM_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(BOOST_CFLAGS) -I$(ROOTINCDIR) $(PNEXUS_CXXFLAGS) $(NEXUS_CFLAGS) $(DKS_CFLAGS)
AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS)
BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources) $(dict_cpp_sources_userFcn) $(dict_h_sources_userFcn)
@@ -100,7 +100,7 @@ libPUserFcnBase_la_LIBADD = $(ROOT_LIBS)
libPUserFcnBase_la_LDFLAGS = -version-info $(MUSR_LIBRARY_VERSION) -release $(MUSR_RELEASE) $(AM_LDFLAGS)
libPMusr_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources)
libPMusr_la_LIBADD = libPUserFcnBase.la $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS) $(FFTW3_LIBS) $(GSL_LIBS) $(ROOT_LIBS)
libPMusr_la_LIBADD = libPUserFcnBase.la $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS) $(FFTW3_LIBS) $(GSL_LIBS) $(DKS_LIBS) $(ROOT_LIBS)
libPMusr_la_LDFLAGS = -version-info $(MUSR_LIBRARY_VERSION) -release $(MUSR_RELEASE) $(AM_LDFLAGS)
pkgconfigdir = $(libdir)/pkgconfig

286
src/classes/PFourier.cpp
View File

@@ -37,8 +37,6 @@ using namespace std;
#include "TF1.h"
#include "TAxis.h"
//#include "TFile.h"
#include "PMusr.h"
#include "PFourier.h"
@@ -59,7 +57,7 @@ using namespace std;
* \param dcCorrected if true, removed DC offset from signal before Fourier transformation, otherwise not
* \param zeroPaddingPower if set to values > 0, there will be zero padding up to 2^zeroPaddingPower
*/
PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, Bool_t dcCorrected, UInt_t zeroPaddingPower) :
PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, Bool_t dcCorrected, UInt_t zeroPaddingPower, Bool_t useFFTW) :
fData(data), fUnitTag(unitTag), fStartTime(startTime), fEndTime(endTime),
fDCCorrected(dcCorrected), fZeroPaddingPower(zeroPaddingPower)
{
@@ -71,8 +69,15 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
}
fValid = true;
fUseFFTW = true;
fIn = 0;
fOut = 0;
#ifdef HAVE_DKS
fInDKS = 0;
fOutDKS = 0;
#endif
SetUseFFTW(useFFTW);
fApodization = F_APODIZATION_NONE;
@@ -81,7 +86,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
// if endTime == 0 set it to the last time slot
if (fEndTime == 0.0) {
Int_t last = fData->GetNbinsX()-1;
Int_t last = fData->GetNbinsX();
fEndTime = fData->GetBinCenter(last);
}
@@ -128,21 +133,67 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
}
// allocate necessary memory
fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
// check if memory allocation has been successful
if ((fIn == 0) || (fOut == 0)) {
fValid = false;
return;
if (fUseFFTW) {
fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
} else { // try DKS
#ifdef HAVE_DKS
fInDKS = new double[fNoOfBins];
unsigned int size=fNoOfBins/2;
if (fNoOfBins % 2 == 1)
size++;
fOutDKS = new complex<double>[size];
#endif
}
// get the FFTW3 plan (see FFTW3 manual)
fFFTwPlan = fftw_plan_dft_1d(fNoOfBins, fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE);
// check if a valid plan has been generated
if (!fFFTwPlan) {
// check if memory allocation has been successful
if (fUseFFTW) {
if ((fIn == 0) || (fOut == 0)) {
fValid = false;
return;
}
} else { // try DKS
#ifdef HAVE_DKS
if ((fInDKS == 0) || (fOutDKS == 0)) {
fValid = false;
return;
}
#else
fValid = false;
return;
#endif
}
if (fUseFFTW) {
// get the FFTW3 plan (see FFTW3 manual)
fFFTwPlan = fftw_plan_dft_1d(fNoOfBins, fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE);
// check if a valid plan has been generated
if (!fFFTwPlan) {
fValid = false;
}
} else { // try DKS
#ifdef HAVE_DKS
// init DKSBase
fDks.setAPI("Cuda", 4);
fDks.setDevice("-gpu", 4);
fDks.initDevice();
int dimsize[3] = {fNoOfBins, 1, 1};
fDks.setupFFT(1, dimsize);
// allocate memory on accelerator
int ierr;
unsigned int size=fNoOfBins/2;
if (fNoOfBins % 2 == 1)
size++;
fReal_ptr = 0;
fComp_ptr = 0;
fReal_ptr = fDks.allocateMemory<double>(fNoOfBins, ierr);
fComp_ptr = fDks.allocateMemory< complex<double> >(size, ierr);
if ((fReal_ptr==0) || (fComp_ptr==0))
fValid = false;
#else
fValid = false;
#endif
}
}
@@ -154,12 +205,27 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
*/
PFourier::~PFourier()
{
if (fFFTwPlan)
fftw_destroy_plan(fFFTwPlan);
if (fIn)
fftw_free(fIn);
if (fOut)
fftw_free(fOut);
if (fUseFFTW) {
if (fFFTwPlan)
fftw_destroy_plan(fFFTwPlan);
if (fIn)
fftw_free(fIn);
if (fOut)
fftw_free(fOut);
} else {
#ifdef HAVE_DKS
// free accelerator memory
fDks.freeMemory<double>(fReal_ptr, (int)fNoOfBins);
int size = fNoOfBins/2;
if (fNoOfBins % 2 == 1)
size++;
fDks.freeMemory< complex<double> >(fComp_ptr, size);
if (fIn)
delete [] fInDKS;
if (fOut)
delete [] fOutDKS;
#endif
}
}
//--------------------------------------------------------------------------
@@ -178,7 +244,28 @@ void PFourier::Transform(UInt_t apodizationTag)
PrepareFFTwInputData(apodizationTag);
fftw_execute(fFFTwPlan);
if (fUseFFTW) {
fftw_execute(fFFTwPlan);
} else {
#ifdef HAVE_DKS
int dimsize[3] = {fNoOfBins, 1, 1};
int status=0, size=fNoOfBins/2;
if (fNoOfBins % 2 == 1)
size++;
// write data to the accelerator
status=fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins);
cout << "debug> status=" << status << ": write" << endl;
// execute the FFT
status = fDks.callR2CFFT(fReal_ptr, fComp_ptr, 1, dimsize);
cout << "debug> status=" << status << ": FFT" << endl;
// read data from accelerator
status = fDks.readData< complex<double> >(fComp_ptr, fOutDKS, size);
cout << "debug> status=" << status << ": read" << endl;
#else
fValid = false;
return;
#endif
}
// correct the phase for tstart != 0.0
// find the first bin >= fStartTime
@@ -191,12 +278,42 @@ void PFourier::Transform(UInt_t apodizationTag)
}
Double_t phase, re, im;
for (UInt_t i=0; i<fNoOfBins; i++) {
phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
re = fOut[i][0] * cos(phase) + fOut[i][1] * sin(phase);
im = -fOut[i][0] * sin(phase) + fOut[i][1] * cos(phase);
fOut[i][0] = re;
fOut[i][1] = im;
if (fUseFFTW) {
for (UInt_t i=0; i<fNoOfBins; i++) {
phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
re = fOut[i][0] * cos(phase) + fOut[i][1] * sin(phase);
im = -fOut[i][0] * sin(phase) + fOut[i][1] * cos(phase);
fOut[i][0] = re;
fOut[i][1] = im;
}
} else { // try DKS
for (UInt_t i=0; i<fNoOfBins; i++) {
phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
#ifdef HAVE_DKS
re = fOutDKS[i].real() * cos(phase) + fOutDKS[i].imag() * sin(phase);
im = -fOutDKS[i].real() * sin(phase) + fOutDKS[i].imag() * cos(phase);
fOutDKS[i].real() = re;
fOutDKS[i].imag() = im;
#endif
}
}
}
//--------------------------------------------------------------------------
// SetUseFFTW
//--------------------------------------------------------------------------
/**
* <p>Set fUseFFTW flag. Checks if DKS support is in place before setting the flag.
*/
void PFourier::SetUseFFTW(const Bool_t flag)
{
if (flag == false) {
#ifndef HAVE_DKS
fUseFFTW = true;
cerr << endl << "PFouier::SetUseFFTW: **ERROR** DKS not in use, will fall back to FFTW" << endl << endl;
#else
fUseFFTW = flag;
#endif
}
}
@@ -251,9 +368,20 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
}
// fill realFourier vector
for (UInt_t i=0; i<noOfFourierBins; i++) {
realFourier->SetBinContent(i+1, scale*fOut[i][0]);
realFourier->SetBinError(i+1, 0.0);
if (fUseFFTW) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
realFourier->SetBinContent(i+1, scale*fOut[i][0]);
realFourier->SetBinError(i+1, 0.0);
}
} else {
for (UInt_t i=0; i<noOfFourierBins; i++) {
#ifdef HAVE_DKS
realFourier->SetBinContent(i+1, scale*fOutDKS[i].real());
#else
realFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif
realFourier->SetBinError(i+1, 0.0);
}
}
return realFourier;
}
@@ -292,11 +420,21 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
}
// fill imaginaryFourier vector
for (UInt_t i=0; i<noOfFourierBins; i++) {
imaginaryFourier->SetBinContent(i+1, scale*fOut[i][1]);
imaginaryFourier->SetBinError(i+1, 0.0);
if (fUseFFTW) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
imaginaryFourier->SetBinContent(i+1, scale*fOut[i][1]);
imaginaryFourier->SetBinError(i+1, 0.0);
}
} else {
for (UInt_t i=0; i<noOfFourierBins; i++) {
#ifdef HAVE_DKS
imaginaryFourier->SetBinContent(i+1, scale*fOutDKS[i].imag());
#else
imaginaryFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif
imaginaryFourier->SetBinError(i+1, 0.0);
}
}
return imaginaryFourier;
}
@@ -334,11 +472,21 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
}
// fill powerFourier vector
for (UInt_t i=0; i<noOfFourierBins; i++) {
pwrFourier->SetBinContent(i+1, scale*sqrt(fOut[i][0]*fOut[i][0]+fOut[i][1]*fOut[i][1]));
pwrFourier->SetBinError(i+1, 0.0);
if (fUseFFTW) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
pwrFourier->SetBinContent(i+1, scale*sqrt(fOut[i][0]*fOut[i][0]+fOut[i][1]*fOut[i][1]));
pwrFourier->SetBinError(i+1, 0.0);
}
} else {
for (UInt_t i=0; i<noOfFourierBins; i++) {
#ifdef HAVE_DKS
pwrFourier->SetBinContent(i+1, scale*sqrt(fOutDKS[i].real()*fOutDKS[i].real()+fOutDKS[i].imag()*fOutDKS[i].imag()));
#else
pwrFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif
pwrFourier->SetBinError(i+1, 0.0);
}
}
return pwrFourier;
}
@@ -377,18 +525,31 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
// fill phaseFourier vector
Double_t value = 0.0;
Double_t re, im;
for (UInt_t i=0; i<noOfFourierBins; i++) {
if (fUseFFTW) {
re = fOut[i][0];
im = fOut[i][1];
} else {
#ifdef HAVE_DKS
re = fOutDKS[i].real();
im = fOutDKS[i].imag();
#else
re = 1.0;
im = 0.0;
#endif
}
// calculate the phase
if (fOut[i][0] == 0) {
if (fOut[i][1] >= 0.0)
if (re == 0) {
if (im >= 0.0)
value = PI_HALF;
else
value = -PI_HALF;
} else {
value = atan(fOut[i][1]/fOut[i][0]);
value = atan(re/im);
// check sector
if (fOut[i][0] < 0.0) {
if (fOut[i][1] > 0.0)
if (re < 0.0) {
if (im > 0.0)
value = PI + value;
else
value = PI - value;
@@ -437,13 +598,26 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
}
// 2nd fill fIn
for (UInt_t i=0; i<fNoOfData; i++) {
fIn[i][0] = fData->GetBinContent(i+start) - mean;
fIn[i][1] = 0.0;
}
for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
fIn[i][0] = 0.0;
fIn[i][1] = 0.0;
if (fUseFFTW) {
for (UInt_t i=0; i<fNoOfData; i++) {
fIn[i][0] = fData->GetBinContent(i+start) - mean;
fIn[i][1] = 0.0;
}
for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
fIn[i][0] = 0.0;
fIn[i][1] = 0.0;
}
} else {
for (UInt_t i=0; i<fNoOfData; i++) {
#ifdef HAVE_DKS
fInDKS[i] = fData->GetBinContent(i+start) - mean;
#endif
}
for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
#ifdef HAVE_DKS
fInDKS[i] = 0.0;
#endif
}
}
// 3rd apodize data (if wished)
@@ -503,6 +677,12 @@ void PFourier::ApodizeData(Int_t apodizationTag) {
for (UInt_t j=1; j<5; j++) {
q += c[j] * pow((Double_t)i/(Double_t)fNoOfData, 2.0*(Double_t)j);
}
fIn[i][0] *= q;
if (fUseFFTW) {
fIn[i][0] *= q;
} else {
#ifdef HAVE_DKS
fInDKS[i] *= q;
#endif
}
}
}

32
src/include/PFourier.h
View File

@@ -30,7 +30,17 @@
#ifndef _PFOURIER_H_
#define _PFOURIER_H_
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifndef HAVE_DKS
#include "fftw3.h"
#else
#include <complex>
using namespace std;
#include "DKSBase.h"
#endif
#include "PMusr.h"
@@ -47,10 +57,13 @@ class PFourier
public:
PFourier(TH1F *data, Int_t unitTag,
Double_t startTime = 0.0, Double_t endTime = 0.0,
Bool_t dcCorrected = false, UInt_t zeroPaddingPower = 0);
Bool_t dcCorrected = false, UInt_t zeroPaddingPower = 0,
Bool_t useFFTW = true);
virtual ~PFourier();
virtual void Transform(UInt_t apodizationTag = 0);
virtual void Transform(UInt_t apodizationTag = F_APODIZATION_NONE);
virtual void SetUseFFTW(const Bool_t flag);
virtual const char* GetDataTitle() { return fData->GetTitle(); }
virtual const Int_t GetUnitTag() { return fUnitTag; }
@@ -62,12 +75,14 @@ class PFourier
virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0);
virtual Bool_t IsValid() { return fValid; }
virtual Bool_t IsUseFFTW() { return fUseFFTW; }
private:
TH1F *fData; ///< data histogram to be Fourier transformed.
Bool_t fValid; ///< true = all boundary conditions fullfilled and hence a Fourier transform can be performed.
Int_t fUnitTag; ///< 1=Field Units (G), 2=Field Units (T), 3=Frequency Units (MHz), 4=Angular Frequency Units (Mc/s)
Bool_t fUseFFTW; ///< true = use FFTW, otherwise use DKS if present
Int_t fApodization; ///< 0=none, 1=weak, 2=medium, 3=strong
@@ -80,9 +95,18 @@ class PFourier
UInt_t fNoOfData; ///< number of bins in the time interval between fStartTime and fStopTime
UInt_t fNoOfBins; ///< number of bins to be Fourier transformed. Might be different to fNoOfData due to zero padding
fftw_plan fFFTwPlan; ///< fftw plan (see FFTW3 User Manual)
fftw_complex *fIn; ///< real part of the Fourier transform
fftw_complex *fOut; ///< imaginary part of the Fourier transform
fftw_complex *fIn; ///< real part of the Fourier transform
fftw_complex *fOut; ///< imaginary part of the Fourier transform
#ifdef HAVE_DKS
double *fInDKS; ///< real part of the Fourier transform
complex<double> *fOutDKS; ///< imaginary part of the Fourier transform
DKSBase fDks; ///< Dynamic Kernel Scheduler
void *fReal_ptr; ///< real part of the Fourier on accelartor
void *fComp_ptr; ///< imaginary part of the Fourier on the acclerator
#endif
virtual void PrepareFFTwInputData(UInt_t apodizationTag);
virtual void ApodizeData(Int_t apodizationTag);

97
src/tests/DKS_FourierTest/Makefile Normal file
View File

@@ -0,0 +1,97 @@
#---------------------------------------------------
# get compilation and library flags from root-config
ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags)
ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs)
ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs)
#---------------------------------------------------
# set compilation and library flags for DKS
DKS_CUDADIR = /usr/local/cuda-6.5/targets/x86_64-linux
DKS_INCDIR = /home/l_suter_a/DKS/src
DKS_CFLAGS = -DHAVE_DKS -DDKS_OPENCL -DDKS_CUDA -I${DKS_INCDIR} -I${DKS_CUDADIR}/include
DKS_LIBDIR = /home/l_suter_a/DKS/build/src
DKS_LIBS = -L${DKS_CUDADIR}/lib -lOpenCL -L${DKS_LIBDIR} -ldksshared
#---------------------------------------------------
# set compilation and library flags for FFTW
FFTW3_INCDIR = /usr/include
FFTW3_LIBS = -L/usr/lib64 -lfftw3
#---------------------------------------------------
# set compilation and library flags for PMusr
PMUSR_INCDIR = $(ROOTSYS)/include
PMUSR_SRCDIR = $(HOME)/musrfit/src/classes
PMUSR_LIBS = -L $(ROOTSYS)/lib -lPMusr
#---------------------------------------------------
# depending on the architecture, choose the compiler,
# linker, and the flags to use
#
ARCH = $(shell $(ROOTSYS)/bin/root-config --arch)
ifeq ($(ARCH),linux)
OS = LINUX
endif
ifeq ($(ARCH),linuxx8664gcc)
OS = LINUX
endif
ifeq ($(ARCH),win32gcc)
OS = WIN32GCC
endif
ifeq ($(ARCH),macosx)
OS = DARWIN
endif
# -- Linux
ifeq ($(OS),LINUX)
CXX = gcc
CXXFLAGS = -g -O2 -Wall -fPIC
INCLUDES = -I$(PMUSR_INCDIR) -I$(DKS_INCDIR)
LD = g++
LDFLAGS =
INSTALLPATH = ./
EXEC = dks_fourierTest
SUFFIX =
endif
# -- MacOSX/Darwin
ifeq ($(OS),DARWIN)
CXX = gcc
CXXFLAGS = -O3 -Wall -fPIC
INCLUDES = -I$(PMUSR_INCDIR) -I$(FFTW3_INCDIR) -I$(DKS_INCDIR)
LD = g++
LDFLAGS = -O
INSTALLPATH = ./
EXEC = dks_fourierTest
SUFFIX =
endif
# the output from the root-config script:
CXXFLAGS += $(ROOTCFLAGS) $(DKS_CFLAGS)
LDFLAGS +=
# the ROOT libraries (G = graphic)
LIBS = $(ROOTLIBS) -lXMLParser $(FFTW3_LIBS) $(DKS_LIBS)
GLIBS = $(ROOTGLIBS) -lXMLParser $(FFTW3_LIBS) $(DKS_LIBS)
# some definitions: headers, sources, objects,...
FOURIER = PFourier
all: $(EXEC)
clean:
rm *.o
$(EXEC): $(FOURIER).o $(EXEC).o
@echo "---> Building $(EXEC) ..."
/bin/rm -f $(SHLIB)
$(LD) $(FOURIER).o $(EXEC).o -o $(EXEC) $(GLIBS)
@echo "done"
$(FOURIER).o: $(PMUSR_SRCDIR)/$(FOURIER).cpp
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<

143
src/tests/DKS_FourierTest/dks_fourierTest.cpp Normal file
View File

@@ -0,0 +1,143 @@
#include <sys/time.h>
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;
#include "TH1F.h"
#include "TFile.h"
#include "PMusr.h"
#include "PFourier.h"
//---------------------------------------------------
double millitime()
{
struct timeval now;
gettimeofday(&now, 0);
return ((double)now.tv_sec * 1.0e6 + (double)now.tv_usec)/1.0e3;
}
//---------------------------------------------------
void dks_fourierTest_syntax()
{
cout << endl << "usage: dks_fourierTest [useFFTW, N, L | --help] ";
cout << endl << " useFFTW : flag, if true -> FFTW, otherwise -> DKS";
cout << endl << " N : number of histos";
cout << endl << " L : histo length";
cout << endl << " if not given, useFFTW=true, N=8, L=2^20=1048576";
cout << endl << " --help: this help";
cout << endl << endl;
}
//---------------------------------------------------
int main(int argc, char *argv[])
{
vector<TH1F*> data;
vector<PFourier*> fourierData;
vector<TH1F*> fourierPowerHistos;
int N=8, L=1048576;
Bool_t useFFTW = true;
switch (argc) {
case 1:
break;
case 4:
if (!strcmp(argv[1], "true") || !strcmp(argv[1], "1"))
useFFTW = true;
else if (!strcmp(argv[1], "false") || !strcmp(argv[1], "0"))
useFFTW = false;
else {
dks_fourierTest_syntax();
return 1;
}
N = strtol(argv[2], 0, 10);
L = strtol(argv[3], 0, 10);
break;
default:
dks_fourierTest_syntax();
return 1;
break;
}
if ((N<=0) || (L<=0)) {
dks_fourierTest_syntax();
return 1;
}
// feed the histos
TH1F *h;
char str[128];
double dt = 10.0 / (double)L;
double dval = 0.0;
for (unsigned int i=0; i<(unsigned int)N; i++) {
snprintf(str, sizeof(str), "h%d", i);
h = new TH1F(str, str, L+1, -dt/2, 10.0+dt/2);
for (unsigned int j=0; j<(unsigned int)h->GetNbinsX(); j++) {
dval = exp(-0.5*dt*j)*cos(0.013554*5000.0*dt*j+6.2832/N*i);
h->SetBinContent(j+1, dval);
}
data.push_back(h);
}
cout << "debug> data.size()=" << data.size() << endl;
// setup FFT
PFourier *fh;
for (unsigned int i=0; i<data.size(); i++) {
fh = new PFourier(data[i], FOURIER_UNIT_GAUSS, 0.0, 0.0, false, 0, useFFTW);
fourierData.push_back(fh);
}
cout << "debug> fourierData.size()=" << fourierData.size() << endl;
// do FFT
double start = millitime();
for (unsigned int i=0; i<fourierData.size(); i++) {
fourierData[i]->Transform();
}
double end = millitime();
cout << "info> overall computation time: " << (end-start)/1.0e3 << " (sec)." << endl << endl;
// get FFT data
TH1F *pf;
unsigned int out=128;
for (unsigned int i=0; i<fourierData.size(); i++) {
pf = fourierData[i]->GetPowerFourier();
fourierPowerHistos.push_back(pf);
}
cout << "debug> fourierPowerHistos.size()=" << fourierPowerHistos.size() << endl;
// dump results
for (unsigned int i=0; i<fourierPowerHistos.size(); i++) {
out = (unsigned int)fourierPowerHistos[i]->GetNbinsX();
cout << "debug> out=" << out << endl;
if (out > 64)
out = 64;
cout << "debug> fourier " << i+1 << ": ";
for (unsigned int j=1; j<=out; j++)
cout << fourierPowerHistos[i]->GetBinContent(j) << ", ";
cout << endl << "---" << endl;
}
TFile fout("dks_fourierTest.root", "recreate");
for (unsigned int i=0; i<fourierPowerHistos.size(); i++) {
fourierPowerHistos[i]->Write();
}
fout.Close();
// clean up
for (unsigned int i=0; i<data.size(); i++) {
delete data[i];
delete fourierData[i];
delete fourierPowerHistos[i];
}
data.clear();
fourierData.clear();
fourierPowerHistos.clear();
return 0;
}