no compile fix

Grid/lattice/Lattice_transfer.h

@@ -191,6 +191,15 @@ accelerator_inline void convertType(vComplexD & out, const ComplexD & in) {
 accelerator_inline void convertType(vComplexD2 & out, const ComplexD & in) {
   ((ComplexD*)&out)[acceleratorSIMTlane(vComplexD::Nsimd()*2)] = in;
 }
+// Lex chart: one lane, so there is no lane to select. Needed because the
+// coalesced accessor hands the kernel a scalar object while the destination
+// is still the vector type.
+accelerator_inline void convertType(sComplexF & out, const ComplexF & in) {
+  out.v = in;
+}
+accelerator_inline void convertType(sComplexD & out, const ComplexD & in) {
+  out.v = in;
+}
 #endif
 
 accelerator_inline void convertType(vComplexF & out, const vComplexD2 & in) {

Grid/simd/Simd.h

@@ -42,8 +42,12 @@ template <typename T>
 struct RealPart {
   typedef T type;
 };
+// Grid's own complex alias, so this matches thrust::complex on CUDA/HIP as
+// well as std::complex on the host. Qualifying it std:: left Real == the
+// complex type itself on device builds, which removes the Grid_simd(Real)
+// constructor and so any assignment of a real to a complex lattice.
 template <typename T>
-struct RealPart<std::complex<T> > {
+struct RealPart<complex<T> > {
   typedef T type;
 };
 

TODO

@@ -1,3 +1,32 @@
+SIMD lane extent audit (2026-08-19)
+
+i) Remove the explicit lattice types from GaugeGroup.h; take a Gimpl or other
+   policy template instead.
+
+ -- LieAlgebraProject(LatticeAlgebraMatrix &, const LatticeMatrix &, int) names the
+    default precision typedefs and launches accelerator_for with vComplex::Nsimd()
+    as the lane extent (GaugeGroup.h:438).
+ -- Correct today ONLY because the signature pins the precision. LatticeMatrixF/D and
+    LatticeAlgebraMatrixF/D exist alongside (GaugeGroup.h:131-145); an overload or a
+    retype would silently mismatch the launch geometry under GRID_SIMT.
+ -- Policy template, and derive the lane extent from the field, not from a named type.
+
+ -- Audit context: 258 accelerator_for* call sites; this is the ONLY library one whose
+    lane extent names a type. 8 more in benchmarks/Benchmark_su3_gpu.cc. 45 derive from
+    vobj::Nsimd() and are clean.
+
+ii) acceleratorSIMTlane(int Nsimd) ignores its argument -- every definition returns the
+    thread index. The lane a kernel reads is set solely by the accelerator_for launch
+    extent, so the argument is decorative and enforces nothing. Remove it.
+    The (threadIdx.x % Nsimd) idea is too constraining.
+
+iii) GeneralCoarsenedMatrix.h:455,642 launch with nbasis lanes and use the SIMT dimension
+     as a BASIS index rather than a SIMD lane. Works only because the launch extent is
+     nbasis. No coalescedRead/coalescedWrite may appear inside such a kernel.
+     Retired with V1: V2 launches nsimd=1 with an explicit j loop.
+
+-------
+
 i)    Refine subspace with HDCG & recompute
 ii)   Block Lanczos in coarse space
 iii)  Batched block project in the operator computation

skills/multigrid-design-notes.md

@@ -0,0 +1,167 @@
+---
+name: multigrid-design-notes
+description: "Design, development and tuning of LQCD multigrid solvers for physical-mass Möbius DWF on Frontier (AMD MI250X); covers HDCG, PVdagM two-level solver, coarse operator performance, and Lüscher deflation of the coarse solve"
+metadata: 
+  node_type: memory
+  type: project
+  originSessionId: cc1844e3-ab6f-4425-bf7e-a091b9554290
+---
+
+# LQCD Multigrid: Design, Development and Tuning
+
+## Physical problem
+
+Physical-mass Möbius DWF: Ls=24, b=1.5, c=0.5, M5=1.8, mass=0.00078.
+48³×96 lattice. MPI geometry 3.6.4.4 (288 ranks / GCDs on Frontier).
+Target: accelerate HMC fermion force and CG solves.
+
+## Two solver paths
+
+### Path 1: HDCG (TwoLevelADEF2 on MdagM)
+
+- File: `examples/Example_mdagm.cc`
+- Operator: MdagM (Hermitian positive definite); outer solver is ADEF2 CG.
+- Subspace: 60 near-null vectors via CG inverse iteration (`CreateSubspace`).
+- Coarse geometry: block {4,4,3,4}, coarse lattice 12×12×16×24, Ls_coarse=1.
+- Coarse operator: `GeneralCoarsenedMatrix` (Petrov-Galerkin, npoint=33, NextToNearestStencil).
+- Smoother: fixed-iteration CG on shifted operator (M†M + lo), lo=hi/80, 20 iters.
+- Coarse solve: CG with DeflatedGuesser using 60 chi deflation vectors.
+- Chi vectors: extracted BEFORE block-GS by diagonalising W_ij=<ψ_i|M†M|ψ_j>, computing
+  chi_k = Σ_i V[i,k] ψ_i. These are global near-null combinations; block-GS destroys this.
+- Deflation effect: 1089 coarse CG iters (undeflated) → 329 with 60 chi vectors.
+- Best result: 274 outer ADEF2 iters, ~400s on 288 ranks Frontier.
+- Coarse MVM performance: 541.7 GFlop/s kernel, 1119.6 GB/s (70% HBM), 97% roofline.
+  MPI latency = 1188 μs = 37% of 3.2 ms per coarse MVM call. Single-RHS is bandwidth-bound.
+
+### Path 2: PVdagM two-level PGCR
+
+- File: `examples/Example_pvdagm.cc` and `examples/Example_pvdagm_defl.cc`
+- Operator: PVdagM = PV†M (non-Hermitian); outer solver is PGCR.
+- PV is the Pauli-Villars (mass=1) Möbius operator. PVdagM has exact zero modes.
+- Subspace: 60 near-null vectors via GCR inverse iteration (`CreateSubspaceGCR`).
+  GCR setup is slow: each vector takes O(600) PGCR steps, total ~4100s setup on Frontier.
+- Coarse geometry: block 2^4 (lattice halved each dim), Ls_coarse=1.
+- Coarse operator: `GeneralCoarsenedMatrix` with non-Hermitian coarsening.
+- Preconditioner: `MGPreconditioner` V-cycle (pre-smooth, project, coarse solve, promote, post-smooth).
+- Smoother: PGCR on ShiftedPVdagM (shift=0.01).
+- Baseline (no deflation, 5e-2 coarse tol): 59 outer iters, 300s solve time.
+  Coarse solve: 4.58s/call, 250 PGCR steps, NEVER converges ("did not converge" every call).
+- Outer iteration count vs coarse tolerance: 5e-2→59 iters, 1e-1→63 iters, 3e-2→34 iters.
+  But at 3e-2 without deflation: 1000+ coarse PGCR steps per call (useless).
+
+## Key implementation work: GeneralCoarsenedMatrix performance
+
+`Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h`:
+
+1. **accelerator_barrier fix**: `acceleratorBarrier()` is not a Grid macro; correct call is
+   `accelerator_barrier(dummy)` (takes a dummy argument). This caused SIGBUS on Frontier.
+
+2. **Coalesced FT kernel**: In `CoarsenOperator`, the loop filling A_v[sss](i,j) was serialised
+   over j. Changed to:
+   ```cpp
+   accelerator_for(sss, osites, nbasis, {
+       int j = acceleratorSIMTlane(nbasis);
+       A_v[sss](i,j) = FT_v[sss](j);
+   });
+   ```
+   This gives coalesced HBM access (nbasis consecutive elements per warp lane).
+
+3. **Batched CoarsenOperator**: Used `MultiRHSBlockProject` to batch all npoint=33 stencil
+   directions in one GEMM call per basis vector, replacing serial blockProject calls.
+   Reduced projection from dominant bottleneck to 12% of CoarsenOperator time.
+   mat (linop applications) now dominates at 83%.
+
+## Lüscher deflation of the coarse solve (Example_pvdagm_defl.cc)
+
+### Theory (Lüscher arXiv:0706.2298, Section A.3)
+
+For near-null vectors {ψ_s} of operator D, the Petrov-Galerkin initial guess is:
+  guess = Ψ W⁻¹ Ψ† src
+where W_st = <ψ_s|D|ψ_t> and Ψ is the matrix of ψ columns.
+
+Condition <ψ_s | src - D*guess> = 0 gives W c = b, b_t = <ψ_t|src>.
+No SVD needed — W is dense, invert directly (LU). The U,V from SVD are unitaries
+within the ψ-basis and cancel in W⁻¹; direct inverse is cleaner.
+
+### Diagnostic results (job 4948520, before deflation)
+
+Fine projected matrix W (60×60):
+- ||W|| = 0.02399 — near-null vectors are genuinely small.
+- Singular values: range [0.00169, 0.00529], ratio ~3:1. Well-conditioned inverse.
+
+Coarse null matrix C_kl = <P ψ_k | A_coarse | P ψ_l>:
+- ||C|| = 0.02399 — **identical to ||W||**. Galerkin property is exact.
+- ||C - C†|| / ||C|| = 2.59e-9 — **C is Hermitian to machine precision**.
+  Despite PVdagM being non-Hermitian, the projected coarse matrix is numerically Hermitian.
+- C singular values match W singular values exactly (coarsening faithful).
+
+This means: the coarse near-null vectors (projections of fine ψ_k) are exact near-null
+vectors of the coarse operator, AND the deflation can use a Hermitian eigensolver.
+
+### Implementation
+
+**CoarseDeflatedGuesser** (in Example_pvdagm_defl.cc, before MGPreconditioner):
+```cpp
+template<class Field>
+class CoarseDeflatedGuesser : public LinearFunction<Field> {
+    const std::vector<Field> &chi;   // coarse eigenvectors of C_sym
+    const std::vector<RealD> &eval;  // eigenvalues
+public:
+    void operator()(const Field &src, Field &guess) {
+        guess = Zero();
+        for (int k = 0; k < chi.size(); k++)
+            axpy(guess, TensorRemove(innerProduct(chi[k], src)) / eval[k], chi[k], guess);
+    }
+};
+```
+
+**Chi vector construction** (in runMG, after CoarsenOperator and C computation):
+1. Project pre-GS fine subspace to coarse grid: psi_coarse[k] = P ψ_k
+2. Compute C_sym = (C + C†)/2
+3. SelfAdjointEigenSolver(C_sym) → eigenvalues lambda[k], eigenvectors V
+4. chi_coarse[k] = Σ_i V[i,k] * psi_coarse[i]  (satisfies <chi_j|A_c|chi_k> = lambda_k δ_jk)
+5. chi_eval[k] = lambda[k]
+
+**MGPreconditioner** modified to hold `CoarseSolver &_CoarseGuesser`:
+- Added as constructor parameter and member reference.
+- In V-cycle: replaced `Csol = Zero()` with `_CoarseGuesser(Csrc, Csol)` before `_CoarseSolve`.
+
+**Note on pre-GS subspace**: The subspace is copied at the top of runMG BEFORE CoarsenOperator
+modifies it via block-GS. This pre-GS copy is essential for computing W (fine projected matrix)
+and psi_coarse. After CoarsenOperator, AggregatesPD.subspace is the block-orthonormal basis.
+
+### Expected benefit
+Without deflation at 3e-2 coarse tolerance: 1000+ PGCR steps per call, never converges.
+With deflation: near-null modes removed by guesser; PGCR sees well-conditioned complement.
+Target: O(100-200) coarse PGCR steps per call, 34 outer iterations, large total time reduction.
+
+## Memory considerations (Frontier, MI250X, 64 GB HBM per GCD)
+
+Each fine fermion field: ~162 MB per rank (Ls=24, 48³×96/288 local vol, spincolour=12 complex).
+60 fine fields (subspace): ~9.7 GB.
+The pre-GS subspace COPY in runMG adds another 9.7 GB.
+Chi/null vector construction: build one at a time into a temp field, project immediately.
+Do NOT allocate 2×nbasis additional fine fields (would add 19 GB → OOM).
+
+Coarse vectors (~2 MB each): 120 coarse vectors (psi_coarse + chi_coarse) = ~240 MB. Fine.
+
+## Roofline analysis (coarse MVM, single-RHS)
+
+Arithmetic intensity: ~0.5 flops/byte. Roofline crossover: 119.7 flops/byte.
+Deeply memory-bandwidth bound. Peak HBM: 1600 GB/s. Achieved: 1119.6 GB/s (70%).
+MPI latency: 1188 μs = 37% of 3.2 ms per call. Irreducible for single-RHS.
+Multi-RHS is the fundamental solution for throughput, but cannot be used in HMC
+(each trajectory has a new gauge field → new coarse operator).
+
+## Three-level perspective
+
+The chi deflation of the coarse solve can be viewed as a third level: coarsening all
+the way to 1⁴ × nbasis. The 60 near-null vectors of the coarse operator span this
+third level's null space. The Lüscher guesser is the exact inverse on this space.
+
+## References
+
+- HDCG paper (2014): arXiv:1409.xxxx (P. Boyle) — ADEF2 CG with coarse deflation.
+- Lüscher (2007): arXiv:0706.2298 — non-Hermitian deflation, Section A.3.
+- Physical DWF multigrid: arXiv:2409.03904.
+- Grid library: github.com/paboyle/Grid.

tests/lexLattice/Makefile.am

@@ -0,0 +1 @@
+include Make.inc

tests/lexLattice/Test_GaugeAction_lex.cc

@@ -0,0 +1,128 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_GaugeAction_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Wilson loops on a lexicographic gauge field.
+//
+// The configuration is generated on the vectorised grid, written as NERSC,
+// and read back into the lexicographic layout: the reader recomputes the
+// plaquette and link trace and checks them against the header written by
+// the vectorised chart, so the crossing is validated inside Grid's own IO.
+// Plaquette, link trace, rectangle and the staple identity are then compared
+// between the two charts.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  typedef WilsonLoops<PeriodicGimplD>    vWL;
+  typedef WilsonLoops<lexPeriodicGimplD> lexWL;
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate lsimd({1,1,1,1});
+  Coordinate mpi   = GridDefaultMpi();
+
+  GridCartesian vGrid(latt,vsimd,mpi);
+  GridCartesian lGrid(latt,lsimd,mpi);
+
+  std::cout << GridLogMessage << "vectorised    Nsimd = " << vGrid.Nsimd() << std::endl;
+  std::cout << GridLogMessage << "lexicographic Nsimd = " << lGrid.Nsimd() << std::endl;
+  GRID_ASSERT( lGrid.Nsimd() == 1 );
+
+  GridParallelRNG pRNG(&vGrid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  LatticeGaugeFieldD Umu(&vGrid);
+  SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  std::string file("./ckpoint_lex.4000");
+  NerscIO::writeConfiguration(Umu,file,0,0);
+
+  lexLatticeGaugeFieldD Ulex(&lGrid);
+  FieldMetaData header;
+  NerscIO::readConfiguration(Ulex,header,file);
+
+  //////////////////////////////////////////////////
+  // Plaquette, link trace, rectangle in both charts
+  //////////////////////////////////////////////////
+  RealD vplaq = vWL::avgPlaquette(Umu);
+  RealD lplaq = lexWL::avgPlaquette(Ulex);
+
+  RealD vlink = vWL::linkTrace(Umu);
+  RealD llink = lexWL::linkTrace(Ulex);
+
+  RealD vrect = vWL::avgRectangle(Umu);
+  RealD lrect = lexWL::avgRectangle(Ulex);
+
+  std::cout << GridLogMessage << "plaquette   simd " << vplaq << " lex " << lplaq
+            << " header " << header.plaquette << std::endl;
+  std::cout << GridLogMessage << "link trace  simd " << vlink << " lex " << llink
+            << " header " << header.link_trace << std::endl;
+  std::cout << GridLogMessage << "rectangle   simd " << vrect << " lex " << lrect << std::endl;
+
+  GRID_ASSERT( fabs(vplaq-lplaq) < tol );
+  GRID_ASSERT( fabs(vlink-llink) < tol );
+  GRID_ASSERT( fabs(vrect-lrect) < tol );
+
+  //////////////////////////////////////////////////
+  // Plaquette via staples, in the lexicographic chart
+  //////////////////////////////////////////////////
+  {
+    RealD vol = lGrid.gSites();
+    RealD stap_plaq = 0.0;
+
+    lexLatticeColourMatrixD stap(&lGrid);
+    lexLatticeColourMatrixD Ul(&lGrid);
+    lexLatticeComplexD      stap_tr(&lGrid);
+
+    for(int mu=0;mu<Nd;mu++){
+      Ul = PeekIndex<LorentzIndex>(Ulex,mu);
+      lexWL::Staple(stap,Ulex,mu);
+      stap_tr = trace(Ul*stap);
+      auto Ts = sum(stap_tr);
+      stap_plaq += real(TensorRemove(Ts));
+    }
+    RealD StapScale = 1.0/vol/6.0/Nc/4.0;
+    RealD plaq_from_staples = stap_plaq*StapScale;
+
+    std::cout << GridLogMessage << "plaquette via staples (lex) " << plaq_from_staples
+              << " direct " << lplaq << std::endl;
+    GRID_ASSERT( fabs(plaq_from_staples - lplaq) < 1.0e-8 );
+  }
+
+  std::cout << GridLogMessage << "Test_GaugeAction_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_RectPlaq_lex.cc

@@ -0,0 +1,216 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_RectPlaq_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Plaquette and 2x1 rectangle built from covariant shifts, checked against
+// WilsonLoops, plus link trace and a blockSum coarsening.
+//
+// The measurement is written once against a gauge implementation and run in
+// both the vectorised and the lexicographic chart; every number must agree.
+// The configuration crosses layouts through NERSC IO, whose reader verifies
+// the header written by the other chart.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+template<class Gimpl>
+struct Measurements
+{
+  RealD plaq_shift;
+  RealD plaq_loops;
+  RealD plaq_staple;
+  RealD rect_shift;
+  RealD rect_loops;
+  RealD link;
+  RealD coarse_plaq;
+};
+
+template<class Gimpl>
+void Measure(typename Gimpl::Field &Umu,GridBase *coarse,Measurements<Gimpl> &m)
+{
+  typedef typename Gimpl::LinkField    LinkField;
+  typedef typename Gimpl::ComplexField ComplexField;
+  typedef WilsonLoops<Gimpl>           WL;
+
+  GridBase *grid = Umu.Grid();
+  RealD vol = grid->gSites();
+
+  std::vector<LinkField> U(Nd,grid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+
+  ///////////////////////////////////////////////////
+  // Link trace
+  ///////////////////////////////////////////////////
+  ComplexField LinkTrace(grid);
+  LinkTrace = Zero();
+  for(int mu=0;mu<Nd;mu++){
+    LinkTrace = LinkTrace + trace(U[mu]);
+  }
+  m.link = real(TensorRemove(sum(LinkTrace)))/vol/Nd/Nc;
+
+  ///////////////////////////////////////////////////
+  // Plaquette from covariant shifts
+  ///////////////////////////////////////////////////
+  ComplexField Plaq(grid);
+  Plaq = Zero();
+  for(int mu=1;mu<Nd;mu++){
+  for(int nu=0;nu<mu;nu++){
+    Plaq = Plaq + trace(PeriodicBC::CovShiftForward(U[mu],mu,U[nu])
+		    *adj(PeriodicBC::CovShiftForward(U[nu],nu,U[mu])));
+  }}
+  m.plaq_shift = real(TensorRemove(sum(Plaq)))/vol/6.0/Nc;
+  m.plaq_loops = WL::avgPlaquette(Umu);
+
+  ///////////////////////////////////////////////////
+  // 2x1 and 1x2 rectangles from covariant shifts
+  ///////////////////////////////////////////////////
+  ComplexField Rect(grid);
+  Rect = Zero();
+  for(int mu=1;mu<Nd;mu++){
+  for(int nu=0;nu<mu;nu++){
+    Rect = Rect + trace(
+	   PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[mu],mu,U[nu]))
+       *adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[mu],mu,U[mu]))) );
+    Rect = Rect + trace(
+	   PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[nu],nu,U[nu]))
+       *adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) );
+  }}
+  m.rect_shift = real(TensorRemove(sum(Rect)))/vol/12.0/Nc;
+  m.rect_loops = WL::avgRectangle(Umu);
+
+  ///////////////////////////////////////////////////
+  // Plaquette through the staples
+  ///////////////////////////////////////////////////
+  {
+    RealD stap = 0.0;
+    LinkField    staple(grid);
+    ComplexField stap_tr(grid);
+    for(int mu=0;mu<Nd;mu++){
+      WL::Staple(staple,Umu,mu);
+      stap_tr = trace(U[mu]*staple);
+      stap += real(TensorRemove(sum(stap_tr)));
+    }
+    m.plaq_staple = stap/vol/6.0/Nc/4.0;
+  }
+
+  ///////////////////////////////////////////////////
+  // Coarsened plaquette; blockSum must conserve the sum
+  ///////////////////////////////////////////////////
+  {
+    ComplexField cPlaq(coarse);
+    blockSum(cPlaq,Plaq);
+    m.coarse_plaq = real(TensorRemove(sum(cPlaq)))/vol/6.0/Nc;
+  }
+}
+
+template<class Gimpl>
+void Report(std::string name,Measurements<Gimpl> &m)
+{
+  std::cout << GridLogMessage << name << ": plaquette shifts " << m.plaq_shift
+            << " loops " << m.plaq_loops << " staples " << m.plaq_staple << std::endl;
+  std::cout << GridLogMessage << name << ": rectangle shifts " << m.rect_shift
+            << " loops " << m.rect_loops << std::endl;
+  std::cout << GridLogMessage << name << ": link trace " << m.link
+            << " coarsened plaquette " << m.coarse_plaq << std::endl;
+}
+
+template<class Gimpl>
+void SelfConsistent(Measurements<Gimpl> &m)
+{
+  GRID_ASSERT( fabs(m.plaq_shift  - m.plaq_loops) < 1.0e-8 );
+  GRID_ASSERT( fabs(m.plaq_staple - m.plaq_loops) < 1.0e-8 );
+  GRID_ASSERT( fabs(m.rect_shift  - m.rect_loops) < 1.0e-8 );
+  GRID_ASSERT( fabs(m.coarse_plaq - m.plaq_shift) < 1.0e-8 );
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate lsimd({1,1,1,1});
+  Coordinate mpi   = GridDefaultMpi();
+
+  Coordinate clatt(Nd);
+  for(int d=0;d<Nd;d++){
+    GRID_ASSERT( (latt[d]%2) == 0 );
+    clatt[d] = latt[d]/2;
+  }
+
+  GridCartesian vGrid(latt,vsimd,mpi);
+  GridCartesian lGrid(latt,lsimd,mpi);
+  GridCartesian vCoarse(clatt,vsimd,mpi);
+  GridCartesian lCoarse(clatt,lsimd,mpi);
+
+  std::cout << GridLogMessage << "vectorised    Nsimd = " << vGrid.Nsimd() << std::endl;
+  std::cout << GridLogMessage << "lexicographic Nsimd = " << lGrid.Nsimd() << std::endl;
+  GRID_ASSERT( lGrid.Nsimd() == 1 );
+
+  GridParallelRNG pRNG(&vGrid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  LatticeGaugeFieldD Umu(&vGrid);
+  SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  std::string file("./ckpoint_rectplaq_lex.4000");
+  NerscIO::writeConfiguration(Umu,file,0,0);
+
+  lexLatticeGaugeFieldD Ulex(&lGrid);
+  FieldMetaData header;
+  NerscIO::readConfiguration(Ulex,header,file);
+
+  Measurements<PeriodicGimplD>    vm;
+  Measurements<lexPeriodicGimplD> lm;
+
+  Measure<PeriodicGimplD>   (Umu ,&vCoarse,vm);
+  Measure<lexPeriodicGimplD>(Ulex,&lCoarse,lm);
+
+  Report("simd",vm);
+  Report("lex ",lm);
+
+  SelfConsistent(vm);
+  SelfConsistent(lm);
+
+  GRID_ASSERT( fabs(vm.plaq_shift  - lm.plaq_shift ) < tol );
+  GRID_ASSERT( fabs(vm.plaq_loops  - lm.plaq_loops ) < tol );
+  GRID_ASSERT( fabs(vm.plaq_staple - lm.plaq_staple) < tol );
+  GRID_ASSERT( fabs(vm.rect_shift  - lm.rect_shift ) < tol );
+  GRID_ASSERT( fabs(vm.rect_loops  - lm.rect_loops ) < tol );
+  GRID_ASSERT( fabs(vm.link        - lm.link       ) < tol );
+  GRID_ASSERT( fabs(vm.coarse_plaq - lm.coarse_plaq) < tol );
+
+  std::cout << GridLogMessage << "Test_RectPlaq_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_cshift_lex.cc

@@ -0,0 +1,175 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_cshift_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Lexicographic (Nsimd()==1) lattice types: LatticeCoordinate, Cshift in all
+// directions and shifts, and predicated where().
+//
+// Cshift is checked by an identity rather than a reference loop: with
+//    d = Cshift(coor_mu,mu,shift) - coor_mu
+// every site has d = shift (no wrap) or d = shift - L (wrapped), so
+//    (d - shift)*(d - shift + L) == 0
+// everywhere.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-5;
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  Coordinate latt   = GridDefaultLatt();
+  Coordinate simd({1,1,1,1});          // lexicographic: no SIMD
+  Coordinate mpi    = GridDefaultMpi();
+
+  GridCartesian grid(latt,simd,mpi);
+
+  std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << grid.Nsimd() << std::endl;
+  GRID_ASSERT( grid.Nsimd() == 1 );
+
+  
+
+  //////////////////////////////////////////////////
+  // Cshift against the wrap identity
+  //////////////////////////////////////////////////
+  for(int mu=0;mu<Nd;mu++){
+
+    lexLatticeComplexD coor(&grid);
+    LatticeCoordinate(coor,mu);
+
+    RealD L = latt[mu];
+
+    for(int shift=0;shift<latt[mu];shift++){
+
+      lexLatticeComplexD shifted(&grid);
+      shifted = Cshift(coor,mu,shift);
+
+      lexLatticeComplexD d(&grid);
+      d = shifted - coor;
+
+      lexLatticeComplexD p(&grid);
+      p = (d - ComplexD(shift)) * (d - ComplexD(shift) + ComplexD(L));
+
+      RealD n = norm2(p);
+      if ( n > tol ) {
+        std::cout << GridLogMessage << "FAIL mu " << mu << " shift " << shift
+                  << " residual " << n << std::endl;
+      }
+      GRID_ASSERT( n < tol );
+    }
+    std::cout << GridLogMessage << "Cshift mu = " << mu << " all shifts pass" << std::endl;
+  }
+
+  //////////////////////////////////////////////////
+  // Predicated where(): zero the second half in time
+  //////////////////////////////////////////////////
+  {
+    int    Tdir = Nd-1;
+    RealD  T    = latt[Tdir];
+    RealD  vol  = grid.gSites();
+
+    lexLatticeInteger  tcoor(&grid);
+    LatticeCoordinate(tcoor,Tdir);
+
+    lexLatticeComplexD f(&grid);   f  = 2.0;
+    lexLatticeComplexD zz(&grid);  zz = Zero();
+    lexLatticeComplexD g(&grid);
+
+    g = where( tcoor < Integer(T/2) , f , zz );
+
+    RealD n_g = norm2(g);
+    RealD expect = 4.0 * vol / 2.0;
+
+    std::cout << GridLogMessage << "where(): norm2 = " << n_g
+              << " expect " << expect << std::endl;
+    DumpSliceNorm("where() slice norm",g,Tdir);
+
+    GRID_ASSERT( fabs(n_g - expect) < tol*expect );
+  }
+
+  //////////////////////////////////////////////////
+  // Cshift against sliceSum of a slice-zeroed random field.
+  // Exhaustive over direction, zeroed slice and shift:
+  //    sliceSum(Cshift(f,mu,s))[t] == sliceSum(f)[(t+s)%L]
+  //////////////////////////////////////////////////
+  {
+    GridParallelRNG pRNG(&grid);
+    pRNG.SeedFixedIntegers(std::vector<int>({7,8,9,10}));
+
+    typedef lexLatticeComplexD::vector_object::scalar_object sobj;
+
+    lexLatticeComplexD rnd(&grid);  random(pRNG,rnd);
+    lexLatticeComplexD zz(&grid);   zz = Zero();
+
+    for(int mu=0;mu<Nd;mu++){
+
+      int L = latt[mu];
+
+      lexLatticeInteger coor(&grid);
+      LatticeCoordinate(coor,mu);
+
+      for(int t0=0;t0<L;t0++){
+
+        lexLatticeComplexD f(&grid);
+        f = where( coor == Integer(t0) , zz , rnd );
+
+        std::vector<sobj> ref;
+        sliceSum(f,ref,mu);
+
+        for(int shift=0;shift<L;shift++){
+
+          lexLatticeComplexD g(&grid);
+          g = Cshift(f,mu,shift);
+
+          std::vector<sobj> res;
+          sliceSum(g,res,mu);
+
+          for(int t=0;t<L;t++){
+            ComplexD got    = TensorRemove(res[t]);
+            ComplexD expect = TensorRemove(ref[(t+shift)%L]);
+            if ( abs(got-expect) > tol*(abs(expect)+1.0) ) {
+              std::cout << GridLogMessage << "FAIL mu " << mu << " t0 " << t0
+                        << " shift " << shift << " t " << t
+                        << " got " << got << " expect " << expect << std::endl;
+            }
+            GRID_ASSERT( abs(got-expect) < tol*(abs(expect)+1.0) );
+          }
+        }
+      }
+      std::cout << GridLogMessage << "sliceSum/Cshift consistency mu = " << mu
+                << " all zeroed slices and shifts pass" << std::endl;
+    }
+  }
+
+  std::cout << GridLogMessage << "Test_cshift_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_dwf_lex.cc

@@ -0,0 +1,196 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_dwf_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Five dimensional fermions on a lexicographic (Nsimd()==1) lattice.
+//
+// Domain wall and Mobius, driven from one gauge field and one source in both
+// charts and compared elementwise through the layout interchange.  Mobius uses
+// b=1.5 c=0.5, the production choice.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+////////////////////////////////////////////////////////////////////////
+// Layout interchange; both lattices share the same scalar_object.
+////////////////////////////////////////////////////////////////////////
+template<class vobjOut,class vobjIn>
+void transfer(Lattice<vobjOut> &out,const Lattice<vobjIn> &in)
+{
+  typedef typename vobjIn::scalar_object sobj;
+  static_assert(std::is_same<sobj,typename vobjOut::scalar_object>::value,
+		"transfer: lattices must share scalar_object");
+
+  GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites());
+
+  std::vector<sobj> buf;
+  unvectorizeToLexOrdArray(buf,in);
+  vectorizeFromLexOrdArray(buf,out);
+}
+
+////////////////////////////////////////////////////////////////////////
+// Compare one operator application between the two charts
+////////////////////////////////////////////////////////////////////////
+template<class vOp,class lexOp>
+void Compare(std::string name,int dag,
+	     vOp &vD,   typename vOp::FermionField   &vsrc,
+	     lexOp &lD, typename lexOp::FermionField &lsrc)
+{
+  typename vOp::FermionField     vres(vsrc.Grid());
+  typename lexOp::FermionField   lres(lsrc.Grid());
+  typename lexOp::FermionField   vres_lex(lsrc.Grid());
+  typename lexOp::FermionField   err(lsrc.Grid());
+
+  if ( dag == DaggerNo ) { vD.M   (vsrc,vres); lD.M   (lsrc,lres); }
+  else                   { vD.Mdag(vsrc,vres); lD.Mdag(lsrc,lres); }
+
+  transfer(vres_lex,vres);
+  err = vres_lex - lres;
+  RealD n = norm2(err)/norm2(lres);
+  std::cout << GridLogMessage << name << " simd vs lex: relative " << n
+	    << "  |simd|^2 " << norm2(vres) << "  |lex|^2 " << norm2(lres) << std::endl;
+  GRID_ASSERT( n < tol );
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls = 8;
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate lsimd({1,1,1,1});
+  Coordinate mpi   = GridDefaultMpi();
+
+  GridCartesian         *vUGrid   = new GridCartesian(latt,vsimd,mpi);
+  GridRedBlackCartesian *vUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(vUGrid);
+  GridCartesian         *vFGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,vUGrid);
+  GridRedBlackCartesian *vFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,vUGrid);
+
+  GridCartesian         *lUGrid   = new GridCartesian(latt,lsimd,mpi);
+  GridRedBlackCartesian *lUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(lUGrid);
+  GridCartesian         *lFGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,lUGrid);
+  GridRedBlackCartesian *lFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,lUGrid);
+
+  std::cout << GridLogMessage << "Ls = " << Ls
+	    << "  vectorised Nsimd(4d) = " << vUGrid->Nsimd()
+	    << " Nsimd(5d) = " << vFGrid->Nsimd() << std::endl;
+  std::cout << GridLogMessage << "     lexicographic Nsimd(4d) = " << lUGrid->Nsimd()
+	    << " Nsimd(5d) = " << lFGrid->Nsimd() << std::endl;
+  GRID_ASSERT( lUGrid->Nsimd() == 1 );
+  GRID_ASSERT( lFGrid->Nsimd() == 1 );
+
+  RealD mass = 0.05;
+  RealD M5   = 1.8;
+
+  GridParallelRNG pRNG4(vUGrid); pRNG4.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+  GridParallelRNG pRNG5(vFGrid); pRNG5.SeedFixedIntegers(std::vector<int>({5,6,7,8}));
+
+  //////////////////////////////////////////////////
+  // One gauge field and one source, both charts
+  //////////////////////////////////////////////////
+  LatticeGaugeFieldD Umu(vUGrid); SU<Nc>::HotConfiguration(pRNG4,Umu);
+  lexLatticeGaugeFieldD Ulex(lUGrid); transfer(Ulex,Umu);
+
+  typedef DomainWallFermion<WilsonImplD>    vDwf;
+  typedef DomainWallFermion<lexWilsonImplD> lexDwf;
+  typedef MobiusFermion<WilsonImplD>        vMob;
+  typedef MobiusFermion<lexWilsonImplD>     lexMob;
+
+  typename vDwf::FermionField   src(vFGrid);   random(pRNG5,src);
+  typename vDwf::FermionField   phi(vFGrid);   random(pRNG5,phi);
+  typename lexDwf::FermionField srclex(lFGrid); transfer(srclex,src);
+  typename lexDwf::FermionField philex(lFGrid); transfer(philex,phi);
+
+  //////////////////////////////////////////////////
+  // Domain wall
+  //////////////////////////////////////////////////
+  {
+    vDwf   Dv (Umu ,*vFGrid,*vFrbGrid,*vUGrid,*vUrbGrid,mass,M5);
+    lexDwf Dl (Ulex,*lFGrid,*lFrbGrid,*lUGrid,*lUrbGrid,mass,M5);
+
+    Compare("DWF M   ",DaggerNo ,Dv,src,Dl,srclex);
+    Compare("DWF Mdag",DaggerYes,Dv,src,Dl,srclex);
+
+    // Adjoint identity within the lexicographic chart
+    typename lexDwf::FermionField Mphi(lFGrid),Mdagchi(lFGrid);
+    Dl.M   (srclex,Mphi);
+    Dl.Mdag(philex,Mdagchi);
+    ComplexD lhs = innerProduct(philex,Mphi);
+    ComplexD rhs = innerProduct(Mdagchi,srclex);
+    RealD n = abs(lhs-rhs)/abs(lhs);
+    std::cout << GridLogMessage << "DWF lex <phi|M src> = " << lhs
+	      << "  <Mdag phi|src> = " << rhs << "  relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+
+    // Checkerboarded hopping term
+    typename vDwf::FermionField   vo(vFrbGrid),ve(vFrbGrid);
+    typename lexDwf::FermionField lo(lFrbGrid),le(lFrbGrid);
+    vo.Checkerboard()=Odd; ve.Checkerboard()=Even;
+    lo.Checkerboard()=Odd; le.Checkerboard()=Even;
+    pickCheckerboard(Odd,vo,src);
+    pickCheckerboard(Odd,lo,srclex);
+    Dv.Meooe(vo,ve);
+    Dl.Meooe(lo,le);
+    RealD vn = norm2(ve), ln = norm2(le);
+    std::cout << GridLogMessage << "DWF Meooe |simd|^2 " << vn << " |lex|^2 " << ln
+	      << " relative " << fabs(vn-ln)/vn << std::endl;
+    GRID_ASSERT( fabs(vn-ln)/vn < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Mobius, production coefficients
+  //////////////////////////////////////////////////
+  {
+    RealD b=1.5, c=0.5;
+    vMob   Dv (Umu ,*vFGrid,*vFrbGrid,*vUGrid,*vUrbGrid,mass,M5,b,c);
+    lexMob Dl (Ulex,*lFGrid,*lFrbGrid,*lUGrid,*lUrbGrid,mass,M5,b,c);
+
+    Compare("Mobius M   ",DaggerNo ,Dv,src,Dl,srclex);
+    Compare("Mobius Mdag",DaggerYes,Dv,src,Dl,srclex);
+
+    // MooeeInv is the Ls-direction solve; check it inverts Mooee in the lex chart
+    typename lexMob::FermionField lo(lFrbGrid),t1(lFrbGrid),t2(lFrbGrid),e(lFrbGrid);
+    lo.Checkerboard()=Odd; t1.Checkerboard()=Odd; t2.Checkerboard()=Odd; e.Checkerboard()=Odd;
+    pickCheckerboard(Odd,lo,srclex);
+    Dl.Mooee(lo,t1);
+    Dl.MooeeInv(t1,t2);
+    e = t2 - lo;
+    RealD n = norm2(e)/norm2(lo);
+    std::cout << GridLogMessage << "Mobius lex MooeeInv(Mooee) - 1: relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  std::cout << GridLogMessage << "Test_dwf_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_innerproduct_norm_lex.cc

@@ -0,0 +1,132 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_innerproduct_norm_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Lexicographic (Nsimd()==1) lattice types: linear combination, expression
+// templates, norm2, innerProduct, innerProductNorm and RNG fill.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-5;
+
+template<class Field>
+void BasicChecks(GridCartesian *grid,GridParallelRNG &pRNG,std::string precision)
+{
+  RealD vol = grid->gSites();
+
+  Field a(grid);  a = 1.0;
+  Field b(grid);  b = 2.0;
+  Field c(grid);  c = 3.0;
+  Field d(grid);
+
+  ///////////////////////////////////////////////////
+  // Linear combination through expression templates
+  ///////////////////////////////////////////////////
+  d = a + b - c;
+  RealD n_zero = norm2(d);
+
+  RealD n_a = norm2(a);
+  RealD n_b = norm2(b);
+
+  std::cout << GridLogMessage << precision << ": norm2(1+2-3) = " << n_zero << std::endl;
+  std::cout << GridLogMessage << precision << ": norm2(1) = " << n_a << " expect " << vol << std::endl;
+  std::cout << GridLogMessage << precision << ": norm2(2) = " << n_b << " expect " << 4.0*vol << std::endl;
+
+  GRID_ASSERT( n_zero          < tol );
+  GRID_ASSERT( fabs(n_a -     vol) < tol*vol );
+  GRID_ASSERT( fabs(n_b - 4.0*vol) < tol*vol );
+
+  ///////////////////////////////////////////////////
+  // innerProduct of constant fields
+  ///////////////////////////////////////////////////
+  ComplexD ip = innerProduct(a,b);
+  std::cout << GridLogMessage << precision << ": innerProduct(1,2) = " << ip
+            << " expect (" << 2.0*vol << ",0)" << std::endl;
+  GRID_ASSERT( fabs(real(ip) - 2.0*vol) < tol*vol );
+  GRID_ASSERT( fabs(imag(ip))           < tol*vol );
+
+  ///////////////////////////////////////////////////
+  // RNG fill; norm2(r+r) == 4 norm2(r)
+  ///////////////////////////////////////////////////
+  Field r(grid); random(pRNG,r);
+  Field s(grid); random(pRNG,s);
+
+  RealD n_r  = norm2(r);
+  Field rr(grid); rr = r + r;
+  RealD n_rr = norm2(rr);
+
+  std::cout << GridLogMessage << precision << ": norm2(r) = " << n_r
+            << " norm2(r+r) = " << n_rr << " ratio " << n_rr/n_r << std::endl;
+  GRID_ASSERT( n_r > tol );                          // RNG actually filled it
+  GRID_ASSERT( fabs(n_rr - 4.0*n_r) < tol*n_rr );
+
+  ///////////////////////////////////////////////////
+  // Fused innerProductNorm against separate calls
+  ///////////////////////////////////////////////////
+  ComplexD ip_ref = innerProduct(r,s);
+  RealD    n2_ref = norm2(r);
+
+  ComplexD ip_res;
+  RealD    n2_res;
+  innerProductNorm(ip_res,n2_res,r,s);
+
+  std::cout << GridLogMessage << precision << ": innerProductNorm ip diff "
+            << abs(ip_ref-ip_res) << " norm2 diff " << fabs(n2_ref-n2_res) << std::endl;
+  GRID_ASSERT( abs(ip_ref-ip_res)   < tol*abs(ip_ref) + tol );
+  GRID_ASSERT( fabs(n2_ref-n2_res)  < tol*n2_ref );
+
+  std::cout << GridLogMessage << precision << ": all checks passed" << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  Coordinate latt   = GridDefaultLatt();
+  Coordinate simd({1,1,1,1});          // lexicographic: no SIMD
+  Coordinate mpi    = GridDefaultMpi();
+
+  GridCartesian grid(latt,simd,mpi);
+
+  std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << grid.Nsimd() << std::endl;
+  GRID_ASSERT( grid.Nsimd() == 1 );
+
+  GridParallelRNG pRNG(&grid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  
+
+  BasicChecks<lexLatticeComplexD>(&grid,pRNG,"Double");
+  BasicChecks<lexLatticeComplexF>(&grid,pRNG,"Single");
+
+  std::cout << GridLogMessage << "Test_innerproduct_norm_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_io.cc

@@ -0,0 +1,194 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_io.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Interchange between vectorised and lexicographic (Nsimd()==1) lattices.
+//
+//   - transfer() in both directions via an unvectorise/vectorise pair
+//   - RNG seeded identically on both layouts produces identical fields
+//   - binary I/O written from one layout and read into the other
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+////////////////////////////////////////////////////////////////////////
+// Layout interchange.  Both lattices share the same scalar_object, so a
+// single lexicographic array mediates; works in either direction.
+////////////////////////////////////////////////////////////////////////
+template<class vobjOut,class vobjIn>
+void transfer(Lattice<vobjOut> &out,const Lattice<vobjIn> &in)
+{
+  typedef typename vobjIn::scalar_object sobj;
+  static_assert(std::is_same<sobj,typename vobjOut::scalar_object>::value,
+		"transfer: lattices must share scalar_object");
+
+  GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites());
+
+  std::vector<sobj> buf;
+  unvectorizeToLexOrdArray(buf,in);
+  vectorizeFromLexOrdArray(buf,out);
+}
+
+////////////////////////////////////////////////////////////////////////
+// Binary write/read of a single lattice object
+////////////////////////////////////////////////////////////////////////
+template<class Field>
+void writeField(const Field &f,std::string file)
+{
+  typedef typename Field::vector_object vobj;
+  typedef typename Field::scalar_object sobj;
+
+  BinarySimpleMunger<sobj,sobj> munge;
+  std::string format = getFormatString<vobj>();
+  uint64_t offset = 0;
+  uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+
+  BinaryIO::writeLatticeObject<vobj,sobj>(const_cast<Field &>(f),file,munge,offset,format,
+					  nersc_csum,scidac_csuma,scidac_csumb);
+}
+
+template<class Field>
+void readField(Field &f,std::string file)
+{
+  typedef typename Field::vector_object vobj;
+  typedef typename Field::scalar_object sobj;
+
+  BinarySimpleMunger<sobj,sobj> munge;
+  std::string format = getFormatString<vobj>();
+  uint64_t offset = 0;
+  uint32_t nersc_csum,scidac_csuma,scidac_csumb;
+
+  BinaryIO::readLatticeObject<vobj,sobj>(f,file,munge,offset,format,
+					 nersc_csum,scidac_csuma,scidac_csumb);
+}
+
+int main(int argc, char** argv)
+{
+  Grid_init(&argc, &argv);
+
+  Coordinate latt = GridDefaultLatt();
+  Coordinate mpi  = GridDefaultMpi();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate ssimd({1,1,1,1});
+
+  GridCartesian vGrid(latt,vsimd,mpi);
+  GridCartesian sGrid(latt,ssimd,mpi);
+
+  std::cout << GridLogMessage << "vectorised   grid Nsimd = " << vGrid.Nsimd() << std::endl;
+  std::cout << GridLogMessage << "lexicographic grid Nsimd = " << sGrid.Nsimd() << std::endl;
+  GRID_ASSERT( sGrid.Nsimd() == 1 );
+
+  typedef LatticeComplexD    vField;
+  typedef lexLatticeComplexD sField;
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG vRNG(&vGrid); vRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG sRNG(&sGrid); sRNG.SeedFixedIntegers(seeds);
+
+  vField v(&vGrid);  random(vRNG,v);
+  sField s(&sGrid);  random(sRNG,s);
+
+  RealD nv = norm2(v);
+  RealD ns = norm2(s);
+  std::cout << GridLogMessage << "norm2 vectorised   = " << nv << std::endl;
+  std::cout << GridLogMessage << "norm2 lexicographic = " << ns << std::endl;
+  GRID_ASSERT( nv > tol );
+
+  ////////////////////////////////////////////////////////
+  // A: transfer round trip v -> s -> v
+  ////////////////////////////////////////////////////////
+  {
+    sField st(&sGrid);
+    vField vt(&vGrid);
+
+    transfer(st,v);
+    transfer(vt,st);
+
+    vField d(&vGrid); d = vt - v;
+    RealD n = norm2(d);
+    std::cout << GridLogMessage << "A: round trip v->s->v  norm2(diff) = " << n << std::endl;
+    GRID_ASSERT( n < tol );
+
+    // and the transferred copy must carry the same norm
+    std::cout << GridLogMessage << "A: norm2 transferred = " << norm2(st) << std::endl;
+    GRID_ASSERT( fabs(norm2(st) - nv) < tol*nv );
+  }
+
+  ////////////////////////////////////////////////////////
+  // B: identically seeded RNGs agree across layouts
+  ////////////////////////////////////////////////////////
+  {
+    vField vs(&vGrid);
+    transfer(vs,s);
+
+    vField d(&vGrid); d = vs - v;
+    RealD n = norm2(d);
+    std::cout << GridLogMessage << "B: same seed, both layouts  norm2(diff) = " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  ////////////////////////////////////////////////////////
+  // C: write vectorised, read lexicographic
+  ////////////////////////////////////////////////////////
+  {
+    sField sref(&sGrid);  transfer(sref,v);   // what the file should contain
+    sField sin(&sGrid);   sin = Zero();
+
+    writeField(v,"nonsimd_io_v.bin");
+    readField(sin,"nonsimd_io_v.bin");
+
+    sField d(&sGrid); d = sin - sref;
+    RealD n = norm2(d);
+    std::cout << GridLogMessage << "C: write simd / read lexicographic  norm2(diff) = " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  ////////////////////////////////////////////////////////
+  // D: write lexicographic, read vectorised
+  ////////////////////////////////////////////////////////
+  {
+    sField sout(&sGrid);  transfer(sout,v);
+    vField vin(&vGrid);   vin = Zero();
+
+    writeField(sout,"nonsimd_io_s.bin");
+    readField(vin,"nonsimd_io_s.bin");
+
+    vField d(&vGrid); d = vin - v;
+    RealD n = norm2(d);
+    std::cout << GridLogMessage << "D: write lexicographic / read simd  norm2(diff) = " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  std::cout << GridLogMessage << "Test_io: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_staggered_lex.cc

@@ -0,0 +1,246 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_staggered_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Naive staggered fermion operator on a lexicographic (Nsimd()==1) lattice.
+//
+// Two independent checks:
+//   - lexicographic Dhop against a covariant-shift reference built in the
+//     same chart
+//   - both charts driven from the same gauge field and source, compared
+//     elementwise through the layout interchange
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+typedef NaiveStaggeredFermion<StaggeredImplD>    vStagOp;
+typedef NaiveStaggeredFermion<lexStaggeredImplD> lexStagOp;
+
+////////////////////////////////////////////////////////////////////////
+// Layout interchange; both lattices share the same scalar_object.
+////////////////////////////////////////////////////////////////////////
+template<class vobjOut,class vobjIn>
+void transfer(Lattice<vobjOut> &out,const Lattice<vobjIn> &in)
+{
+  typedef typename vobjIn::scalar_object sobj;
+  static_assert(std::is_same<sobj,typename vobjOut::scalar_object>::value,
+		"transfer: lattices must share scalar_object");
+
+  GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites());
+
+  std::vector<sobj> buf;
+  unvectorizeToLexOrdArray(buf,in);
+  vectorizeFromLexOrdArray(buf,out);
+}
+
+////////////////////////////////////////////////////////////////////////
+// Dhop from covariant shifts, in whichever chart Gimpl names
+////////////////////////////////////////////////////////////////////////
+template<class Impl>
+void ReferenceDhop(typename Impl::FermionField &ref,
+		   const typename Impl::GaugeField &Umu,
+		   const typename Impl::FermionField &src,
+		   RealD c1,RealD u0)
+{
+  typedef typename Impl::GaugeLinkField LinkField;
+  typedef typename Impl::ComplexField   ComplexField;
+  typedef typename Impl::FermionField   FermionField;
+
+  GridBase *grid = src.Grid();
+  RealD c1tad = 0.5*c1/u0;
+
+  std::vector<LinkField> U(Nd,grid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+
+  typedef Lattice<iScalar<typename GridTypeMapper<typename Impl::Simd>::Integerified> > IntField;
+
+  IntField x(grid); LatticeCoordinate(x,0);
+  IntField y(grid); LatticeCoordinate(y,1);
+  IntField z(grid); LatticeCoordinate(z,2);
+  IntField lin_z(grid); lin_z = x+y;
+  IntField lin_t(grid); lin_t = x+y+z;
+
+  FermionField tmp(grid);
+  ref = Zero();
+
+  for(int mu=0;mu<Nd;mu++){
+
+    ComplexField phases(grid); phases=1.0;
+    if ( mu == 1 ) phases = where( mod(x    ,2)==(Integer)0, phases,-phases);
+    if ( mu == 2 ) phases = where( mod(lin_z,2)==(Integer)0, phases,-phases);
+    if ( mu == 3 ) phases = where( mod(lin_t,2)==(Integer)0, phases,-phases);
+
+    tmp = PeriodicBC::CovShiftForward(U[mu],mu,src);
+    ref = ref + c1tad*tmp*phases;
+
+    tmp = PeriodicBC::CovShiftBackward(U[mu],mu,src);
+    ref = ref - c1tad*tmp*phases;
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate lsimd({1,1,1,1});
+  Coordinate mpi   = GridDefaultMpi();
+
+  GridCartesian         vGrid(latt,vsimd,mpi);
+  GridRedBlackCartesian vRBGrid(&vGrid);
+  GridCartesian         lGrid(latt,lsimd,mpi);
+  GridRedBlackCartesian lRBGrid(&lGrid);
+
+  std::cout << GridLogMessage << "vectorised    Nsimd = " << vGrid.Nsimd() << std::endl;
+  std::cout << GridLogMessage << "lexicographic Nsimd = " << lGrid.Nsimd() << std::endl;
+  GRID_ASSERT( lGrid.Nsimd() == 1 );
+
+  RealD mass = 0.1;
+  RealD c1   = 9.0/8.0;
+  RealD u0   = 1.0;
+
+  GridParallelRNG pRNG(&vGrid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  //////////////////////////////////////////////////
+  // One gauge field and one source, both charts
+  //////////////////////////////////////////////////
+  typename vStagOp::GaugeField   Umu(&vGrid);   SU<Nc>::HotConfiguration(pRNG,Umu);
+  typename vStagOp::FermionField src(&vGrid);   random(pRNG,src);
+  typename vStagOp::FermionField phi(&vGrid);   random(pRNG,phi);
+
+  typename lexStagOp::GaugeField   Ulex(&lGrid);  transfer(Ulex,Umu);
+  typename lexStagOp::FermionField srclex(&lGrid); transfer(srclex,src);
+  typename lexStagOp::FermionField philex(&lGrid); transfer(philex,phi);
+
+  typename vStagOp::ImplParams   vparams;
+  typename lexStagOp::ImplParams lparams;
+
+  vStagOp   Ds   (Umu ,vGrid,vRBGrid,mass,c1,u0,vparams);
+  lexStagOp Dslex(Ulex,lGrid,lRBGrid,mass,c1,u0,lparams);
+
+  //////////////////////////////////////////////////
+  // Lexicographic Dhop against covariant shifts
+  //////////////////////////////////////////////////
+  {
+    typename lexStagOp::FermionField ref(&lGrid);
+    typename lexStagOp::FermionField res(&lGrid);
+    typename lexStagOp::FermionField err(&lGrid);
+
+    ReferenceDhop<lexStaggeredImplD>(ref,Ulex,srclex,c1,u0);
+    Dslex.Dhop(srclex,res,DaggerNo);
+
+    err = res - ref;
+    RealD n = norm2(err)/norm2(ref);
+    std::cout << GridLogMessage << "lex Dhop vs covariant shifts: relative " << n
+	      << " (|Dhop|^2 " << norm2(res) << ")" << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Chart equivalence, elementwise
+  //////////////////////////////////////////////////
+  {
+    typename vStagOp::FermionField   vres(&vGrid);
+    typename lexStagOp::FermionField lres(&lGrid);
+    typename lexStagOp::FermionField vres_lex(&lGrid);
+    typename lexStagOp::FermionField err(&lGrid);
+
+    struct { const char *name; int dag; } ops[2] = { {"M",DaggerNo}, {"Mdag",DaggerYes} };
+
+    for(int o=0;o<2;o++){
+      if ( ops[o].dag == DaggerNo ) { Ds.M   (src,vres); Dslex.M   (srclex,lres); }
+      else                          { Ds.Mdag(src,vres); Dslex.Mdag(srclex,lres); }
+
+      transfer(vres_lex,vres);
+      err = vres_lex - lres;
+      RealD n = norm2(err)/norm2(lres);
+      std::cout << GridLogMessage << ops[o].name << " simd vs lex: relative " << n
+		<< "  |simd|^2 " << norm2(vres) << "  |lex|^2 " << norm2(lres) << std::endl;
+      GRID_ASSERT( n < tol );
+    }
+
+    Ds.Dhop(src,vres,DaggerNo);  Dslex.Dhop(srclex,lres,DaggerNo);
+    transfer(vres_lex,vres);
+    err = vres_lex - lres;
+    RealD n = norm2(err)/norm2(lres);
+    std::cout << GridLogMessage << "Dhop simd vs lex: relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Adjoint identity within the lexicographic chart
+  //////////////////////////////////////////////////
+  {
+    typename lexStagOp::FermionField Mphi(&lGrid);
+    typename lexStagOp::FermionField Mdagchi(&lGrid);
+
+    Dslex.M   (srclex,Mphi);
+    Dslex.Mdag(philex,Mdagchi);
+
+    ComplexD lhs = innerProduct(philex,Mphi);
+    ComplexD rhs = innerProduct(Mdagchi,srclex);
+    RealD n = abs(lhs-rhs)/abs(lhs);
+    std::cout << GridLogMessage << "lex <phi|M src> = " << lhs
+	      << "  <Mdag phi|src> = " << rhs << "  relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Checkerboarded hopping term, both charts
+  //////////////////////////////////////////////////
+  {
+    typename vStagOp::FermionField   vsrc_o(&vRBGrid);   vsrc_o.Checkerboard()=Odd;
+    typename vStagOp::FermionField   vres_e(&vRBGrid);   vres_e.Checkerboard()=Even;
+    typename lexStagOp::FermionField lsrc_o(&lRBGrid);   lsrc_o.Checkerboard()=Odd;
+    typename lexStagOp::FermionField lres_e(&lRBGrid);   lres_e.Checkerboard()=Even;
+
+    pickCheckerboard(Odd,vsrc_o,src);
+    pickCheckerboard(Odd,lsrc_o,srclex);
+
+    Ds.Meooe   (vsrc_o,vres_e);
+    Dslex.Meooe(lsrc_o,lres_e);
+
+    RealD vn = norm2(vres_e);
+    RealD ln = norm2(lres_e);
+    RealD n  = fabs(vn-ln)/vn;
+    std::cout << GridLogMessage << "Meooe |simd|^2 " << vn << "  |lex|^2 " << ln
+	      << "  relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  std::cout << GridLogMessage << "Test_staggered_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_stencil_lex.cc

@@ -0,0 +1,122 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_stencil_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// CartesianStencil on a lexicographic (Nsimd()==1) lattice, checked against
+// Cshift for every direction and displacement of either sign.
+//
+// On a unit simd grid no stencil entry may request a permute; the Grid_simd1
+// permute is an assert, so a request would abort rather than pass silently.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+int main(int argc, char ** argv)
+{
+  Grid_init(&argc, &argv);
+
+  typedef lexLatticeComplexD Field;
+  typedef Field::vector_object vobj;
+
+  Coordinate latt = GridDefaultLatt();
+  Coordinate simd({1,1,1,1});
+  Coordinate mpi  = GridDefaultMpi();
+
+  GridCartesian Fine(latt,simd,mpi);
+
+  std::cout << GridLogMessage << "Lexicographic grid, Nsimd = " << Fine.Nsimd() << std::endl;
+  GRID_ASSERT( Fine.Nsimd() == 1 );
+
+  GridParallelRNG fRNG(&Fine);
+  fRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  Field Foo(&Fine);   random(fRNG,Foo);
+  Field Bar(&Fine);
+  Field Check(&Fine);
+  Field Diff(&Fine);
+
+  typedef CartesianStencil<vobj,vobj,SimpleStencilParams> Stencil;
+  SimpleStencilParams p;
+
+  for(int dir=0;dir<Nd;dir++){
+
+    int L = Fine._fdimensions[dir];
+
+    for(int disp=-(L-1);disp<=L-1;disp++){
+
+      int npoint=1;
+      std::vector<int> directions(npoint,dir);
+      std::vector<int> displacements(npoint,disp);
+
+      Stencil myStencil(&Fine,npoint,0,directions,displacements,p);
+
+      SimpleCompressor<vobj> compress;
+      myStencil.HaloExchange(Foo,compress);
+
+      Bar = Cshift(Foo,dir,disp);
+
+      {
+	autoView( check , Check, AcceleratorWrite);
+	autoView( foo   , Foo,   AcceleratorRead);
+	autoView( st_v  , myStencil, AcceleratorRead);
+	auto CBp=myStencil.CommBuf();
+	accelerator_for(i,Check.Grid()->oSites(), 1, {
+
+	    int permute_type;
+	    StencilEntry *SE;
+	    SE = st_v.GetEntry(permute_type,0,i);
+
+	    if ( SE->_is_local && SE->_permute )
+	      permute(check[i],foo[SE->_offset],permute_type);
+	    else if (SE->_is_local)
+	      check[i] = foo[SE->_offset];
+	    else
+	      check[i] = CBp[SE->_offset];
+	});
+      }
+
+      Diff = Check-Bar;
+      RealD nrm = norm2(Diff);
+
+      if ( nrm > tol ) {
+	std::cout << GridLogMessage << "FAIL dir " << dir << " disp " << disp
+		  << " norm2(stencil-cshift) = " << nrm << std::endl;
+      }
+      GRID_ASSERT( nrm < tol );
+    }
+    std::cout << GridLogMessage << "Stencil == Cshift, dir " << dir
+	      << ", all displacements " << -(L-1) << " .. " << L-1 << std::endl;
+  }
+
+  std::cout << GridLogMessage << "Test_stencil_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/lexLattice/Test_wilson_lex.cc

@@ -0,0 +1,239 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/lexLattice/Test_wilson_lex.cc
+
+    Copyright (C) 2026
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//
+// Wilson fermion operator on a lexicographic (Nsimd()==1) lattice.
+//
+// Exercises the spin-projected half spinor path -- WilsonCompressor and
+// WilsonStencil -- which the staggered and stencil tests do not reach.
+//
+//   - lexicographic Dhop against a covariant-shift reference in the same chart
+//   - both charts from one gauge field and source, compared elementwise
+//   - gamma5 hermiticity and the even-odd hopping term in the lex chart
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const RealD tol = 1.0e-10;
+
+Gamma::Algebra Gmu [] = {
+  Gamma::Algebra::GammaX,
+  Gamma::Algebra::GammaY,
+  Gamma::Algebra::GammaZ,
+  Gamma::Algebra::GammaT
+};
+
+typedef WilsonFermion<WilsonImplD>    vWilsonOp;
+typedef WilsonFermion<lexWilsonImplD> lexWilsonOp;
+
+////////////////////////////////////////////////////////////////////////
+// Layout interchange; both lattices share the same scalar_object.
+////////////////////////////////////////////////////////////////////////
+template<class vobjOut,class vobjIn>
+void transfer(Lattice<vobjOut> &out,const Lattice<vobjIn> &in)
+{
+  typedef typename vobjIn::scalar_object sobj;
+  static_assert(std::is_same<sobj,typename vobjOut::scalar_object>::value,
+		"transfer: lattices must share scalar_object");
+
+  GRID_ASSERT(out.Grid()->gSites() == in.Grid()->gSites());
+
+  std::vector<sobj> buf;
+  unvectorizeToLexOrdArray(buf,in);
+  vectorizeFromLexOrdArray(buf,out);
+}
+
+////////////////////////////////////////////////////////////////////////
+// Wilson hopping term from covariant shifts, in whichever chart Impl names
+////////////////////////////////////////////////////////////////////////
+template<class Impl>
+void ReferenceDhop(typename Impl::FermionField &ref,
+		   const typename Impl::GaugeField &Umu,
+		   const typename Impl::FermionField &src)
+{
+  typedef typename Impl::GaugeLinkField LinkField;
+  typedef typename Impl::FermionField   FermionField;
+
+  GridBase *grid = src.Grid();
+
+  std::vector<LinkField> U(Nd,grid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+
+  FermionField tmp(grid);
+  ref = Zero();
+
+  for(int mu=0;mu<Nd;mu++){
+    tmp = U[mu]*Cshift(src,mu,1);
+    ref = ref + tmp - Gamma(Gmu[mu])*tmp;
+
+    tmp = adj(U[mu])*src;
+    tmp = Cshift(tmp,mu,-1);
+    ref = ref + tmp + Gamma(Gmu[mu])*tmp;
+  }
+  ref = -0.5*ref;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate vsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate lsimd({1,1,1,1});
+  Coordinate mpi   = GridDefaultMpi();
+
+  GridCartesian         vGrid(latt,vsimd,mpi);
+  GridRedBlackCartesian vRBGrid(&vGrid);
+  GridCartesian         lGrid(latt,lsimd,mpi);
+  GridRedBlackCartesian lRBGrid(&lGrid);
+
+  std::cout << GridLogMessage << "vectorised    Nsimd = " << vGrid.Nsimd() << std::endl;
+  std::cout << GridLogMessage << "lexicographic Nsimd = " << lGrid.Nsimd() << std::endl;
+  GRID_ASSERT( lGrid.Nsimd() == 1 );
+
+  RealD mass = 0.1;
+
+  GridParallelRNG pRNG(&vGrid);
+  pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  //////////////////////////////////////////////////
+  // One gauge field and one source, both charts
+  //////////////////////////////////////////////////
+  typename vWilsonOp::GaugeField   Umu(&vGrid);  SU<Nc>::HotConfiguration(pRNG,Umu);
+  typename vWilsonOp::FermionField src(&vGrid);  random(pRNG,src);
+  typename vWilsonOp::FermionField phi(&vGrid);  random(pRNG,phi);
+
+  typename lexWilsonOp::GaugeField   Ulex(&lGrid);   transfer(Ulex,Umu);
+  typename lexWilsonOp::FermionField srclex(&lGrid); transfer(srclex,src);
+  typename lexWilsonOp::FermionField philex(&lGrid); transfer(philex,phi);
+
+  vWilsonOp   Dw   (Umu ,vGrid,vRBGrid,mass);
+  lexWilsonOp Dwlex(Ulex,lGrid,lRBGrid,mass);
+
+  //////////////////////////////////////////////////
+  // Lexicographic Dhop against covariant shifts
+  //////////////////////////////////////////////////
+  {
+    typename lexWilsonOp::FermionField ref(&lGrid);
+    typename lexWilsonOp::FermionField res(&lGrid);
+    typename lexWilsonOp::FermionField err(&lGrid);
+
+    ReferenceDhop<lexWilsonImplD>(ref,Ulex,srclex);
+    Dwlex.Dhop(srclex,res,DaggerNo);
+
+    err = res - ref;
+    RealD n = norm2(err)/norm2(ref);
+    std::cout << GridLogMessage << "lex Dhop vs covariant shifts: relative " << n
+	      << " (|Dhop|^2 " << norm2(res) << ")" << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Chart equivalence, elementwise
+  //////////////////////////////////////////////////
+  {
+    typename vWilsonOp::FermionField   vres(&vGrid);
+    typename lexWilsonOp::FermionField lres(&lGrid);
+    typename lexWilsonOp::FermionField vres_lex(&lGrid);
+    typename lexWilsonOp::FermionField err(&lGrid);
+
+    Dw.M(src,vres);   Dwlex.M(srclex,lres);
+    transfer(vres_lex,vres);
+    err = vres_lex - lres;
+    RealD n = norm2(err)/norm2(lres);
+    std::cout << GridLogMessage << "M simd vs lex: relative " << n
+	      << "  |simd|^2 " << norm2(vres) << "  |lex|^2 " << norm2(lres) << std::endl;
+    GRID_ASSERT( n < tol );
+
+    Dw.Mdag(src,vres);   Dwlex.Mdag(srclex,lres);
+    transfer(vres_lex,vres);
+    err = vres_lex - lres;
+    n = norm2(err)/norm2(lres);
+    std::cout << GridLogMessage << "Mdag simd vs lex: relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+
+    Dw.Dhop(src,vres,DaggerNo);   Dwlex.Dhop(srclex,lres,DaggerNo);
+    transfer(vres_lex,vres);
+    err = vres_lex - lres;
+    n = norm2(err)/norm2(lres);
+    std::cout << GridLogMessage << "Dhop simd vs lex: relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // gamma5 hermiticity in the lexicographic chart
+  //   g5 M g5 = Mdag
+  //////////////////////////////////////////////////
+  {
+    typename lexWilsonOp::FermionField tmp(&lGrid);
+    typename lexWilsonOp::FermionField g5Mg5(&lGrid);
+    typename lexWilsonOp::FermionField Mdag(&lGrid);
+    typename lexWilsonOp::FermionField err(&lGrid);
+
+    tmp = Gamma(Gamma::Algebra::Gamma5)*srclex;
+    Dwlex.M(tmp,g5Mg5);
+    g5Mg5 = Gamma(Gamma::Algebra::Gamma5)*g5Mg5;
+
+    Dwlex.Mdag(srclex,Mdag);
+
+    err = g5Mg5 - Mdag;
+    RealD n = norm2(err)/norm2(Mdag);
+    std::cout << GridLogMessage << "lex g5 M g5 - Mdag: relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  //////////////////////////////////////////////////
+  // Checkerboarded hopping term, both charts
+  //////////////////////////////////////////////////
+  {
+    typename vWilsonOp::FermionField   vsrc_o(&vRBGrid);   vsrc_o.Checkerboard()=Odd;
+    typename vWilsonOp::FermionField   vres_e(&vRBGrid);   vres_e.Checkerboard()=Even;
+    typename lexWilsonOp::FermionField lsrc_o(&lRBGrid);   lsrc_o.Checkerboard()=Odd;
+    typename lexWilsonOp::FermionField lres_e(&lRBGrid);   lres_e.Checkerboard()=Even;
+
+    pickCheckerboard(Odd,vsrc_o,src);
+    pickCheckerboard(Odd,lsrc_o,srclex);
+
+    Dw.Meooe   (vsrc_o,vres_e);
+    Dwlex.Meooe(lsrc_o,lres_e);
+
+    RealD vn = norm2(vres_e);
+    RealD ln = norm2(lres_e);
+    RealD n  = fabs(vn-ln)/vn;
+    std::cout << GridLogMessage << "Meooe |simd|^2 " << vn << "  |lex|^2 " << ln
+	      << "  relative " << n << std::endl;
+    GRID_ASSERT( n < tol );
+  }
+
+  std::cout << GridLogMessage << "Test_wilson_lex: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}