LMU/musrfit
5
1
Fork 0
Code Issues 2 Pull Requests Actions Packages Projects Releases 6 Activity

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:

Browse Source 
- 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 <omp.h> and the pragma were
    needed.
This commit is contained in:
suter_a
2026-08-25 10:00:41 +02:00
parent 66472c979a
commit 40769be238
4 changed files with 101 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+1 -1
View File
@@ -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)
src/CMakeLists.txt
+5 -1
View File
@@ -101,7 +101,11 @@ target_include_directories(musrFT
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/include>
)
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}")
src/classes/PMusrCanvas.cpp
+85 -41
View File
@@ -27,6 +27,8 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <omp.h>
#include <iostream>
#include <iomanip>
#include <fstream>
@@ -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<PFourier*> fourierData(fData.size(), nullptr);
std::vector<PFourier*> fourierTheory(fData.size(), nullptr);
for (UInt_t i=0; i<fData.size(); i++) {
// calculate fourier transform of the data
PFourier fourierData(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
fourierData[i] = new PFourier(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData[i]->IsValid()) {
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<Int_t>(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; i<noOfFT; i++) {
fourierData[i]->Transform(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; i<fData.size(); i++) {
double scale;
scale = sqrt(fData[0].data->GetBinWidth(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<PFourier*> fourierData(fData.size(), nullptr);
for (UInt_t i=0; i<fData.size(); i++) {
// calculate fourier transform of the data
PFourier fourierData(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
fourierData[i] = new PFourier(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData[i]->IsValid()) {
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<Int_t>(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; i<noOfFT; i++) {
fourierData[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; i<fData.size(); i++) {
double scale;
scale = sqrt(fData[0].diff->GetBinWidth(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
src/musrFT.cpp
+10 -1
View File
@@ -31,6 +31,8 @@
#include "config.h"
#endif
#include <omp.h>
#include <sys/time.h>
#include <iostream>
@@ -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(); i++) {
Int_t noOfFourier = static_cast<Int_t>(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; i<noOfFourier; i++) {
fourier[i]->Transform(apodTag);
}
Double_t end = millitime();