From 75578f1977ac3f2b0e904e80fc565a2e6c49ae62 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 7 Apr 2015 16:44:20 +0200 Subject: [PATCH] first work to add GPU support via DKS for Fourier. --- configure.ac | 44 ++- src/Makefile.am | 6 +- src/classes/Makefile.am | 4 +- src/classes/PFourier.cpp | 286 ++++++++++++++---- src/include/PFourier.h | 32 +- src/tests/DKS_FourierTest/Makefile | 97 ++++++ src/tests/DKS_FourierTest/dks_fourierTest.cpp | 143 +++++++++ 7 files changed, 545 insertions(+), 67 deletions(-) create mode 100644 src/tests/DKS_FourierTest/Makefile create mode 100644 src/tests/DKS_FourierTest/dks_fourierTest.cpp diff --git a/configure.ac b/configure.ac index d91c0836..d021dd81 100644 --- a/configure.ac +++ b/configure.ac @@ -845,6 +845,33 @@ if test "x$enable_omp" != "xno"; then [CXXFLAGS="$SAVED_CXXFLAGS" LIBS="$SAVED_LIBS"], []) fi +dnl ----------------------------------------------- +dnl Ask user if DKS support should be enabled (DKS = Dynamic Kernel Scheduler. Used to interface GPU's, MIC, etc) +dnl ----------------------------------------------- + +AC_ARG_ENABLE([dks], [AS_HELP_STRING([--enable-dks],[build musrfit with DKS (GPU/MIC) support [default=no]])], [DKS_ENABLED=1], [DKS_ENABLED=0]) + +DKS_LIBDIR="" +DKS_INCDIR="" +DKS_LIBS="" +DKS_CFLAGS="" +if test "${DKS_ENABLED}" == "1"; then + dnl --------------------------------------------- + dnl Eventually some strict testing about the availibilty of DKS should go in here + dnl Most probably it will be most sensible to ship DKS with the musrfit source code + dnl All hard wired stuff needs to disappear + dnl --------------------------------------------- + AC_DEFINE([HAVE_DKS], [1], [Define to 1 if DKS is available]) + AC_DEFINE([DKS_CUDA], [1], [Define to 1 if DKS Cuda is available]) + DKS_CUDADIR="/usr/local/cuda-6.5/targets/x86_64-linux" + DKS_INCDIR="/home/l_suter_a/DKS/src" + DKS_CFLAGS="-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" + AC_SUBST(DKS_CFLAGS) + AC_SUBST(DKS_LIBS) +fi + dnl ----------------------------------------------- dnl Ask user if the building of musredit/musrgui should be disabled dnl ----------------------------------------------- @@ -1126,6 +1153,7 @@ AM_CONDITIONAL([BUILD_CUBALIB], [test "${BUILD_CUBA}" = "1"]) AM_CONDITIONAL([BUILD_BMWLIBS], [test "${BUILD_BMW_LIBS}" = "1"]) AM_CONDITIONAL([BUILD_ASLIBS], [test "${BUILD_AS_LIBS}" = "1"]) AM_CONDITIONAL([BUILD_BNMRLIBS], [test "${BUILD_BNMR_LIBS}" = "1"]) +AM_CONDITIONAL([DKS_ENABLED], [test "${DKS_ENABLED}" = "1"]) AC_CONFIG_FILES([Makefile \ src/Makefile \ @@ -1133,8 +1161,8 @@ AC_CONFIG_FILES([Makefile \ src/classes/PMusr.pc \ src/classes/PUserFcnBase.pc \ src/external/Makefile \ - src/external/MusrRoot/Makefile \ - src/external/MusrRoot/TMusrRunHeader.pc \ + src/external/MusrRoot/Makefile \ + src/external/MusrRoot/TMusrRunHeader.pc \ src/external/TLemRunHeader/Makefile \ src/external/TLemRunHeader/TLemRunHeader.pc \ src/external/MuSR_software/Makefile \ @@ -1162,8 +1190,8 @@ AC_CONFIG_FILES([Makefile \ src/external/libCalcMeanFieldsLEM/Makefile \ src/external/Nonlocal/Makefile \ src/external/MagProximity/Makefile \ - src/external/libSpinValve/Makefile \ - src/external/libSpinValve/classes/Makefile \ + src/external/libSpinValve/Makefile \ + src/external/libSpinValve/classes/Makefile \ src/external/libPhotoMeissner/Makefile \ src/external/libPhotoMeissner/classes/Makefile \ src/external/libBNMR/Makefile \ @@ -1261,11 +1289,17 @@ else echo " Qt not needed (Qt editors disabled)" fi echo "" +if test "${DKS_ENABLED}" -eq 1; then + echo " DKS enabled" +else + echo " DKS disabled" +fi +echo "" echo "" echo " Features:" echo " ---------" echo "" -echo " musrfit (including musrfit, musrview, musrt0," +echo " musrfit (including musrfit, musrview, musrFT, musrt0," echo " msr2msr, msr2data, any2many, dump_header," echo " musrRootValidation, write_musrRoot_runHeader): yes" echo "" diff --git a/src/Makefile.am b/src/Makefile.am index 733f3445..8d6887ad 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -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)" diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index f6c13d5d..5e81e809 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -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 diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index 38f94d12..3bb48825 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -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[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(fNoOfBins, ierr); + fComp_ptr = fDks.allocateMemory< complex >(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(fReal_ptr, (int)fNoOfBins); + int size = fNoOfBins/2; + if (fNoOfBins % 2 == 1) + size++; + fDks.freeMemory< complex >(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(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 >(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; iSet 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; iSetBinContent(i+1, scale*fOut[i][0]); - realFourier->SetBinError(i+1, 0.0); + if (fUseFFTW) { + for (UInt_t i=0; iSetBinContent(i+1, scale*fOut[i][0]); + realFourier->SetBinError(i+1, 0.0); + } + } else { + for (UInt_t i=0; iSetBinContent(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; iSetBinContent(i+1, scale*fOut[i][1]); - imaginaryFourier->SetBinError(i+1, 0.0); + if (fUseFFTW) { + for (UInt_t i=0; iSetBinContent(i+1, scale*fOut[i][1]); + imaginaryFourier->SetBinError(i+1, 0.0); + } + } else { + for (UInt_t i=0; iSetBinContent(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; iSetBinContent(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; iSetBinContent(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; iSetBinContent(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= 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; iGetBinContent(i+start) - mean; - fIn[i][1] = 0.0; - } - for (UInt_t i=fNoOfData; iGetBinContent(i+start) - mean; + fIn[i][1] = 0.0; + } + for (UInt_t i=fNoOfData; iGetBinContent(i+start) - mean; +#endif + } + for (UInt_t i=fNoOfData; i +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 *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); diff --git a/src/tests/DKS_FourierTest/Makefile b/src/tests/DKS_FourierTest/Makefile new file mode 100644 index 00000000..23b12d1e --- /dev/null +++ b/src/tests/DKS_FourierTest/Makefile @@ -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 $< + diff --git a/src/tests/DKS_FourierTest/dks_fourierTest.cpp b/src/tests/DKS_FourierTest/dks_fourierTest.cpp new file mode 100644 index 00000000..2074b576 --- /dev/null +++ b/src/tests/DKS_FourierTest/dks_fourierTest.cpp @@ -0,0 +1,143 @@ +#include + +#include +#include +#include +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 data; + vector fourierData; + vector 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 fourierData.size()=" << fourierData.size() << endl; + + // do FFT + double start = millitime(); + for (unsigned int i=0; iTransform(); + } + 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; iGetPowerFourier(); + fourierPowerHistos.push_back(pf); + } + + cout << "debug> fourierPowerHistos.size()=" << fourierPowerHistos.size() << endl; + + // dump results + for (unsigned int i=0; iGetNbinsX(); + 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; iWrite(); + } + fout.Close(); + + // clean up + for (unsigned int i=0; i