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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
leonarski_fandClaude Opus 5 committed 2026-09-20 18:45:03 +02:00
1 parent 13de0e7a87
commit 7dea9eaf9c
6 files changed
+55 -26

No files matched your search

+14
View File
@@ -0,0 +1,14 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <mutex>
// 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;
}
@@ -11,13 +11,9 @@
#include <fftw3.h>
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<size_t>(ny * nx));
spectrum.resize(static_cast<size_t>(ny * nxc));
std::unique_lock lock(FftwPlanMutex());
std::unique_lock lock(FFTWPlannerMutex());
fwd = fftwf_plan_dft_r2c_2d(static_cast<int>(ny), static_cast<int>(nx), real.data(),
reinterpret_cast<fftwf_complex *>(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<float>
fftwf_plan fwd, bwd;
{
std::unique_lock lock(FftwPlanMutex());
std::unique_lock lock(FFTWPlannerMutex());
fwd = fftwf_plan_dft_r2c_1d(static_cast<int>(nfft), in.data(),
reinterpret_cast<fftwf_complex *>(A.data()), FFTW_ESTIMATE);
bwd = fftwf_plan_dft_c2r_1d(static_cast<int>(nfft),
@@ -190,7 +186,7 @@ BeamCenterConvSurfaces1D BeamCenterFFTCPU::LineSurfaces(const std::vector<float>
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);
}
+3 -3
View File
@@ -10,7 +10,7 @@
#include <mutex>
#include <vector>
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);
}
@@ -15,6 +15,8 @@
#include <fftw3.h>
#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<MergedReflectio
result.refusal = "could not allocate the Patterson grid";
return result;
}
fftwf_plan plan = fftwf_plan_dft_c2r_3d(nx, ny, nz, reinterpret_cast<fftwf_complex *>(coeff), map,
FFTW_ESTIMATE);
fftwf_plan plan;
{
std::lock_guard lock(FFTWPlannerMutex());
plan = fftwf_plan_dft_c2r_3d(nx, ny, nz, reinterpret_cast<fftwf_complex *>(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<MergedReflectio
}
}
fftwf_destroy_plan(plan);
{
std::lock_guard lock(FFTWPlannerMutex());
fftwf_destroy_plan(plan);
}
fftwf_free(coeff);
fftwf_free(map);
result.measurable = true;
+6 -4
View File
@@ -13,15 +13,17 @@
#include "gemmi/fail.hpp"
#include "../common/FFTWPlannerLock.h"
namespace {
// The r2c plan for an (nu, nv, nw) map stored u fastest, halving w - the layout gemmi's own transform
// reads and writes. FFTW's planner is not thread-safe, so plans are made under a lock; executing one
// on new arrays (fftwf_execute_dft_r2c) is, which is how every caller uses it.
// reads and writes. FFTW's planner is not thread-safe, so plans are made under the process-wide planner
// lock (which also guards this cache); executing one on new arrays (fftwf_execute_dft_r2c) is, which is
// how every caller uses it.
fftwf_plan PlanFor(int nu, int nv, int nw) {
static std::mutex m;
static std::map<std::array<int, 3>, fftwf_plan> plans;
std::lock_guard lock(m);
std::lock_guard lock(FFTWPlannerMutex());
const std::array<int, 3> key{nu, nv, nw};
if (const auto it = plans.find(key); it != plans.end())
return it->second;
+14 -6
View File
@@ -11,6 +11,7 @@
#include <cmath>
#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<fftwf_complex> out(static_cast<size_t>(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<size_t>(nComplex));