mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

Merge FFT indexer GPU peak search rewrite into rc175: ~45% faster per call, output unchanged

Browse Source 
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-10-14 13:15:17 +02:00
co-authored by Claude Opus 5.5
parent 02ea56b2d3 ab645847ed
commit 51397c15f4
2 changed files with 59 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/ACKNOWLEDGEMENT.md
+5
View File
@@ -518,6 +518,11 @@ components with a parallel union-find, following D. P. Playne and K. Hawick, "A
Parallel Connected-Component Labelling on GPUs" (2018), IEEE Trans. Parallel Distrib. Syst. 29,
1217-1230 [doi:10.1109/TPDS.2018.2799216](https://doi.org/10.1109/TPDS.2018.2799216).
The FFT indexer's GPU peak search divides by a fixed window width with a reciprocal multiplication
corrected once by the exact remainder, which returns the correctly rounded quotient, following
P. W. Markstein, "Computation of elementary functions on the IBM RISC System/6000 processor"
(1990), IBM J. Res. Dev. 34, 111-119 (DOI not verified here).
This software uses the Viridis, Magma and Inferno colormaps from Matplotlib under its
BSD-compatible license. J. D. Hunter, "Matplotlib: A 2D graphics environment" (2007), Comput. Sci.
Eng. 9, 90-95 [doi:10.1109/MCSE.2007.55](https://doi.org/10.1109/MCSE.2007.55).
image_analysis/indexing/FFTIndexerGPU.cu
+54 -24
View File
@@ -61,41 +61,71 @@ __global__ void calculate_fft_result(
// exceed the true lattice peaks, so a plain argmax|spec| returns a short envelope
// vector on weak/pink-beam frames. Subtracting the local mean removes that envelope
// while sharp lattice peaks keep their height (mirrors FFTIndexerCPU).
double winsum = 0.0;
int wlo = 0;
int whi = min(bg_half, out_len - 1);
for (int k = wlo; k <= whi; ++k)
winsum += mag_at(k);
float best_prom = 0.0f;
FFTResult result{.magnitude = 0.0f, .direction = i, .length = -1};
for (int j = 0; j < out_len; ++j) {
// Constant-width window, slid inward at the ends instead of truncated (see
// FFTIndexerCPU): a peak within bg_half of either end - where the LONGEST cells sit -
// otherwise gets a one-sided background and a biased prominence. Both bounds stay
// monotonically non-decreasing in j, so the running sum below is still valid.
int want_lo = j - bg_half;
int want_hi = j + bg_half;
if (want_lo < 0) { want_hi = min(out_len - 1, want_hi - want_lo); want_lo = 0; }
if (want_hi > out_len - 1) {
want_lo = max(0, want_lo - (want_hi - (out_len - 1)));
want_hi = out_len - 1;
}
while (whi < want_hi) { ++whi; winsum += mag_at(whi); }
while (wlo < want_lo) { winsum -= mag_at(wlo); ++wlo; }
// bg is consider()'s background, the window mean - formed by whichever loop below applies.
const auto consider = [&](int j, float bg) {
const float len = len_coeff * static_cast<float>(j);
if (len <= min_length_A) continue;
if (len <= min_length_A) return;
const float mag = mag_at(j);
const float bg = static_cast<float>(winsum / static_cast<double>(whi - wlo + 1));
const float prom = mag - bg;
if (prom > best_prom) {
best_prom = prom;
result.magnitude = prom;
result.length = len;
}
};
if (out_len - 1 >= 2 * bg_half) {
// Where the spectrum holds a whole window, the window of the general loop below is
// [0, 2 bg_half] until it can be centred, centred while it fits, and the last 2 bg_half + 1
// bins after that: the running sum takes the same additions and subtractions in the same
// order, so it is bit for bit the same number. Spelled out, the loads no longer wait on
// data-dependent loops and the compiler can issue them ahead, which is what the search -
// latency-bound, one dependent chain per direction - is short of.
// The window width is fixed here, so the mean divides by one number throughout. A double
// division is a long iteration on the FP64 units a consumer card has few of, and it was a
// quarter of this kernel's time; multiplying by the correctly rounded reciprocal and correcting the
// quotient once with its exact remainder gives the correctly rounded quotient - the same
// double the division returns (Markstein (1990) IBM J. Res. Dev. 34, 111-119).
const double count = 2 * bg_half + 1;
const double inv_count = 1.0 / count;
double winsum = 0.0;
for (int k = 0; k <= 2 * bg_half; ++k)
winsum += mag_at(k);
#pragma unroll 4
for (int j = 0; j < out_len; ++j) {
if (j > bg_half && j + bg_half <= out_len - 1) {
winsum += mag_at(j + bg_half);
winsum -= mag_at(j - bg_half - 1);
}
const double q = __dmul_rn(winsum, inv_count);
consider(j, static_cast<float>(__fma_rn(__fma_rn(-q, count, winsum), inv_count, q)));
}
} else {
double winsum = 0.0;
int wlo = 0;
int whi = min(bg_half, out_len - 1);
for (int k = wlo; k <= whi; ++k)
winsum += mag_at(k);
for (int j = 0; j < out_len; ++j) {
// Constant-width window, slid inward at the ends instead of truncated (see
// FFTIndexerCPU): a peak within bg_half of either end - where the LONGEST cells sit -
// otherwise gets a one-sided background and a biased prominence. Both bounds stay
// monotonically non-decreasing in j, so the running sum below is still valid.
int want_lo = j - bg_half;
int want_hi = j + bg_half;
if (want_lo < 0) { want_hi = min(out_len - 1, want_hi - want_lo); want_lo = 0; }
if (want_hi > out_len - 1) {
want_lo = max(0, want_lo - (want_hi - (out_len - 1)));
want_hi = out_len - 1;
}
while (whi < want_hi) { ++whi; winsum += mag_at(whi); }
while (wlo < want_lo) { winsum -= mag_at(wlo); ++wlo; }
consider(j, static_cast<float>(winsum / static_cast<double>(whi - wlo + 1)));
}
}
d_results[i] = result; // Store the result
}