From 7dea9eaf9c88d034218302e6aa89eccf7e8816d9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 19 Sep 2026 16:04:23 +0200 Subject: [PATCH] One process-wide lock for FFTW planning FFTW's planner (every fftwf_plan_* and fftwf_destroy_plan) shares global state and is not thread-safe; executing a plan is. FFTIndexerCPU and BeamCenterFFTCPU each guarded their planning with a lock of their own, TranslationalNCS and the viewer's spectrum with none, and ModelFFT with a third - which does not stop two of them planning at once. common/FFTWPlannerLock.h holds the one mutex they all now take. No numerical change. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- common/FFTWPlannerLock.h | 14 +++++++++++++ .../geom_refinement/BeamCenterFFTCPU.cpp | 16 ++++++--------- image_analysis/indexing/FFTIndexerCPU.cpp | 6 +++--- .../scale_merge/TranslationalNCS.cpp | 15 +++++++++++--- rugnux/ModelFFT.cpp | 10 ++++++---- viewer/charts/JFJochDatasetInfoChartView.cpp | 20 +++++++++++++------ 6 files changed, 55 insertions(+), 26 deletions(-) create mode 100644 common/FFTWPlannerLock.h diff --git a/common/FFTWPlannerLock.h b/common/FFTWPlannerLock.h new file mode 100644 index 000000000..a432b7971 --- /dev/null +++ b/common/FFTWPlannerLock.h @@ -0,0 +1,14 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +// FFTW's planner - every fftwf_plan_* and fftwf_destroy_plan call - shares global state and is not +// thread-safe; executing a plan is. So every piece of code that makes or destroys a plan takes this +// one lock, not a lock of its own: two private locks do not stop two planners running at once. +inline std::mutex &FFTWPlannerMutex() { + static std::mutex m; + return m; +} diff --git a/image_analysis/geom_refinement/BeamCenterFFTCPU.cpp b/image_analysis/geom_refinement/BeamCenterFFTCPU.cpp index 6c59f77d6..a9be3f21f 100644 --- a/image_analysis/geom_refinement/BeamCenterFFTCPU.cpp +++ b/image_analysis/geom_refinement/BeamCenterFFTCPU.cpp @@ -11,13 +11,9 @@ #include -namespace { +#include "../../common/FFTWPlannerLock.h" -// FFTW planning is not thread-safe (execution is); same idiom as FFTIndexerCPU. -std::mutex &FftwPlanMutex() { - static std::mutex m; - return m; -} +namespace { // One real-to-complex 2D FFT workspace of fixed padded size, reused for every forward transform. class Rfft2 { @@ -32,7 +28,7 @@ public: Rfft2(int64_t ny_, int64_t nx_) : ny(ny_), nx(nx_), nxc(nx_ / 2 + 1) { real.resize(static_cast(ny * nx)); spectrum.resize(static_cast(ny * nxc)); - std::unique_lock lock(FftwPlanMutex()); + std::unique_lock lock(FFTWPlannerMutex()); fwd = fftwf_plan_dft_r2c_2d(static_cast(ny), static_cast(nx), real.data(), reinterpret_cast(spectrum.data()), FFTW_ESTIMATE); @@ -45,7 +41,7 @@ public: Rfft2(const Rfft2 &) = delete; Rfft2 &operator=(const Rfft2 &) = delete; ~Rfft2() { - std::unique_lock lock(FftwPlanMutex()); + std::unique_lock lock(FFTWPlannerMutex()); if (fwd) fftwf_destroy_plan(fwd); if (bwd) @@ -136,7 +132,7 @@ BeamCenterConvSurfaces1D BeamCenterFFTCPU::LineSurfaces(const std::vector fftwf_plan fwd, bwd; { - std::unique_lock lock(FftwPlanMutex()); + std::unique_lock lock(FFTWPlannerMutex()); fwd = fftwf_plan_dft_r2c_1d(static_cast(nfft), in.data(), reinterpret_cast(A.data()), FFTW_ESTIMATE); bwd = fftwf_plan_dft_c2r_1d(static_cast(nfft), @@ -190,7 +186,7 @@ BeamCenterConvSurfaces1D BeamCenterFFTCPU::LineSurfaces(const std::vector out.Q = inverse(accQ); out.D = inverse(accD); { - std::unique_lock lock(FftwPlanMutex()); + std::unique_lock lock(FFTWPlannerMutex()); fftwf_destroy_plan(fwd); fftwf_destroy_plan(bwd); } diff --git a/image_analysis/indexing/FFTIndexerCPU.cpp b/image_analysis/indexing/FFTIndexerCPU.cpp index 8061b62f5..218c95cef 100644 --- a/image_analysis/indexing/FFTIndexerCPU.cpp +++ b/image_analysis/indexing/FFTIndexerCPU.cpp @@ -10,7 +10,7 @@ #include #include -static std::mutex fftw_plan_mutex; +#include "../../common/FFTWPlannerLock.h" static inline double dot_abs(const Coord& a, const Coord& b) { return std::fabs(a.x * b.x + a.y * b.y + a.z * b.z); @@ -35,7 +35,7 @@ FFTIndexerCPU::FFTIndexerCPU(const IndexingSettings& settings) int n[1] = { H }; { - std::unique_lock ul(fftw_plan_mutex); + std::unique_lock ul(FFTWPlannerMutex()); plan = fftwf_plan_many_dft_r2c( 1, n, nDirections, h_input_fft.data(), nullptr, 1, H, @@ -49,7 +49,7 @@ FFTIndexerCPU::FFTIndexerCPU(const IndexingSettings& settings) FFTIndexerCPU::~FFTIndexerCPU() { - std::unique_lock ul(fftw_plan_mutex); + std::unique_lock ul(FFTWPlannerMutex()); fftwf_destroy_plan(plan); } diff --git a/image_analysis/scale_merge/TranslationalNCS.cpp b/image_analysis/scale_merge/TranslationalNCS.cpp index 5fa3b1a96..7b5a91e07 100644 --- a/image_analysis/scale_merge/TranslationalNCS.cpp +++ b/image_analysis/scale_merge/TranslationalNCS.cpp @@ -15,6 +15,8 @@ #include +#include "../../common/FFTWPlannerLock.h" + namespace { // The Patterson is measured in a band, not over the whole dataset: below ~20 A the bulk solvent // and the beam stop dominate, above ~5 A the map answers a question about atoms rather than about @@ -231,8 +233,12 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector(coeff), map, - FFTW_ESTIMATE); + fftwf_plan plan; + { + std::lock_guard lock(FFTWPlannerMutex()); + plan = fftwf_plan_dft_c2r_3d(nx, ny, nz, reinterpret_cast(coeff), map, + FFTW_ESTIMATE); + } // The near-origin exclusion, as a set of grid points rather than a per-point distance test: only // the points inside a 15 A ball around each origin-equivalent point can be excluded, and there @@ -368,7 +374,10 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector, fftwf_plan> plans; - std::lock_guard lock(m); + std::lock_guard lock(FFTWPlannerMutex()); const std::array key{nu, nv, nw}; if (const auto it = plans.find(key); it != plans.end()) return it->second; diff --git a/viewer/charts/JFJochDatasetInfoChartView.cpp b/viewer/charts/JFJochDatasetInfoChartView.cpp index 068699d7d..abdfe4305 100644 --- a/viewer/charts/JFJochDatasetInfoChartView.cpp +++ b/viewer/charts/JFJochDatasetInfoChartView.cpp @@ -11,6 +11,7 @@ #include #include "JFJochDatasetInfoChartView.h" +#include "../../common/FFTWPlannerLock.h" namespace { // Qt Charts draws its curves 2 px wide, which is hard to follow across a desk. Take the weight @@ -965,18 +966,25 @@ void JFJochDatasetInfoChartView::buildFFTChart() { std::vector out(static_cast(nComplex)); - fftwf_plan plan = fftwf_plan_dft_r2c_1d( - n, - in.data(), - out.data(), - FFTW_ESTIMATE); + fftwf_plan plan; + { + std::lock_guard lock(FFTWPlannerMutex()); + plan = fftwf_plan_dft_r2c_1d( + n, + in.data(), + out.data(), + FFTW_ESTIMATE); + } if (!plan) { return; } fftwf_execute(plan); - fftwf_destroy_plan(plan); + { + std::lock_guard lock(FFTWPlannerMutex()); + fftwf_destroy_plan(plan); + } // Compute amplitude spectrum and frequencies (0 .. Nyquist) m_fftMagnitudes.resize(static_cast(nComplex));