From 40769be23822d8b4c39e1c8ddba0475bdcf725e7 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 25 Aug 2026 09:29:45 +0200 Subject: [PATCH] Implemented OpenMP parallelization of the independent-transform loops in both files, following the existing #ifdef HAVE_GOMP / #pragma omp parallel for pattern already used elsewhere in the codebase: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - src/musrFT.cpp: the timed Transform() loop over all PFourier objects is now parallelized. Added OpenMP::OpenMP_CXX linkage to the musrFT CMake target (it wasn't linked before). - src/classes/PMusrCanvas.cpp (HandleFourier() and HandleDifferenceFourier()): split each loop into three phases — serial PFourier construction (FFTW plan creation isn't thread-safe), a parallel Transform() loop (pure FFTW compute, safe to parallelize), and serial ROOT histogram creation/styling (ROOT globals aren't thread-safe). PMusr already links OpenMP, so only the #include and the pragma were needed. --- CMakeLists.txt | 2 +- src/CMakeLists.txt | 6 +- src/classes/PMusrCanvas.cpp | 126 ++++++++++++++++++++++++------------ src/musrFT.cpp | 11 +++- 4 files changed, 101 insertions(+), 44 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 058138f8..a4c8e868 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,7 @@ if (POLICY CMP0167) cmake_policy(SET CMP0167 NEW) endif () -project(musrfit VERSION 1.12.0 LANGUAGES C CXX) +project(musrfit VERSION 1.12.1 LANGUAGES C CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 51b75ec8..23ad2017 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -101,7 +101,11 @@ target_include_directories(musrFT $ $ ) -target_link_libraries(musrFT FFTW3::FFTW3 ${ROOT_LIBRARIES} ${MUSRFIT_LIBS} Boost::headers) +if (OpenMP_CXX_FOUND) + target_link_libraries(musrFT FFTW3::FFTW3 ${ROOT_LIBRARIES} ${MUSRFIT_LIBS} Boost::headers OpenMP::OpenMP_CXX) +else () + target_link_libraries(musrFT FFTW3::FFTW3 ${ROOT_LIBRARIES} ${MUSRFIT_LIBS} Boost::headers) +endif (OpenMP_CXX_FOUND) add_executable(musrRootValidation musrRootValidation.cpp) target_compile_options(musrRootValidation BEFORE PRIVATE "-DHAVE_CONFIG_H" "${HAVE_GIT_REV_H}") diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 80214a98..2db9b679 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -27,6 +27,8 @@ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ +#include + #include #include #include @@ -3374,24 +3376,60 @@ void PMusrCanvas::HandleFourier() bin = fHistoFrame->GetXaxis()->GetLast(); endTime = fHistoFrame->GetBinLowEdge(bin)+fHistoFrame->GetBinWidth(bin); } + // 1st: create the data/theory PFourier objects. FFTW plan creation is not + // thread-safe, hence this is done serially. + std::vector fourierData(fData.size(), nullptr); + std::vector fourierTheory(fData.size(), nullptr); for (UInt_t i=0; iIsValid()) { std::cerr << std::endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier data ..." << std::endl; + for (auto *ft : fourierData) delete ft; + for (auto *ft : fourierTheory) delete ft; return; } - fourierData.Transform(fFourier.fApodization); + + if (fTheoAsData) { // theory only at the data points + fourierTheory[i] = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower); + } else { + Int_t powerPad = fFourier.fFourierPower+5; // +5 means 8 times more points on theo (+3) + 4 times more points in fourier (+2) + fourierTheory[i] = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad); + } + if (!fourierTheory[i]->IsValid()) { + std::cerr << std::endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier theory ..." << std::endl; + for (auto *ft : fourierData) delete ft; + for (auto *ft : fourierTheory) delete ft; + return; + } + } + + // 2nd: execute the FFTs. Each data set is transformed independently, hence + // this can safely be parallelized. + Int_t noOfFT = static_cast(fData.size()); + #ifdef HAVE_GOMP + Int_t chunk = noOfFT/omp_get_num_procs(); + if (chunk < 1) + chunk = 1; + #pragma omp parallel for default(shared) schedule(dynamic,chunk) + #endif + for (Int_t i=0; iTransform(fFourier.fApodization); + fourierTheory[i]->Transform(fFourier.fApodization); + } + + // 3rd: extract the results into ROOT histograms and style them. ROOT + // object creation is not thread-safe, hence this is done serially. + for (UInt_t i=0; iGetBinWidth(1)/(endTime-startTime)); // get real part of the data - fData[i].dataFourierRe = fourierData.GetRealFourier(scale); + fData[i].dataFourierRe = fourierData[i]->GetRealFourier(scale); // get imaginary part of the data - fData[i].dataFourierIm = fourierData.GetImaginaryFourier(scale); + fData[i].dataFourierIm = fourierData[i]->GetImaginaryFourier(scale); // get power part of the data - fData[i].dataFourierPwr = fourierData.GetPowerFourier(scale); + fData[i].dataFourierPwr = fourierData[i]->GetPowerFourier(scale); // get phase part of the data - fData[i].dataFourierPhase = fourierData.GetPhaseFourier(); + fData[i].dataFourierPhase = fourierData[i]->GetPhaseFourier(); // set marker and line color fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); @@ -3415,42 +3453,25 @@ void PMusrCanvas::HandleFourier() fData[i].dataFourierPwr->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPhase->SetMarkerStyle(fData[i].data->GetMarkerStyle()); - // calculate fourier transform of the theory - Bool_t useFFTW = true; - PFourier *fourierTheory = nullptr; - if (fTheoAsData) { // theory only at the data points - fourierTheory = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower, useFFTW); - } else { - Int_t powerPad = fFourier.fFourierPower+5; // +5 means 8 times more points on theo (+3) + 4 times more points in fourier (+2) -#ifdef HAVE_DKS - if ((powerPad >= 20) && fUseDKS) - useFFTW = false; // i.e. use DKS -#endif - fourierTheory = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad, useFFTW); - } - if (!fourierTheory->IsValid()) { - std::cerr << std::endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier theory ..." << std::endl; - return; - } - fourierTheory->Transform(fFourier.fApodization); scale = sqrt(fData[0].theory->GetBinWidth(1)/(endTime-startTime)*fData[0].theory->GetBinWidth(1)/fData[0].data->GetBinWidth(1)); // get real part of the data - fData[i].theoryFourierRe = fourierTheory->GetRealFourier(scale); + fData[i].theoryFourierRe = fourierTheory[i]->GetRealFourier(scale); // get imaginary part of the data - fData[i].theoryFourierIm = fourierTheory->GetImaginaryFourier(scale); + fData[i].theoryFourierIm = fourierTheory[i]->GetImaginaryFourier(scale); // get power part of the data - fData[i].theoryFourierPwr = fourierTheory->GetPowerFourier(scale); + fData[i].theoryFourierPwr = fourierTheory[i]->GetPowerFourier(scale); // get phase part of the data - fData[i].theoryFourierPhase = fourierTheory->GetPhaseFourier(); - - // clean up - delete fourierTheory; + fData[i].theoryFourierPhase = fourierTheory[i]->GetPhaseFourier(); // set line colors for the theory fData[i].theoryFourierRe->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierIm->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPwr->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPhase->SetLineColor(fData[i].theory->GetLineColor()); + + // clean up + delete fourierData[i]; + delete fourierTheory[i]; } // phase opt. real FT requested initially in the msr-file, hence calculate it here @@ -3525,24 +3546,44 @@ void PMusrCanvas::HandleDifferenceFourier() bin = fHistoFrame->GetXaxis()->GetLast(); double endTime = fHistoFrame->GetBinCenter(bin); + // 1st: create the PFourier objects. FFTW plan creation is not + // thread-safe, hence this is done serially. + std::vector fourierData(fData.size(), nullptr); for (UInt_t i=0; iIsValid()) { std::cerr << std::endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier diff ..." << std::endl; + for (auto *ft : fourierData) delete ft; return; } - fourierData.Transform(fFourier.fApodization); + } + + // 2nd: execute the FFTs. Each data set is transformed independently, hence + // this can safely be parallelized. + Int_t noOfFT = static_cast(fData.size()); + #ifdef HAVE_GOMP + Int_t chunk = noOfFT/omp_get_num_procs(); + if (chunk < 1) + chunk = 1; + #pragma omp parallel for default(shared) schedule(dynamic,chunk) + #endif + for (Int_t i=0; iTransform(fFourier.fApodization); + } + + // 3rd: extract the results into ROOT histograms and style them. ROOT + // object creation is not thread-safe, hence this is done serially. + for (UInt_t i=0; iGetBinWidth(1)/(endTime-startTime)); // get real part of the data - fData[i].diffFourierRe = fourierData.GetRealFourier(scale); + fData[i].diffFourierRe = fourierData[i]->GetRealFourier(scale); // get imaginary part of the data - fData[i].diffFourierIm = fourierData.GetImaginaryFourier(scale); + fData[i].diffFourierIm = fourierData[i]->GetImaginaryFourier(scale); // get power part of the data - fData[i].diffFourierPwr = fourierData.GetPowerFourier(scale); + fData[i].diffFourierPwr = fourierData[i]->GetPowerFourier(scale); // get phase part of the data - fData[i].diffFourierPhase = fourierData.GetPhaseFourier(); + fData[i].diffFourierPhase = fourierData[i]->GetPhaseFourier(); // set marker and line color fData[i].diffFourierRe->SetMarkerColor(fData[i].diff->GetMarkerColor()); @@ -3567,6 +3608,9 @@ void PMusrCanvas::HandleDifferenceFourier() // set diffFourierTag fData[i].diffFourierTag = 1; // d-f + + // clean up + delete fourierData[i]; } // apply phase diff --git a/src/musrFT.cpp b/src/musrFT.cpp index b00c3ebd..a8bb997a 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -31,6 +31,8 @@ #include "config.h" #endif +#include + #include #include @@ -1463,7 +1465,14 @@ Int_t main(Int_t argc, Char_t *argv[]) apodTag = F_APODIZATION_STRONG; Double_t start = millitime(); - for (UInt_t i=0; i(fourier.size()); + #ifdef HAVE_GOMP + Int_t chunk = noOfFourier/omp_get_num_procs(); + if (chunk < 1) + chunk = 1; + #pragma omp parallel for default(shared) schedule(dynamic,chunk) + #endif + for (Int_t i=0; iTransform(apodTag); } Double_t end = millitime();