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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB

CMakeLists.txt

@@ -41,6 +41,18 @@ IF (MSVC)
     ADD_COMPILE_OPTIONS($<$<COMPILE_LANGUAGE:CUDA>:-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($<$<COMPILE_LANGUAGE:C,CXX>:-ffp-contract=off>)
+    ADD_COMPILE_OPTIONS($<$<COMPILE_LANGUAGE:CUDA>:-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")

image_analysis/indexing/PostIndexingRefinement.cpp

@@ -82,7 +82,6 @@ namespace {
         const unsigned nspots = spots.rows();
         const unsigned ncells = scores.rows();
         VectorX<bool> below{nspots};
-        MatrixX3<bool> 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<Mx3> 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<float>();
             }
 
             resid = CalculateResiduals(spots, cell);