From ab645847ede7b184b73e53f5ae3ea62f0524d1c8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 14 Oct 2026 10:32:03 +0200 Subject: [PATCH] FFT indexer GPU peak search ~45% faster per call, output unchanged calculate_fft_result (the per-direction prominence peak search) was 61-65% of the GPU kernel time of a de-novo stills run, and the card is what limits that loop. It is one dependent chain per direction, latency-bound. Two exact changes, taken where the spectrum holds a whole background window (any FFT max cell above ~30 A; the general loop stays for the rest): - the sliding window written out as its three phases (clamped at the start, centred, clamped at the end): the same additions and subtractions in the same order, so the running double sum is bit for bit the same, but the loads no longer sit behind data-dependent while loops and are issued ahead; - the window mean divides by a constant width, so the double division (a long FP64 iteration on a consumer card) is replaced by the reciprocal multiply plus one correction by the exact remainder, which returns the correctly rounded quotient (Markstein). Checked on the CPU against the division for 4e8 random dividends over all odd widths 3-401: 0 differences. Measured (RTX 5080, box loaded by other builds, load 12-25): kernel minimum 1.50 -> 0.83 ms, mean under the 8-stream pool 4.7 -> 2.8 ms; de-novo stills 2000 frames ABAB 113/120 -> 127/138 frames/s. p.hkl byte-identical (2000 and 3000 frames, base and new min-pix schemes). Co-Authored-By: Claude Opus 5.5 (1M context) --- docs/ACKNOWLEDGEMENT.md | 5 ++ image_analysis/indexing/FFTIndexerGPU.cu | 78 ++++++++++++++++-------- 2 files changed, 59 insertions(+), 24 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 10b01e011..33feb7705 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -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). diff --git a/image_analysis/indexing/FFTIndexerGPU.cu b/image_analysis/indexing/FFTIndexerGPU.cu index 36bbbe827..99ecbd4d2 100644 --- a/image_analysis/indexing/FFTIndexerGPU.cu +++ b/image_analysis/indexing/FFTIndexerGPU.cu @@ -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(j); - if (len <= min_length_A) continue; + if (len <= min_length_A) return; const float mag = mag_at(j); - const float bg = static_cast(winsum / static_cast(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(__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(winsum / static_cast(whi - wlo + 1))); + } } d_results[i] = result; // Store the result }