From 20ab92fcaba21293209c67daefb5fc97832e235f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 00:05:36 +0200 Subject: [PATCH] Build: no FMA contraction; first-pass cell refinement independent of SIMD width The same source built with and without -march=x86-64-v3 gave different results on battery sets (9zmu axis-harmonic supercell arbiter fired in one build only; 8rud resolution cut 1.69 vs 1.70 A; myob_x06da_powder_1 cut 1.07 vs 1.42 A; 7mzt short-axis first pass; lcystine_x10sa_20keV indexed vs no lattice), against the rule that no decision may depend on compiler flags. Two mechanisms, found by building the merged tree four ways (baseline, x86-64-v2, x86-64-v3, x86-64-v3 -mno-fma) with and without -ffp-contract=off and comparing p.mtz: - FMA contraction. GCC contracts a*b+c whenever the target has FMA. With -ffp-contract=off the x86-64-v3 build gives a p.mtz byte-identical to the baseline build on 8 of 9 sets (7mzt, 8rud, 9zmu, insu_I_x06da_5keV_2, myob_x06da_powder_1, myob_x10sa, cytc_x10sa, thau_x10sa_16keV), on both the GPU and the CPU build. Cost: none measurable (user core-s, CPU build, x86-64-v3 vs the same with -ffp-contract=off: 1877/1874, 3285/3253, 2212/2195 on myob/cytc/thau; GPU likewise within noise). Set project-wide for C, C++ and CUDA host code; MSVC does not contract under /fp:precise. - SIMD width. lcystine still differed: baseline and x86-64-v2 (128-bit) agreed, x86-64-v3 with or without FMA (256-bit) disagreed - Eigen's HouseholderQR in the FFT indexer's candidate refinement (PostIndexingRefinement.cpp) reduces column norms over all spots in packets of the target width. Replaced by the 3x3 normal equations summed in spot order in double. All four builds now agree on all nine sets. Changes results of the default x86-64-v3 build (contraction off); lcystine_x10sa_20keV now gives no lattice in every build (its first pass is a knife-edge: 0/60 vs 9/60 validation frames before). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- CMakeLists.txt | 12 +++++++++++ .../indexing/PostIndexingRefinement.cpp | 20 +++++++++++++++---- 2 files changed, 28 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1caea9809..dfc4da29e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,18 @@ IF (MSVC) ADD_COMPILE_OPTIONS($<$:-Xcompiler=/Zc:preprocessor>) ENDIF() +# No contraction of a*b+c into a fused multiply-add. GCC contracts whenever the target has FMA, so a +# -march=x86-64-v3 build (CI, production) and a baseline build rounded differently, and knife-edge +# decisions - which candidate cell a first pass indexes with, a resolution cut, a symmetry gate - +# came out differently from the same source and data. Off, every x86-64 level computes the same bits +# (measured: identical p.mtz on nine battery sets, no measurable time cost), and aarch64, where FMA is +# always there and GCC contracts by default, follows the same arithmetic. MSVC does not contract under +# its default /fp:precise. Host code of .cu sources gets it through nvcc; device code is not affected. +IF (NOT MSVC) + ADD_COMPILE_OPTIONS($<$:-ffp-contract=off>) + ADD_COMPILE_OPTIONS($<$:-Xcompiler=-ffp-contract=off>) +ENDIF() + SET(JFJOCH_WRITER_ONLY OFF CACHE BOOL "Compile HDF5 writer only") SET(JFJOCH_INSTALL_DRIVER_SOURCE OFF CACHE BOOL "Install kernel driver source (ignored if building writer only; necessary for RPM building)") SET(JFJOCH_USE_CUDA ON CACHE BOOL "Compile Jungfraujoch with CUDA") diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index d8d210eb0..a31bb7de3 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -82,7 +82,6 @@ namespace { const unsigned nspots = spots.rows(); const unsigned ncells = scores.rows(); VectorX below{nspots}; - MatrixX3 sel{nspots, 3u}; Mx3 resid{nspots, 3u}; Mx3 miller{nspots, 3u}; M3 cell; @@ -111,9 +110,22 @@ namespace { break; threshold *= cifssr.threshold_contraction; - sel.colwise() = below; - HouseholderQR qr{sel.select(spots, .0f)}; - cell = qr.solve(sel.select(miller, .0f)); + // Least squares spots * cell = miller over the spots below the threshold, by the normal + // equations summed in spot order in double. A QR over all the spots reduces its column + // norms in SIMD packets of the target's width, so its cell - and which candidate a + // borderline first pass indexed with - depended on -march; these sums do not. + Matrix3d ata = Matrix3d::Zero(); + Matrix3d atb = Matrix3d::Zero(); + for (unsigned i = 0; i < nspots; i++) { + if (!below[i]) + continue; + for (int r = 0; r < 3; r++) + for (int c = 0; c < 3; c++) { + ata(r, c) += double(spots(i, r)) * double(spots(i, c)); + atb(r, c) += double(spots(i, r)) * double(miller(i, c)); + } + } + cell = (ata.inverse() * atb).cast(); } resid = CalculateResiduals(spots, cell);