From 482f3cbaa21b2b91987455b147f9a5478427d889 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Sat, 5 Sep 2026 08:05:03 -0400
Subject: [PATCH] Real part comparisons

---
 Grid/lattice/Lattice_comparison.h | 63 +++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)

diff --git a/Grid/lattice/Lattice_comparison.h b/Grid/lattice/Lattice_comparison.h
index 91d683da8..d1067a968 100644
--- a/Grid/lattice/Lattice_comparison.h
+++ b/Grid/lattice/Lattice_comparison.h
@@ -175,5 +175,68 @@ template<class lobj,class robj>
 inline Lattice<vPredicate<robj> > operator != (const lobj & lhs, const Lattice<robj> & rhs) {
   return SLComparison(vne<lobj,robj>(),lhs,rhs);
 }
+
+//////////////////////////////////////////////////////////////////////////
+// Real-part relational comparison for COMPLEX lattices.
+//
+// Complex has no ordering, so operator<,>,<=,>= are (deliberately) undefined for
+// complex operands -- Comparison() in Lattice_comparison_utils.h is IfNotComplex --
+// and where() cannot be driven by a complex lattice directly.  It is however often
+// useful to threshold on the REAL PART of a complex (scalar/singlet) field, e.g. a
+// momentum-magnitude mask phat^2 > pc^2.  These free functions compare the real part
+// to a real threshold and return the matching IntegerPredicate, so the result feeds
+// where() exactly like the built-in relationals.
+//
+// Written with the per-lane getlane/putlane accessors (as in FFT.h, PaddedCell.h),
+// not extract/merge buffers.  One wrinkle: the predicate type IntegerPredicate<CComplex>
+// is vInteger, whose Nsimd (>= the widest real type's) exceeds the complex operand's
+// Nsimd by s = Npred/Nsimd.  where() only reads the ii=0 representative of each group,
+// at physical lane lane*s (the "s-fold skip" -- cf. extract()'s getlane(i*s)), so
+// filling the rest of the group is not functionally required; we replicate the value
+// across all s lanes anyway because some Grid code asserts the s replicas are equal
+// (and it mirrors what merge() does).  When vInteger is reworked to carry the operand
+// Nsimd (scope later) s becomes 1 and the inner loop drops out.
+//////////////////////////////////////////////////////////////////////////
+template<class scalar> class sRealLt { public:
+  accelerator_inline Integer operator()(const scalar &a, RealD b) const { return a.real() <  b ? 1 : 0; } };
+template<class scalar> class sRealLe { public:
+  accelerator_inline Integer operator()(const scalar &a, RealD b) const { return a.real() <= b ? 1 : 0; } };
+template<class scalar> class sRealGt { public:
+  accelerator_inline Integer operator()(const scalar &a, RealD b) const { return a.real() >  b ? 1 : 0; } };
+template<class scalar> class sRealGe { public:
+  accelerator_inline Integer operator()(const scalar &a, RealD b) const { return a.real() >= b ? 1 : 0; } };
+
+template<class sfunctor,class CComplex>
+inline Lattice<vPredicate<CComplex> > RealPartComparison(sfunctor op,const Lattice<CComplex> &lhs, RealD thr)
+{
+  Lattice<vPredicate<CComplex> > ret(lhs.Grid());
+  autoView( lv, lhs, AcceleratorRead);
+  autoView( rv, ret, AcceleratorWrite);
+  typedef typename CComplex::vector_type vsimd;
+  const int Nsimd = vsimd::Nsimd();
+  const int s     = IntegerPredicate<CComplex>::Nsimd() / Nsimd;   // lane-replication factor (see note)
+  accelerator_for(ss, lhs.Grid()->oSites(), Nsimd, {
+    vsimd v = TensorRemove(lv[ss]);                 // strip iScalar nest to the bare complex SIMD word
+#ifdef GRID_SIMT
+    { int lane = acceleratorSIMTlane(Nsimd);        // GPU: this thread == this operand lane
+#else
+    for(int lane=0;lane<Nsimd;lane++){              // CPU: walk the packed operand lanes
+#endif
+      Integer p = op(v.getlane(lane), thr);
+      for(int ii=0;ii<s;ii++) rv[ss]._internal.putlane(p, lane*s+ii);   // replicate across the group
+#ifdef GRID_SIMT
+    }
+#else
+    }
+#endif
+  });
+  return ret;
+}
+
+template<class CComplex> inline Lattice<vPredicate<CComplex> > RealPartLessThan    (const Lattice<CComplex> &a, RealD b){ return RealPartComparison(sRealLt<typename CComplex::vector_type::scalar_type>(),a,b); }
+template<class CComplex> inline Lattice<vPredicate<CComplex> > RealPartLessEqual    (const Lattice<CComplex> &a, RealD b){ return RealPartComparison(sRealLe<typename CComplex::vector_type::scalar_type>(),a,b); }
+template<class CComplex> inline Lattice<vPredicate<CComplex> > RealPartGreaterThan  (const Lattice<CComplex> &a, RealD b){ return RealPartComparison(sRealGt<typename CComplex::vector_type::scalar_type>(),a,b); }
+template<class CComplex> inline Lattice<vPredicate<CComplex> > RealPartGreaterEqual (const Lattice<CComplex> &a, RealD b){ return RealPartComparison(sRealGe<typename CComplex::vector_type::scalar_type>(),a,b); }
+
 NAMESPACE_END(Grid);
 #endif