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 }