From aca2b2fdc3d23050d0e949707eae27af4bb92d17 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 15 Oct 2026 15:14:51 +0200 Subject: [PATCH] CPU profile fit: the Gaussian's cells in a loop that vectorises, output unchanged Building each reflection's Gaussian profile was a third of the CPU pass B's samples (cytochrome C), one scalar BraggExp per cell. BraggExp branched to its two ends before the arithmetic, and the branches - with GCC's -ftrapping-math, which forbids running float arithmetic on a path that does not need it - kept every loop over it scalar. It now tells the ends apart on the bits of x, an integer comparison, runs the same arithmetic on x (or on 0 out of range) and masks the result: the same bits for every float but NaN, checked over all 2^32 inputs against the old function, scalar and vectorised. The cells are filled a row at a time and summed afterwards in the same order, so the loop vectorises (8 lanes, AVX2) and the profile is unchanged. BraggScale2 chooses its operands instead of a branch; scaling by 1 is exact. The Metal kernel shares the header and gets the same arithmetic. cytochrome C --no-gpu: p.hkl and MTZ data identical; user cycles -7% against the previous commit on the same (loaded, load 16-32) machine. Co-Authored-By: Claude Opus 5.5 --- image_analysis/bragg_integration/BraggExp.h | 31 ++++++++++++++----- .../BraggIntegrationEngineCPU.cpp | 13 ++++---- 2 files changed, 30 insertions(+), 14 deletions(-) diff --git a/image_analysis/bragg_integration/BraggExp.h b/image_analysis/bragg_integration/BraggExp.h index 818afd4dc..f4e4da92f 100644 --- a/image_analysis/bragg_integration/BraggExp.h +++ b/image_analysis/bragg_integration/BraggExp.h @@ -18,24 +18,38 @@ #ifdef __METAL_VERSION__ inline float BraggBitsToFloat(int i) { return as_type(i); } +inline int BraggFloatToBits(float f) { return as_type(f); } #else #include #include inline float BraggBitsToFloat(int32_t i) { return std::bit_cast(i); } +inline int32_t BraggFloatToBits(float f) { return std::bit_cast(f); } #endif -// p * 2^k, rounded once where the result is subnormal. +// p * 2^k, rounded once where the result is subnormal: below k = -125 the scaling is split in two. The +// two cases differ only in their operands - the normal one scales by 1 the second time, which is exact - +// so a loop over it has no branch and vectorises. inline float BraggScale2(float p, int k) { - if (k < -125) - return (p * BraggBitsToFloat((k + 100 + 127) << 23)) * BraggBitsToFloat((-100 + 127) << 23); - return p * BraggBitsToFloat((k + 127) << 23); + const bool subnormal = k < -125; + const float first = BraggBitsToFloat(((subnormal ? k + 100 : k) + 127) << 23); + const float second = subnormal ? BraggBitsToFloat((-100 + 127) << 23) : 1.0f; + return (p * first) * second; } // Following Cody & Waite (1980) Software Manual for the Elementary Functions; Dekker (1971) Numer. Math. // 18, 224-242; Knuth, The Art of Computer Programming Vol. 2, 4.2.2 (2Sum). -inline float BraggExp(float x) { - if (x < -104.0f) return 0.0f; // below half the smallest subnormal - if (x > 88.72283f) return BraggBitsToFloat(0x7f800000); // overflows +// The two ends are told apart on the bits of x and applied to the bits of the result, with integer masks +// rather than branches or selects: the compiler moves arithmetic behind a branch or a select that only +// one side needs, and may not run float arithmetic on a path that does not need it, which either way +// keeps a loop over BraggExp from vectorising. Out of range the arithmetic runs on 0 and its result is +// masked off; in range it is the same arithmetic on the same x. (A negative float's bits are a negative +// integer that grows as the float falls; a positive one's grow with it.) +inline float BraggExp(float x_in) { + const int bits = BraggFloatToBits(x_in); + const bool under = (bits < 0) & (bits > BraggFloatToBits(-104.0f)); // below half the smallest subnormal + const bool over = bits > BraggFloatToBits(88.72283f); // overflows + const int in_range = -static_cast(!(under | over)); // every bit set, or none + const float x = BraggBitsToFloat(bits & in_range); const float shifter = 12582912.0f; // 1.5 * 2^23: adding it rounds to an integer const float kf = (x * 1.44269504f + shifter) - shifter; const int k = static_cast(kf); @@ -59,5 +73,6 @@ inline float BraggExp(float x) { const float e1 = (rr - (s1 - bv1)) + (hz - bv1); const float t1 = 1.0f + s1; // 1 + s1 = t1 + e3 (Fast2Sum, |s1| < 1) const float e3 = (1.0f - t1) + s1; - return BraggScale2(t1 + (e3 + (e1 + c)), k); + const float e = BraggScale2(t1 + (e3 + (e1 + c)), k); + return BraggBitsToFloat((BraggFloatToBits(e) & in_range) | (static_cast(over) * 0x7f800000)); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 46ba4fefe..8f96fe036 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -486,13 +486,14 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne prefetch_row(img.address(static_cast(y) * W + gx0), img.address(static_cast(y) * W + gx1)); cells.resize(static_cast(Gf) * Gf); + // The cells a row at a time, a loop that vectorises; then their sum, in the same order. + for (int dy = -Rf; dy <= Rf; ++dy) { + float *row = cells.data() + static_cast(dy + Rf) * Gf; + for (int dx = -Rf; dx <= Rf; ++dx) + row[dx + Rf] = BraggGaussian(dx, dy, s.fx, s.fy, s.ux, s.uy, s.inv_2s2r, s.inv_2s2t); + } float gs = 0.0f; - for (int dy = -Rf; dy <= Rf; ++dy) - for (int dx = -Rf; dx <= Rf; ++dx) { - const float g = BraggGaussian(dx, dy, s.fx, s.fy, s.ux, s.uy, s.inv_2s2r, s.inv_2s2t); - cells[(dy + Rf) * Gf + (dx + Rf)] = g; - gs += g; - } + for (const float g : cells) gs += g; const float inv_gs = 1.0f / gs; for (float &p : cells) p *= inv_gs; P = cells.data();