Merge branch 'dq-flag-stability' into overnight-cand

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);