4D pseudo fermion

Worth an effort to understand if the different stochastic estimator is cleaner.

CLAUDE.md

@@ -31,12 +31,18 @@ Key configure options:
 | `--with-gmp=`, `--with-mpfr=`, `--with-fftw=`, `--with-lime=` | paths to libs |
 | `--enable-hdf5`, `--enable-mkl`, `--enable-lapack` | optional features |
 
+GPU builds additionally need `--enable-gen-simd-width=64` (sets 512-bit SIMD width for GPU warp/wavefront sizing) and `--enable-unified=no --enable-shm=nvlink` for multi-GPU runs.
+
+To speed up compilation, `--disable-fermion-reps --disable-gparity` skips instantiating G-parity and higher-representation fermion operators.
+
 Platform recipes from `README.md`:
 - **KNL**: `--enable-simd=KNL --enable-comms=mpi3-auto --enable-mkl`
 - **Skylake/Haswell**: `--enable-simd=AVX512` or `AVX2` + `--enable-comms=mpi3-auto`
 - **AMD EPYC**: `--enable-simd=AVX2 --enable-comms=mpi3`
 - **A64FX (Fugaku)**: `--enable-simd=A64FX --enable-comms=mpi3 --enable-shm=shmget` (see `SVE_README.txt`)
 
+Complete, working `configure` invocations for specific HPC systems (Frontier/ROCm, Perlmutter/CUDA, Summit, SDCC-A100, etc.) live in `systems/<platform>/config-command`. These are the canonical references for production builds.
+
 Required external libs: GMP, MPFR, OpenSSL, zlib.
 
 ### Use `systems/` for real machines
@@ -59,9 +65,9 @@ Note the GPU builds use `--enable-simd=GPU --enable-gen-simd-width=64`, so `Nsim
 
 ```bash
 # From build directory
-make check                          # runs only Test_simd, Test_cshift, Test_stencil, Test_dwf_mixedcg_prec
-make -C tests/<subdir> tests        # build (not run) the tests in a subdirectory
-make bench                          # build benchmarks
+make check                          # root-level tests (Test_simd, Test_cshift, etc.)
+make -C tests/<subdir> tests        # build tests in a subdirectory
+make tests                          # build all tests across all subdirectories
 ./tests/core/Test_simd              # run a single test binary directly
 mpirun -n 4 ./tests/core/Test_cshift --grid 16.16.16.16 --mpi 1.1.1.4
 ```
@@ -84,7 +90,7 @@ Tests and benchmarks that need optional fermion representations are guarded by `
 
 4. **Cartesian/comms layer** (`Grid/cartesian/`, `Grid/communicator/`) — `GridCartesian` holds the MPI topology and local/global geometry. `Grid/cshift/` implements nearest-neighbour halo exchange; `Grid/stencil/` is the optimised multi-hop stencil used by Dirac operators.
 
-5. **Algorithm layer** (`Grid/algorithms/`) — iterative solvers (CG, GMRES, BiCGSTAB, mixed-precision), eigensolvers (Lanczos, LAPACK), FFT, smearing.
+5. **Algorithm layer** (`Grid/algorithms/`) — iterative solvers (CG, GMRES, BiCGSTAB, mixed-precision), eigensolvers (Lanczos, LAPACK), FFT, smearing, and multigrid.
 
 6. **QCD layer** (`Grid/qcd/`) — gauge and fermion actions, HMC integrators, observables.
 
@@ -98,6 +104,11 @@ Tests and benchmarks that need optional fermion representations are guarded by `
 - `observables/` — Polyakov loop, plaquette, topological charge
 
 ### GPU acceleration and the view/memory-manager discipline
+### Multigrid (`Grid/algorithms/multigrid/`)
+
+Aggregation-based algebraic multigrid for Wilson-type fermions. Key files: `CoarsenedMatrix.h` (coarse operator), `GeneralCoarsenedMatrix.h` and `GeneralCoarsenedMatrixMultiRHS.h` (general coarsening supporting multi-RHS solves), `Aggregates.h` (near-null vector construction), `Geometry.h` (coarse-grid geometry). `MultiGrid.h` is the top-level include.
+
+### GPU acceleration
 
 GPU support is injected via macros in `Grid/threads/Accelerator.h` — `accelerator_for(i, n, nsimd, {...})`, `accelerator_forNB` (non-blocking, must be followed by `accelerator_barrier()`), `accelerator_for2dNB`, and `accelerator_inline`. On a CPU build these degrade to `thread_for` (OpenMP). Unified virtual memory is on by default (`--enable-unified=yes`); device-aware MPI (`--enable-accelerator-aware-mpi`) avoids device→host copies on transfers.
 
@@ -127,6 +138,20 @@ Modes are `AcceleratorRead/Write/WriteDiscard` and `CpuRead/Write/WriteDiscard`.
 | `correctness-verification.md` | Double-run fingerprinting, per-packet checksums, flight recorder |
 | `hang-diagnosis.md` | Diagnosing MPI/accelerator hangs |
 
+The key loop macros (defined in `Grid/threads/Accelerator.h`) are:
+- `accelerator_for(iter, num, nsimd, {...})` — maps to CUDA/HIP kernel or OpenMP loop; `nsimd` is the innermost SIMD lane count
+- `accelerator_forNB(...)` — non-blocking variant (no implicit barrier)
+- `accelerator_for2dNB(iter1, num1, iter2, num2, nsimd, {...})` — 2D kernel launch
+- `thread_for(iter, num, {...})` — CPU OpenMP loop (never dispatches to GPU)
+
+On CPU builds, `accelerator_for` aliases to `thread_for`.
+
+### Solver patterns
+
+`SchurRedBlack` (`Grid/algorithms/iterative/SchurRedBlack.h`) implements red-black (even/odd) preconditioning for fermion operators. Most production fermion solves use `SchurRedBlackDiagMooeeSolve` or similar wrappers that internally call a `ConjugateGradient` on the Schur complement.
+
+Mixed-precision solvers (`ConjugateGradientMixedPrec`, `BiCGSTABMixedPrec`) drive a double-precision outer loop with single-precision inner solves.
+
 ### Memory and I/O
 
 - `Grid/allocator/` — aligned/NUMA-aware allocators; caching allocator via `--enable-alloc-cache`
@@ -146,7 +171,9 @@ Every programme is wrapped in `Grid_init(&argc, &argv)` / `Grid_finalize()` (`Gr
 ## Key Conventions
 
 - **C++17** is required throughout.
-- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`; the explicit instantiation `.cc` files live under `Grid/qcd/action/fermion/instantiation/<Impl>/` and are selected by `scripts/filelist`.
+- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`.
+- **Tensor indices are positional, not labelled.** The `Grid/tensors/` arithmetic recurses structurally over the `iScalar`/`iVector`/`iMatrix` nest: each level defines only the {scalar,vector,matrix}² products at its own level, with element types resolved by automatic type deduction, so every colour/spin/lorentz combination composes from ~200 lines (versus the pre-C++11 QDP++/PETE approach of machine-generating every case). An index's meaning derives entirely from its nesting depth counted from the outside; `iScalar` is the identity/broadcast case at every level. Never insert or remove a nesting level casually — the multiplication tables contract by position.
+- **Multigrid coarsening deepens the tensor nest by one level.** A coarse site vector is `iVector<CComplex,nbasis>`, and `innerProduct` on it returns `iScalar<CComplex>` — one level deeper than the fine block scalar. So the block-inner-product scalar type gains one `iScalar` wrapper per MG level (fine: `vTComplex`; level 2: `iScalar<vTComplex>`; see `examples/Example_pvdagm_3level.cc`). When calling `blockInnerProduct`/`blockZAXPY`/`blockOrthogonalise` on coarse fields, the coarse scalar type must match `decltype(innerProduct(siteVector(),siteVector()))` exactly; a wrong depth fails to compile (no viable `operator=` deep in the instantiation chain) rather than mis-contracting.
 - The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`.
 - Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI.
 - Logging is stream-based, not macro-based: `std::cout << GridLogMessage << ... << std::endl;`. Channels declared in `Grid/log/Log.h` include `GridLogError`, `GridLogWarning`, `GridLogDebug`, `GridLogPerformance`, `GridLogIterative`, `GridLogSolver`, `GridLogHMC`, `GridLogComms`, `GridLogMemory`, `GridLogDslash`, `GridLogIRL`, `GridLogMG`. A subset is switched on at runtime with e.g. `--log Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours` (names given without the `GridLog` prefix).

Grid/algorithms/blas/BatchedBlas.h

@@ -695,6 +695,228 @@ public:
      RealD bytes = 1.0*sizeof(ComplexF)*(m*k+k*n+m*n)*batchCount;
   }
   
+  ///////////////////////////////////////////////////////////////////////////////////
+  // Explicit-leading-dimension complex single GEMM.
+  //
+  // A,B,C may be SLICES of larger parent allocations: lda/ldb/ldc are the
+  // PARENT strides (>= the compact values the ten-argument overload derives).
+  // Motivating use: software split-K for tiny-output/huge-K dense multiplies
+  // (arXiv:2409.03904 fig 11; cf MultiRHSBlockCGLinalg) -- batch over K-chunks
+  // of a dense slab by pointer offset j*Kchunk with lda = the full K extent,
+  // then reduce the partial C's.  Backends pass lda straight through; only
+  // the compact overload invented them.
+  ///////////////////////////////////////////////////////////////////////////////////
+  void gemmBatched(GridBLASOperation_t OpA,
+		   GridBLASOperation_t OpB,
+		   int m,int n, int k,
+		   ComplexF alpha,
+		   deviceVector<ComplexF*> &Amk, int lda,
+		   deviceVector<ComplexF*> &Bkn, int ldb,
+		   ComplexF beta,
+		   deviceVector<ComplexF*> &Cmn, int ldc,
+		   GridBLASPrecision_t precision = GridBLAS_PRECISION_DEFAULT)
+  {
+    RealD t2=usecond();
+    int32_t batchCount = Amk.size();
+
+    GRID_ASSERT( lda >= ((OpA==GridBLAS_OP_N) ? m : k) );
+    GRID_ASSERT( ldb >= ((OpB==GridBLAS_OP_N) ? k : n) );
+    GRID_ASSERT( ldc >= m );
+
+    static deviceVector<ComplexF> alpha_p(1);
+    static deviceVector<ComplexF> beta_p(1);
+    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(ComplexF));
+    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(ComplexF));
+    RealD t0=usecond();
+
+    GRID_ASSERT(Bkn.size()==batchCount);
+    GRID_ASSERT(Cmn.size()==batchCount);
+#ifdef GRID_HIP
+    GRID_ASSERT(precision == GridBLAS_PRECISION_DEFAULT);
+    hipblasOperation_t hOpA;
+    hipblasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
+#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
+    auto err = hipblasCgemmBatched(gridblasHandle,
+				   hOpA,
+				   hOpB,
+				   m,n,k,
+				   (hipComplex *) &alpha_p[0],
+				   (hipComplex **)&Amk[0], lda,
+				   (hipComplex **)&Bkn[0], ldb,
+				   (hipComplex *) &beta_p[0],
+				   (hipComplex **)&Cmn[0], ldc,
+				   batchCount);
+#else
+    auto err = hipblasCgemmBatched(gridblasHandle,
+                                   hOpA,
+                                   hOpB,
+                                   m,n,k,
+                                   (hipblasComplex *) &alpha_p[0],
+                                   (hipblasComplex **)&Amk[0], lda,
+                                   (hipblasComplex **)&Bkn[0], ldb,
+                                   (hipblasComplex *) &beta_p[0],
+                                   (hipblasComplex **)&Cmn[0], ldc,
+                                   batchCount);
+#endif
+    GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_CUDA
+    cublasOperation_t hOpA;
+    cublasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
+    cublasStatus_t err;
+    if (precision == GridBLAS_PRECISION_DEFAULT) {
+      err = cublasCgemmBatched(gridblasHandle,
+			       hOpA,
+			       hOpB,
+			       m,n,k,
+			       (cuComplex *) &alpha_p[0],
+			       (cuComplex **)&Amk[0], lda,
+			       (cuComplex **)&Bkn[0], ldb,
+			       (cuComplex *) &beta_p[0],
+			       (cuComplex **)&Cmn[0], ldc,
+			       batchCount);
+    } else {
+      cublasComputeType_t compute_precision = toDataType(precision);
+      err = cublasGemmBatchedEx(gridblasHandle,
+				hOpA,
+				hOpB,
+				m,n,k,
+				(void *) &alpha_p[0],
+				(void **)&Amk[0], CUDA_C_32F, lda,
+				(void **)&Bkn[0], CUDA_C_32F, ldb,
+				(void *) &beta_p[0],
+				(void **)&Cmn[0], CUDA_C_32F, ldc,
+				batchCount, compute_precision, CUBLAS_GEMM_DEFAULT);
+    }
+    GRID_ASSERT(err==CUBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_SYCL
+    GRID_ASSERT(precision == GridBLAS_PRECISION_DEFAULT);
+    int64_t m64=m;
+    int64_t n64=n;
+    int64_t k64=k;
+    int64_t lda64=lda;
+    int64_t ldb64=ldb;
+    int64_t ldc64=ldc;
+    int64_t batchCount64=batchCount;
+
+    oneapi::mkl::transpose iOpA;
+    oneapi::mkl::transpose iOpB;
+
+    if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
+    if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
+    if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
+    if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
+    if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
+    if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
+
+    oneapi::mkl::blas::column_major::gemm_batch(*gridblasHandle,
+						&iOpA,
+						&iOpB,
+						&m64,&n64,&k64,
+						(ComplexF *) &alpha_p[0],
+						(const ComplexF **)&Amk[0], (const int64_t *)&lda64,
+						(const ComplexF **)&Bkn[0], (const int64_t *)&ldb64,
+						(ComplexF *) &beta_p[0],
+						(ComplexF **)&Cmn[0], (const int64_t *)&ldc64,
+						(int64_t)1,&batchCount64,std::vector<sycl::event>());
+    synchronise();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    GRID_ASSERT(precision == GridBLAS_PRECISION_DEFAULT);
+    // Reference implementation: Eigen with explicit outer stride
+    typedef Eigen::Map<Eigen::MatrixXcf,0,Eigen::OuterStride<> > eMat;
+    if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+	else
+	  eCmn = alpha * eAmk * eBkn ;
+	});
+    } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn ;
+	else
+	  eCmn = alpha * eAmk.adjoint() * eBkn ;
+	});
+    } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],k,n,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn ;
+	else
+	  eCmn = alpha * eAmk.transpose() * eBkn ;
+	});
+    } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint() ;
+	else
+	  eCmn = alpha * eAmk * eBkn.adjoint() ;
+	});
+    } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_T) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],m,k,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk * eBkn.transpose() ;
+	else
+	  eCmn = alpha * eAmk * eBkn.transpose() ;
+	});
+    } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint() ;
+	else
+	  eCmn = alpha * eAmk.adjoint() * eBkn.adjoint() ;
+	} );
+    } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_T) ) {
+      thread_for (p, batchCount, {
+	eMat eAmk(Amk[p],k,m,Eigen::OuterStride<>(lda));
+	eMat eBkn(Bkn[p],n,k,Eigen::OuterStride<>(ldb));
+	eMat eCmn(Cmn[p],m,n,Eigen::OuterStride<>(ldc));
+	if (std::abs(beta) != 0.0)
+	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn.transpose() ;
+	else
+	  eCmn = alpha * eAmk.transpose() * eBkn.transpose() ;
+	} );
+    } else {
+      assert(0);
+    }
+#endif
+    RealD t1=usecond();
+    RealD flops = 8.0*m*n*k*batchCount;
+    RealD bytes = 1.0*sizeof(ComplexF)*(m*k+k*n+m*n)*batchCount;
+  }
+
   ///////////////////////////////////////////////////////////////////////////
   // Single precision real GEMM
   ///////////////////////////////////////////////////////////////////////////

Grid/algorithms/blas/BatchedInverse.h

@@ -0,0 +1,282 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: BatchedInverse.h
+
+    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 */
+#pragma once
+
+#include <Grid/algorithms/blas/BatchedBlas.h>
+
+#ifdef GRID_HIP
+#include <rocsolver/rocsolver.h>
+#endif
+// GRID_CUDA: batched LU inversion lives in cuBLAS (getrfBatched/getriBatched);
+//            cublas_v2.h already included via BatchedBlas.h.
+// GRID_SYCL: oneapi/mkl.hpp already included via BatchedBlas.h (lapack::getrf/getri).
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// GridBLASInverse: cross-platform batched dense matrix inversion.
+//
+// HIGH LEVEL contract (deliberately NOT a getrf/getrs interface): invert a
+// batch of dense N x N matrices IN PLACE,
+//
+//     A[i] <- A[i]^{-1}       i = 0 .. batchCount-1
+//
+// Layout: column major, lda = N, contiguous per batch element; pointer list
+// exactly as GridBLAS::gemmBatched (deviceVector<T*> of device pointers).
+// Each backend chooses HOW:
+//   HIP  : rocSOLVER getrf_batched + getri_batched
+//   CUDA : cuBLAS getrfBatched + getriBatched (out-of-place getri; workspace
+//          hidden here, result copied back so the surface stays in-place)
+//   SYCL : oneMKL LAPACK getrf + getri per batch element (USM, in-order queue)
+//   CPU  : Eigen PartialPivLU (the correctness oracle for all of the above)
+//
+// The int32 vendor-batched entry points bound N < 2^31 (asserted); the huge
+// single-matrix ILP64 path (getrf_64 + blocked identity-getrs harvest, proven
+// in the dense coarse-coarse setup at N=69120) migrates here as a batch==1
+// large-N dispatch in a follow-up -- the recursive Schur leaves are the
+// batched consumers this surface is shaped for.
+//
+// NB GPU-backend call signatures are written to vendor documentation but the
+// air-gapped development loop compiles only the CPU/Eigen path; verify the
+// rocSOLVER/cuBLAS/oneMKL calls against headers on first device compile.
+// Semantics are locked by the CPU unit test (Test_batched_blas).
+///////////////////////////////////////////////////////////////////////////////
+class GridBLASInverse {
+public:
+
+#ifdef GRID_HIP
+  // rocSOLVER runs on a rocblas_handle (distinct type from hipblasHandle_t)
+  static rocblas_handle & Handle(void) {
+    static rocblas_handle h;
+    static int init = 0;
+    if ( !init ) {
+      auto st = rocblas_create_handle(&h);
+      GRID_ASSERT(st == rocblas_status_success);
+      init = 1;
+    }
+    return h;
+  }
+#endif
+#ifdef GRID_CUDA
+  // cuBLAS batched LU shares the GridBLAS handle
+  static cublasHandle_t & Handle(void) {
+    GridBLAS::Init();
+    return GridBLAS::gridblasHandle;
+  }
+#endif
+#ifdef GRID_SYCL
+  static sycl::queue * & Handle(void) {
+    GridBLAS::Init();
+    return GridBLAS::gridblasHandle;
+  }
+#endif
+
+  GridBLASInverse() {};
+  ~GridBLASInverse() {};
+
+  void inverseBatched(int64_t N, deviceVector<ComplexF*> &Amat)
+  {
+    int32_t batchCount = Amat.size();
+    GRID_ASSERT(batchCount > 0);
+
+#ifdef GRID_HIP
+    GRID_ASSERT( N < 2147483647L );
+    rocblas_int n   = (rocblas_int)N;
+    rocblas_int lda = (rocblas_int)N;
+
+    deviceVector<rocblas_int> ipiv((uint64_t)batchCount*N);
+    deviceVector<rocblas_int> info(batchCount);
+
+    auto st1 = rocsolver_cgetrf_batched(Handle(), n, n,
+					(rocblas_float_complex *const *)&Amat[0], lda,
+					&ipiv[0], (rocblas_stride)N,
+					&info[0], batchCount);
+    GRID_ASSERT(st1 == rocblas_status_success);
+    auto st2 = rocsolver_cgetri_batched(Handle(), n,
+					(rocblas_float_complex *const *)&Amat[0], lda,
+					&ipiv[0], (rocblas_stride)N,
+					&info[0], batchCount);
+    GRID_ASSERT(st2 == rocblas_status_success);
+    accelerator_barrier();
+    std::vector<rocblas_int> info_h(batchCount);
+    acceleratorCopyFromDevice(&info[0],&info_h[0],batchCount*sizeof(rocblas_int));
+    for(int i=0;i<batchCount;i++) GRID_ASSERT(info_h[i]==0); // singular pivot => abort loudly
+#endif
+#ifdef GRID_CUDA
+    GRID_ASSERT( N < 2147483647L );
+    int n = (int)N;
+
+    deviceVector<int> ipiv((uint64_t)batchCount*N);
+    deviceVector<int> info(batchCount);
+
+    auto st1 = cublasCgetrfBatched(Handle(), n,
+				   (cuComplex **)&Amat[0], n,
+				   &ipiv[0], &info[0], batchCount);
+    GRID_ASSERT(st1 == CUBLAS_STATUS_SUCCESS);
+
+    // getri is OUT of place: hidden workspace keeps the surface in-place
+    deviceVector<ComplexF>  work((uint64_t)batchCount*N*N);
+    deviceVector<ComplexF*> Cptr(batchCount);
+    std::vector<ComplexF*>  Cptr_h(batchCount);
+    std::vector<ComplexF*>  Aptr_h(batchCount);
+    for(int i=0;i<batchCount;i++) Cptr_h[i] = &work[(uint64_t)i*N*N];
+    acceleratorCopyToDevice(&Cptr_h[0],&Cptr[0],batchCount*sizeof(ComplexF*));
+    acceleratorCopyFromDevice(&Amat[0],&Aptr_h[0],batchCount*sizeof(ComplexF*));
+
+    auto st2 = cublasCgetriBatched(Handle(), n,
+				   (const cuComplex *const *)&Amat[0], n,
+				   &ipiv[0],
+				   (cuComplex **)&Cptr[0], n,
+				   &info[0], batchCount);
+    GRID_ASSERT(st2 == CUBLAS_STATUS_SUCCESS);
+    accelerator_barrier();
+    std::vector<int> info_h(batchCount);
+    acceleratorCopyFromDevice(&info[0],&info_h[0],batchCount*sizeof(int));
+    for(int i=0;i<batchCount;i++) GRID_ASSERT(info_h[i]==0);
+    for(int i=0;i<batchCount;i++)
+      acceleratorCopyDeviceToDevice(Cptr_h[i],Aptr_h[i],(uint64_t)N*N*sizeof(ComplexF));
+#endif
+#ifdef GRID_SYCL
+    // Per-element oneMKL LAPACK on the in-order queue; group API optimisation later.
+    sycl::queue *q = Handle();
+    std::vector<ComplexF*> Aptr_h(batchCount);
+    acceleratorCopyFromDevice(&Amat[0],&Aptr_h[0],batchCount*sizeof(ComplexF*));
+
+    int64_t lwf = oneapi::mkl::lapack::getrf_scratchpad_size<std::complex<float> >(*q,N,N,N);
+    int64_t lwi = oneapi::mkl::lapack::getri_scratchpad_size<std::complex<float> >(*q,N,N);
+    deviceVector<ComplexF> scratchf(lwf);
+    deviceVector<ComplexF> scratchi(lwi);
+    deviceVector<int64_t>  ipiv(N);
+    for(int i=0;i<batchCount;i++){
+      oneapi::mkl::lapack::getrf(*q,N,N,(std::complex<float>*)Aptr_h[i],N,&ipiv[0],
+				 (std::complex<float>*)&scratchf[0],lwf);
+      oneapi::mkl::lapack::getri(*q,N,  (std::complex<float>*)Aptr_h[i],N,&ipiv[0],
+				 (std::complex<float>*)&scratchi[0],lwi);
+    }
+    q->wait();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    // Reference implementation; the oracle the unit test locks semantics with.
+    thread_for (p, batchCount, {
+      Eigen::Map<Eigen::MatrixXcf> eA(Amat[p],N,N);
+      Eigen::PartialPivLU<Eigen::MatrixXcf> lu(eA);
+      eA = lu.inverse();
+      });
+#endif
+  }
+
+  void inverseBatched(int64_t N, deviceVector<ComplexD*> &Amat)
+  {
+    int32_t batchCount = Amat.size();
+    GRID_ASSERT(batchCount > 0);
+
+#ifdef GRID_HIP
+    GRID_ASSERT( N < 2147483647L );
+    rocblas_int n   = (rocblas_int)N;
+    rocblas_int lda = (rocblas_int)N;
+
+    deviceVector<rocblas_int> ipiv((uint64_t)batchCount*N);
+    deviceVector<rocblas_int> info(batchCount);
+
+    auto st1 = rocsolver_zgetrf_batched(Handle(), n, n,
+					(rocblas_double_complex *const *)&Amat[0], lda,
+					&ipiv[0], (rocblas_stride)N,
+					&info[0], batchCount);
+    GRID_ASSERT(st1 == rocblas_status_success);
+    auto st2 = rocsolver_zgetri_batched(Handle(), n,
+					(rocblas_double_complex *const *)&Amat[0], lda,
+					&ipiv[0], (rocblas_stride)N,
+					&info[0], batchCount);
+    GRID_ASSERT(st2 == rocblas_status_success);
+    accelerator_barrier();
+    std::vector<rocblas_int> info_h(batchCount);
+    acceleratorCopyFromDevice(&info[0],&info_h[0],batchCount*sizeof(rocblas_int));
+    for(int i=0;i<batchCount;i++) GRID_ASSERT(info_h[i]==0);
+#endif
+#ifdef GRID_CUDA
+    GRID_ASSERT( N < 2147483647L );
+    int n = (int)N;
+
+    deviceVector<int> ipiv((uint64_t)batchCount*N);
+    deviceVector<int> info(batchCount);
+
+    auto st1 = cublasZgetrfBatched(Handle(), n,
+				   (cuDoubleComplex **)&Amat[0], n,
+				   &ipiv[0], &info[0], batchCount);
+    GRID_ASSERT(st1 == CUBLAS_STATUS_SUCCESS);
+
+    deviceVector<ComplexD>  work((uint64_t)batchCount*N*N);
+    deviceVector<ComplexD*> Cptr(batchCount);
+    std::vector<ComplexD*>  Cptr_h(batchCount);
+    std::vector<ComplexD*>  Aptr_h(batchCount);
+    for(int i=0;i<batchCount;i++) Cptr_h[i] = &work[(uint64_t)i*N*N];
+    acceleratorCopyToDevice(&Cptr_h[0],&Cptr[0],batchCount*sizeof(ComplexD*));
+    acceleratorCopyFromDevice(&Amat[0],&Aptr_h[0],batchCount*sizeof(ComplexD*));
+
+    auto st2 = cublasZgetriBatched(Handle(), n,
+				   (const cuDoubleComplex *const *)&Amat[0], n,
+				   &ipiv[0],
+				   (cuDoubleComplex **)&Cptr[0], n,
+				   &info[0], batchCount);
+    GRID_ASSERT(st2 == CUBLAS_STATUS_SUCCESS);
+    accelerator_barrier();
+    std::vector<int> info_h(batchCount);
+    acceleratorCopyFromDevice(&info[0],&info_h[0],batchCount*sizeof(int));
+    for(int i=0;i<batchCount;i++) GRID_ASSERT(info_h[i]==0);
+    for(int i=0;i<batchCount;i++)
+      acceleratorCopyDeviceToDevice(Cptr_h[i],Aptr_h[i],(uint64_t)N*N*sizeof(ComplexD));
+#endif
+#ifdef GRID_SYCL
+    sycl::queue *q = Handle();
+    std::vector<ComplexD*> Aptr_h(batchCount);
+    acceleratorCopyFromDevice(&Amat[0],&Aptr_h[0],batchCount*sizeof(ComplexD*));
+
+    int64_t lwf = oneapi::mkl::lapack::getrf_scratchpad_size<std::complex<double> >(*q,N,N,N);
+    int64_t lwi = oneapi::mkl::lapack::getri_scratchpad_size<std::complex<double> >(*q,N,N);
+    deviceVector<ComplexD> scratchf(lwf);
+    deviceVector<ComplexD> scratchi(lwi);
+    deviceVector<int64_t>  ipiv(N);
+    for(int i=0;i<batchCount;i++){
+      oneapi::mkl::lapack::getrf(*q,N,N,(std::complex<double>*)Aptr_h[i],N,&ipiv[0],
+				 (std::complex<double>*)&scratchf[0],lwf);
+      oneapi::mkl::lapack::getri(*q,N,  (std::complex<double>*)Aptr_h[i],N,&ipiv[0],
+				 (std::complex<double>*)&scratchi[0],lwi);
+    }
+    q->wait();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    thread_for (p, batchCount, {
+      Eigen::Map<Eigen::MatrixXcd> eA(Amat[p],N,N);
+      Eigen::PartialPivLU<Eigen::MatrixXcd> lu(eA);
+      eA = lu.inverse();
+      });
+#endif
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/multigrid/DenseCoarseMatrix.h

@@ -0,0 +1,702 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/algorithms/multigrid/DenseCoarseMatrix.h
+
+    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 */
+#pragma once
+
+#include <Grid/algorithms/blas/BatchedBlas.h>
+#include <Grid/algorithms/blas/BatchedInverse.h>
+
+#include <unordered_map>
+
+NAMESPACE_BEGIN(Grid);
+
+//////////////////////////////////////////////////////////////////////////////////////
+// DenseCoarseMatrix: a coarsened operator treated as a DENSE matrix -- explicit,
+// row-distributed A^{-1} of a GeneralCoarsenedMatrix.  Library-grade successor of
+// the example-local DistributedDenseInverse (Example_pvdagm_mrhs_3level_dense.cc,
+// FROZEN as the regression baseline).
+//
+// What is new versus the example class:
+//  - Stencil -> dense DIRECT IMPORT.  The coarse operator IS the dense matrix
+//    unrolled: Dense[(s,a),(s+shift_p,b)] += A[p][s]_{a,b}.  Rows of my sites are
+//    assembled from purely LOCAL _A[p] data: no operator applies, no comms -- the
+//    O(N) probe assembly (93 s at N=69120) is retired.  ACCUMULATE (+=) because on
+//    short axes distinct shifts wrap to the same neighbour.  An IMPORT CERTIFICATE
+//    compares the dense apply against Op.M on a NON-CONSTANT vector (a constant one
+//    cannot see a shift-sign error); DENSE_IMPORT_SIGN=-1 flips the convention
+//    without recompiling.
+//  - Split-K apply through GridBLAS.gemmBatched with EXPLICIT leading dimensions
+//    (arXiv:2409.03904 fig 11): the tiny-output/huge-K GEMM Y = slab^T X becomes
+//    DENSE_SPLITK chunk-GEMMs by pointer offset into the resident slab (lda = N),
+//    partials reduced in one accelerator_for.  Platform-agnostic: deviceVector +
+//    GridBLAS run the SAME code on HIP/CUDA/SYCL and CPU(Eigen).
+//  - deviceVector everywhere in the apply path; the ONE surviving naked-HIP block
+//    is the boss inversion buffer (quarantined below, documented).
+//
+// Setup: SLAB_FILE=<stem> loads per-rank <stem>.<rank> (header-guarded N/nrows/
+// nbasis -- the interchange format shared with the frozen example; the STEM must
+// encode cfg/mass/blocking/nbasis, only the header is guarded).  Absent: direct
+// import -> import certificate -> chunked zero-fill+GlobalSum gather streamed to
+// the boss GCD -> cgetrf_64 (ILP64) -> rows of A^{-1} via blocked identity
+// cgetrs_64 + broadcast, each rank keeping the rows of its own sites -> save.
+// VERIFY ||A Ainv x - x||/||x|| runs in BOTH paths (and now certifies the DEVICE
+// slab + split-K path, since the single-RHS apply routes through the same core).
+//
+// Env: SLAB_FILE  DENSE_SPLITK (default 32, snapped to a divisor of N)
+//      DENSE_DEVICE_SUM  DENSE_IMPORT_SIGN  DENSE_APPLY_PROFILE  DENSE_CC_CHECK
+//
+// Eventual internal upgrade (unchanged surface): RecursiveSchur distributed
+// factorisation replacing the single-GCD gather/invert, lifting BOTH the fp32
+// N ~ 90k boss-HBM ceiling AND the CC-grid 256-rank SIMD cap; leaves land on
+// GridBLASInverse::inverseBatched.
+//
+// Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD
+// (iScalar wrappers add no data), so any MG level's coarse operator imports.
+//////////////////////////////////////////////////////////////////////////////////////
+template<class Fobj,class CComplex,int nbasis>
+class DenseCoarseMatrix : public LinearFunction<typename GeneralCoarsenedMatrix<Fobj,CComplex,nbasis>::CoarseVector> {
+public:
+  typedef GeneralCoarsenedMatrix<Fobj,CComplex,nbasis> GeneralCoarseOp;
+  typedef typename GeneralCoarseOp::CoarseVector       Field;
+  typedef typename GeneralCoarseOp::CoarseMatrix       CoarseMatrix;
+  using LinearFunction<Field>::operator();
+  typedef typename Field::vector_object vobj;
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename CoarseMatrix::vector_object Mvobj;
+  typedef typename Mvobj::scalar_object        Msobj;
+
+  GeneralCoarseOp &_Op;             // the coarse operator: stencil source + certificate oracle
+  GridBase *grid;
+  int      nd;
+  int64_t  N;                       // dense rank = gSites * nbasis
+  int      lsites;                  // my local sites
+  int64_t  nrows;                   // my rows = lsites * nbasis
+  std::vector<Coordinate> myLcoor;  // local coordinate of my site ss
+  std::vector<int64_t>    myGsite;  // global lex site index of my site ss
+  std::vector<ComplexF>   slab;     // nrows x N row-major: A during setup, rows of A^{-1} after
+
+  static const int64_t CHUNKROWS = 1024;  // getrs harvest block (trsm efficiency + fewer broadcasts)
+  static const int MRHS_MAX = 32;
+
+  // Apply machinery: resident slab + persistent buffers + AOT split-K pointers.
+  GridBLAS BLAS;
+  deviceVector<ComplexF>  dSlab;
+  deviceVector<ComplexF>  dX;       // N x MRHS_MAX
+  deviceVector<ComplexF>  dY;       // nrows x MRHS_MAX
+  deviceVector<ComplexF>  dPartial; // NK x (nrows x MRHS_MAX)
+  deviceVector<ComplexF*> aptrs;    // slab K-chunk pointers   (lda = N)
+  deviceVector<ComplexF*> xptrs;    // X    K-chunk pointers   (ldb = N)
+  deviceVector<ComplexF*> cptrs;    // partial buffers         (ldc = nrows)
+  std::vector<ComplexF>   hX;
+  std::vector<ComplexF>   hY;
+  int NK;                           // split-K chunk count (divides N)
+  int devSum;
+
+  DenseCoarseMatrix(GeneralCoarseOp &Op, GridBase *g)
+    : _Op(Op), grid(g)
+  {
+    GRID_ASSERT( sizeof(sobj)  == nbasis*sizeof(ComplexD) );
+    GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(ComplexD) );
+    GRID_ASSERT( grid == Op.Grid() );
+    nd     = grid->_ndimension;
+    N      = grid->gSites() * nbasis;
+    lsites = grid->lSites();
+    nrows  = (int64_t)lsites * nbasis;
+
+    std::cout << GridLogMessage << "DenseCoarseMatrix: N = " << N
+              << " (" << grid->gSites() << " sites x " << nbasis << ")"
+              << "  rows/rank = " << nrows
+              << "  slab = " << (double)nrows*N*sizeof(ComplexF)/1024./1024. << " MB/rank"
+              << std::endl;
+
+    ////////////////////////////////////////////////////////////////////
+    // Enumerate my sites: local coords and global lexicographic indices
+    ////////////////////////////////////////////////////////////////////
+    Coordinate ldims = grid->LocalDimensions();
+    Coordinate gdims = grid->GlobalDimensions();
+    myLcoor.resize(lsites);
+    myGsite.resize(lsites);
+    for(int ss=0; ss<lsites; ss++){
+      Coordinate lcoor(nd);
+      Lexicographic::CoorFromIndex(lcoor, ss, ldims);
+      Coordinate gcoor(nd);
+      for(int d=0; d<nd; d++) gcoor[d] = grid->_lstart[d] + lcoor[d];
+      int64_t gsite;
+      Lexicographic::IndexFromCoor(gcoor, gsite, gdims);
+      myLcoor[ss] = lcoor;
+      myGsite[ss] = gsite;
+    }
+
+    slab.resize((uint64_t)nrows * N);
+
+    double t0 = usecond();
+    ////////////////////////////////////////////////////////////////////
+    // 0. Slab cache: SLAB_FILE=<stem> -> per-rank raw file <stem>.<rank>.
+    //    SAME format as the frozen example (interchange compatible).
+    ////////////////////////////////////////////////////////////////////
+    bool loaded = false;
+    char *sfile = getenv("SLAB_FILE");
+    std::string slabfile;
+    if (sfile) {
+      slabfile = std::string(sfile) + "." + std::to_string(grid->ThisRank());
+      FILE *f = fopen(slabfile.c_str(),"rb");
+      if (f) {
+        int64_t hdr[4] = {0,0,0,0};
+        GRID_ASSERT( fread(hdr,sizeof(int64_t),4,f) == 4 );
+        GRID_ASSERT( hdr[0] == (int64_t)0x44454E5345 );  // magic "DENSE"
+        GRID_ASSERT( hdr[1] == N && hdr[2] == (int64_t)nrows && hdr[3] == (int64_t)nbasis );
+        uint64_t nelem = (uint64_t)nrows * N;
+        GRID_ASSERT( fread(&slab[0], sizeof(ComplexF), nelem, f) == nelem );
+        fclose(f);
+        loaded = true;
+        std::cout << GridLogMessage << "DenseCoarseMatrix: slab loaded from "
+                  << slabfile << " -- skipping import/factor/solve" << std::endl;
+      } else {
+        std::cout << GridLogMessage << "DenseCoarseMatrix: slab cache " << slabfile
+                  << " absent -- full setup, will write it" << std::endl;
+      }
+    }
+    if (!loaded) {
+      ImportDense();          // slab <- my rows of A   (LOCAL, no comms)
+      ImportCertificate();    // dense apply == Op.M on a non-constant vector
+      InvertDense();          // slab <- my rows of A^{-1}
+      double t1 = usecond();
+      if (sfile) {
+        FILE *f = fopen(slabfile.c_str(),"wb");
+        GRID_ASSERT(f != nullptr);
+        int64_t hdr[4] = { (int64_t)0x44454E5345, N, (int64_t)nrows, (int64_t)nbasis };
+        GRID_ASSERT( fwrite(hdr,sizeof(int64_t),4,f) == 4 );
+        uint64_t nelem = (uint64_t)nrows * N;
+        GRID_ASSERT( fwrite(&slab[0], sizeof(ComplexF), nelem, f) == nelem );
+        fclose(f);
+        std::cout << GridLogMessage << "DenseCoarseMatrix: slab written to " << slabfile << std::endl;
+      }
+      std::cout << GridLogMessage << "DenseCoarseMatrix: import+invert took "
+                << (t1-t0)/1.0e6 << " s" << std::endl;
+    }
+
+    ////////////////////////////////////////////////////////////////////
+    // Device residency + persistent apply buffers + AOT split-K pointers
+    ////////////////////////////////////////////////////////////////////
+    {
+      uint64_t sbytes = (uint64_t)nrows * N * sizeof(ComplexF);
+      dSlab.resize((uint64_t)nrows*N);
+      acceleratorCopyToDevice(&slab[0],&dSlab[0],sbytes);
+
+      // DENSE_SPLITK: requested chunk count, snapped DOWN to a divisor of N.
+      int req = getenv("DENSE_SPLITK") ? atoi(getenv("DENSE_SPLITK")) : 32;
+      if (req < 1) req = 1;
+      NK = 1;
+      for(int j=1;j<=req;j++) if ( (N % j) == 0 ) NK = j;
+      int64_t Kc = N / NK;
+
+      dX.resize((uint64_t)N*MRHS_MAX);
+      dY.resize((uint64_t)nrows*MRHS_MAX);
+      dPartial.resize((uint64_t)NK*nrows*MRHS_MAX);
+      hX.resize((uint64_t)N*MRHS_MAX);
+      hY.resize((uint64_t)nrows*MRHS_MAX);
+
+      aptrs.resize(NK); xptrs.resize(NK); cptrs.resize(NK);
+      std::vector<ComplexF*> h(NK);
+      for(int j=0;j<NK;j++) h[j] = &dSlab[0]    + (uint64_t)j*Kc;              // K-offset, lda=N
+      acceleratorCopyToDevice(&h[0],&aptrs[0],NK*sizeof(ComplexF*));
+      for(int j=0;j<NK;j++) h[j] = &dX[0]       + (uint64_t)j*Kc;              // K-offset, ldb=N
+      acceleratorCopyToDevice(&h[0],&xptrs[0],NK*sizeof(ComplexF*));
+      for(int j=0;j<NK;j++) h[j] = &dPartial[0] + (uint64_t)j*nrows*MRHS_MAX;  // compact, ldc=nrows
+      acceleratorCopyToDevice(&h[0],&cptrs[0],NK*sizeof(ComplexF*));
+
+      devSum = getenv("DENSE_DEVICE_SUM") ? 1 : 0;
+      std::cout << GridLogMessage << "DenseCoarseMatrix: slab resident on device ("
+                << sbytes/1024./1024. << " MB/rank), split-K NK=" << NK << " (Kc=" << Kc << "); "
+                << (devSum ? "DEVICE-buffer allreduce (GPU-aware MPI)" : "host allreduce")
+                << std::endl;
+    }
+
+    ////////////////////////////////////////////////////////////////////
+    // VERIFY: || A (Ainv x) - x || / ||x|| through the DEVICE split-K core.
+    ////////////////////////////////////////////////////////////////////
+    {
+      Field x(grid); Field y(grid); Field z(grid);
+      x = ComplexD(1.0,0.0);
+      double ta = usecond();
+      (*this)(x, y);
+      double tb = usecond();
+      _Op.M(y, z);
+      z = z - x;
+      RealD rel = std::sqrt(norm2(z)/norm2(x));
+      std::cout << GridLogMessage << "DenseCoarseMatrix: VERIFY ||A Ainv x - x||/||x|| = "
+                << rel << "   (one apply took " << (tb-ta)/1000.0 << " ms)" << std::endl;
+      GRID_ASSERT(rel < 1.0e-2);
+    }
+    std::cout << GridLogMessage << "DenseCoarseMatrix: setup complete, total "
+              << (usecond()-t0)/1.0e6 << " s" << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 1. Direct stencil -> dense import of MY ROWS of A (no comms):
+  //      Dense[(s,a),(wrap(s+shift_p),b)] += A[p][s]_{a,b}
+  ////////////////////////////////////////////////////////////////////
+  void ImportDense(void)
+  {
+    double t = -usecond();
+    Coordinate gdims = grid->GlobalDimensions();
+    int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1;
+    GRID_ASSERT( sign==1 || sign==-1 );
+
+    uint64_t nelem = (uint64_t)nrows * N;
+    thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); });
+
+    for(int p=0; p<_Op.geom.npoint; p++){
+      Coordinate shift = _Op.geom.shifts[p];
+      autoView(Av, _Op._A[p], CpuRead);
+      thread_for(ss, lsites, {
+        Coordinate ncoor(nd);
+        for(int d=0; d<nd; d++){
+          int64_t g = grid->_lstart[d] + myLcoor[ss][d] + sign*shift[d];
+          ncoor[d] = (int)((g % gdims[d] + gdims[d]) % gdims[d]);
+        }
+        int64_t nsite;
+        Lexicographic::IndexFromCoor(ncoor, nsite, gdims);
+        Msobj m;
+        peekLocalSite(m, Av, myLcoor[ss]);
+        ComplexD *md = (ComplexD *)&m;
+        for(int a=0; a<nbasis; a++){
+          ComplexF *row = &slab[(uint64_t)(ss*nbasis+a)*N + nsite*nbasis];
+          for(int b=0; b<nbasis; b++)
+            row[b] += ComplexF(md[a*nbasis+b]);   // += : wrapped shifts may collide
+        }
+      });
+    }
+    t += usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: stencil->dense import took "
+              << t/1.0e6 << " s  (" << _Op.geom.npoint << " points, local, no comms)" << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 2. IMPORT CERTIFICATE: dense rows vs Op.M on a NON-CONSTANT vector.
+  //    (Constant x has x[s+d]==x[s-d]: blind to a shift-sign error.)
+  ////////////////////////////////////////////////////////////////////
+  void ImportCertificate(void)
+  {
+    Field x(grid); Field Ax(grid); Field Dx(grid);
+    for(int ss=0; ss<lsites; ss++){
+      sobj s;
+      for(int b=0; b<nbasis; b++){
+        double ph = 0.37*(double)(myGsite[ss]*nbasis+b);
+        ((ComplexD *)&s)[b] = ComplexD(std::cos(ph),std::sin(0.61*ph));
+      }
+      pokeLocalSite(s, x, myLcoor[ss]);
+    }
+    // gather full x (zero-fill + exact GlobalSum), dense rows on host
+    std::vector<ComplexD> xh((uint64_t)N, ComplexD(0.0,0.0));
+    for(int ss=0; ss<lsites; ss++){
+      sobj s;
+      peekLocalSite(s, x, myLcoor[ss]);
+      for(int b=0; b<nbasis; b++) xh[ myGsite[ss]*nbasis + b ] = ((ComplexD *)&s)[b];
+    }
+    grid->GlobalSumVector(&xh[0], (int)N);
+    std::vector<ComplexD> yh(nrows);
+    thread_for(r, nrows, {
+      ComplexD acc(0.0,0.0);
+      const ComplexF *row = &slab[(uint64_t)r * N];
+      for(int64_t j=0; j<N; j++) acc += ComplexD(row[j]) * xh[j];
+      yh[r] = acc;
+    });
+    for(int ss=0; ss<lsites; ss++){
+      sobj s;
+      for(int b=0; b<nbasis; b++) ((ComplexD *)&s)[b] = yh[ss*nbasis+b];
+      pokeLocalSite(s, Dx, myLcoor[ss]);
+    }
+    _Op.M(x, Ax);
+    Field d(grid); d = Dx - Ax;
+    RealD rel = std::sqrt(norm2(d)/norm2(Ax));
+    std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE ||Dense x - A x||/||A x|| = "
+              << rel << std::endl;
+    if ( rel >= 1.0e-3 ) {
+      std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE FAILED. If O(1), the "
+                << "stencil shift-sign convention is opposite: rerun with DENSE_IMPORT_SIGN=-1"
+                << std::endl;
+    }
+    GRID_ASSERT(rel < 1.0e-3);
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 3. Invert: chunked zero-fill+GlobalSum gather of A streamed to the
+  //    boss GCD, cgetrf_64 (ILP64), rows of A^{-1} by blocked identity
+  //    cgetrs_64 + broadcast; each rank keeps its own rows (in `slab`,
+  //    overwriting A).  Ported proven path from the frozen example.
+  ////////////////////////////////////////////////////////////////////
+  void InvertDense(void)
+  {
+    double t1 = usecond();
+    int boss = grid->IsBoss();
+    std::vector<ComplexF> Afull;
+#ifdef GRID_HIP
+    // QUARANTINED naked HIP: the boss-only N^2 inversion buffer (34GB at
+    // N=65536) must come from raw HBM; EvictAll flushes the device-copy
+    // layer to make the window.  (FreePool of the allocator free-list
+    // awaits the type-dispatched fix.)  Confined to setup; the apply path
+    // is pure Grid primitives.
+    rocblas_float_complex *dA = nullptr;
+    rocblas_float_complex *dB = nullptr;
+    int64_t *dIpiv = nullptr;
+    uint64_t Abytes = (uint64_t)N * N * sizeof(ComplexF);
+    MemoryManager::EvictAll();
+    if (boss) {
+      auto aerr = hipMalloc((void **)&dA, Abytes);
+      if (aerr != hipSuccess) {
+        std::cout << GridLogMessage << "DenseCoarseMatrix: hipMalloc of "
+                  << Abytes/1024./1024./1024. << " GB FAILED -- reduce --device-mem" << std::endl;
+        GRID_ASSERT(aerr == hipSuccess);
+      }
+      std::cout << GridLogMessage << "DenseCoarseMatrix: device inversion buffer allocated ("
+                << Abytes/1024./1024./1024. << " GB)" << std::endl;
+    }
+#else
+    if (boss) Afull.resize((uint64_t)N * N);
+#endif
+    {
+      std::unordered_map<int64_t,int> rowmap;   // global row -> my slab row
+      for(int ss=0; ss<lsites; ss++)
+        for(int a=0; a<nbasis; a++)
+          rowmap[ myGsite[ss]*nbasis + a ] = ss*nbasis + a;
+
+      std::vector<ComplexF> chunk((uint64_t)CHUNKROWS * N);
+      for(int64_t row0=0; row0<N; row0+=CHUNKROWS){
+        int64_t nrow = std::min(CHUNKROWS, N-row0);
+        uint64_t nelem = (uint64_t)nrow * N;
+        for(uint64_t i=0;i<nelem;i++) chunk[i]=ComplexF(0.0,0.0);
+        for(int64_t r=row0; r<row0+nrow; r++){
+          auto it = rowmap.find(r);
+          if (it != rowmap.end()) {
+            uint64_t src = (uint64_t)(it->second) * N;
+            uint64_t dst = (uint64_t)(r-row0) * N;
+            for(int64_t j=0;j<N;j++) chunk[dst+j] = slab[src+j];
+          }
+        }
+        grid->GlobalSumVector(&chunk[0], (int)nelem);
+        if (boss) {
+#ifdef GRID_HIP
+          GRID_ASSERT( hipMemcpy((char *)dA + (uint64_t)row0*N*sizeof(ComplexF),
+                                 &chunk[0], nelem*sizeof(ComplexF),
+                                 hipMemcpyHostToDevice) == hipSuccess );
+#else
+          uint64_t dst = (uint64_t)row0 * N;
+          for(uint64_t i=0;i<nelem;i++) Afull[dst+i] = chunk[i];
+#endif
+        }
+      }
+    }
+    double t2 = usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: gather to boss took "
+              << (t2-t1)/1.0e6 << " s" << std::endl;
+
+    ////////////////////////////////////////////////////////////////////
+    // Factor in place on the boss (fp32, ILP64).  Row-major buffer handed
+    // to column-major LAPACK => LU of A^T.
+    ////////////////////////////////////////////////////////////////////
+    if (boss) {
+#ifdef GRID_HIP
+      std::cout << GridLogMessage << "DenseCoarseMatrix: rocSOLVER cgetrf_64 (ILP64 LU) N=" << N
+                << " in place on resident device buffer" << std::endl;
+      rocblas_handle handle = GridBLASInverse::Handle();
+      int64_t *dInfo;
+      GRID_ASSERT( hipMalloc((void **)&dIpiv, N*sizeof(int64_t)) == hipSuccess );
+      GRID_ASSERT( hipMalloc((void **)&dInfo, sizeof(int64_t))   == hipSuccess );
+      auto st1 = rocsolver_cgetrf_64(handle, (int64_t)N, (int64_t)N, dA, (int64_t)N, dIpiv, dInfo);
+      hipDeviceSynchronize();
+      int64_t info_h = -1;
+      hipMemcpy(&info_h, dInfo, sizeof(int64_t), hipMemcpyDeviceToHost);
+      std::cout << GridLogMessage << "DenseCoarseMatrix: cgetrf_64 status " << (int)st1
+                << " info = " << (int)info_h << std::endl;
+      GRID_ASSERT(st1 == rocblas_status_success);
+      GRID_ASSERT(info_h == 0);
+      hipFree(dInfo);
+      GRID_ASSERT( hipMalloc((void **)&dB, (uint64_t)CHUNKROWS*N*sizeof(ComplexF)) == hipSuccess );
+      // dA holds the LU of A^T; rows of A^{-1} are produced blockwise below via
+      // cgetrs_64 on identity-column blocks: A^T X = E => X columns = rows of
+      // A^{-1}, in exactly the linear layout the harvest expects.
+#else
+      // Eigen fallback: small local CPU tests only.
+      std::cout << GridLogMessage << "DenseCoarseMatrix: Eigen fallback inversion N=" << N
+                << (N > 10000 ? "  (WARNING: SLOW; use the HIP/rocSOLVER path)" : "")
+                << std::endl;
+      typedef Eigen::Matrix<std::complex<float>,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> MatF;
+      Eigen::Map<MatF> A(reinterpret_cast<std::complex<float>*>(&Afull[0]), N, N);
+      MatF Ainv = A.inverse();
+      A = Ainv;
+#endif
+    }
+    double t3 = usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: factorisation took "
+              << (t3-t2)/1.0e6 << " s" << std::endl;
+
+    ////////////////////////////////////////////////////////////////////
+    // Blocked solve + broadcast: rows of A^{-1} chunk by chunk; each
+    // rank keeps the rows of its own sites (ownership-aligned).
+    ////////////////////////////////////////////////////////////////////
+    {
+      std::unordered_map<int64_t,int> rowmap;
+      for(int ss=0; ss<lsites; ss++)
+        for(int a=0; a<nbasis; a++)
+          rowmap[ myGsite[ss]*nbasis + a ] = ss*nbasis + a;
+
+      std::vector<ComplexF> chunk((uint64_t)CHUNKROWS * N);
+      for(int64_t row0=0; row0<N; row0+=CHUNKROWS){
+        int64_t nrow = std::min(CHUNKROWS, N-row0);
+        uint64_t nelem = (uint64_t)nrow * N;
+        if (boss) {
+#ifdef GRID_HIP
+          // Identity block E: column j = e_{row0+j}; solve A^T X = E so X's
+          // columns are rows [row0,row0+nrow) of A^{-1}.
+          for(uint64_t i=0;i<nelem;i++) chunk[i] = ComplexF(0.0,0.0);
+          for(int64_t j=0;j<nrow;j++) chunk[(uint64_t)j*N + (uint64_t)(row0+j)] = ComplexF(1.0,0.0);
+          GRID_ASSERT( hipMemcpy(dB, &chunk[0], nelem*sizeof(ComplexF), hipMemcpyHostToDevice) == hipSuccess );
+          auto strs = rocsolver_cgetrs_64(GridBLASInverse::Handle(), rocblas_operation_none,
+                                          (int64_t)N, (int64_t)nrow,
+                                          dA, (int64_t)N, dIpiv, dB, (int64_t)N);
+          GRID_ASSERT(strs == rocblas_status_success);
+          hipDeviceSynchronize();
+          GRID_ASSERT( hipMemcpy(&chunk[0], dB, nelem*sizeof(ComplexF), hipMemcpyDeviceToHost) == hipSuccess );
+#else
+          uint64_t src = (uint64_t)row0 * N;
+          for(uint64_t i=0;i<nelem;i++) chunk[i] = Afull[src+i];
+#endif
+        }
+        grid->Broadcast(0, &chunk[0], nelem*sizeof(ComplexF));
+        for(int64_t r=row0; r<row0+nrow; r++){
+          auto it = rowmap.find(r);
+          if (it != rowmap.end()) {
+            uint64_t dst = (uint64_t)(it->second) * N;
+            uint64_t src = (uint64_t)(r-row0) * N;
+            for(int64_t j=0;j<N;j++) slab[dst+j] = chunk[src+j];
+          }
+        }
+      }
+    }
+#ifdef GRID_HIP
+    if (boss) {
+      if (dA)    hipFree(dA);
+      if (dB)    hipFree(dB);
+      if (dIpiv) hipFree(dIpiv);
+    }
+#endif
+    double t4 = usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: blocked getrs solve+scatter took "
+              << (t4-t3)/1.0e6 << " s" << std::endl;
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // CORE apply on packed data: hX[N x nr] (zero-filled, my sites only)
+  // -> allreduce -> split-K GEMM against the resident slab -> reduce
+  // partials -> hY[nrows x nr] (column major).  Platform-agnostic:
+  // deviceVector + GridBLAS (Eigen fallback on CPU builds).
+  // fp32 allreduce is EXACT: zero-fill assembly gives every element
+  // exactly one contributing rank.
+  ////////////////////////////////////////////////////////////////////
+  void SlabApplyPacked(int nr, double *tprof)
+  {
+    GRID_ASSERT(nr <= MRHS_MAX);
+    uint64_t nX = (uint64_t)N * nr;
+    uint64_t nY = (uint64_t)nrows * nr;
+    int64_t  Kc = N / NK;
+    double t1 = usecond();
+    double t2, t3;
+    if (devSum) {
+      acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      t2 = usecond();
+      grid->GlobalSumVector((ComplexF *)&dX[0], (int)nX);
+      t3 = usecond();
+    } else {
+      grid->GlobalSumVector(&hX[0], (int)nX);
+      t2 = usecond();
+      acceleratorCopyToDevice(&hX[0],&dX[0],nX*sizeof(ComplexF));
+      t3 = usecond();
+    }
+    // Y = op(slab,T) . X : row-major slab (nrows x N) == col-major A^T
+    // (N x nrows, lda=N) => transpose gives the nrows x N operator.
+    // Split-K: NK chunk-GEMMs by pointer offset (AOT lists), then reduce.
+    ComplexF one (1.0,0.0);
+    ComplexF zero(0.0,0.0);
+    BLAS.gemmBatched(GridBLAS_OP_T, GridBLAS_OP_N,
+                     (int)nrows, nr, (int)Kc,
+                     one,  aptrs, (int)N,
+                           xptrs, (int)N,
+                     zero, cptrs, (int)nrows);
+    BLAS.synchronise();
+    {
+      ComplexF *pp = &dPartial[0];
+      ComplexF *py = &dY[0];
+      uint64_t stride = (uint64_t)nrows*MRHS_MAX;
+      int nk = NK;
+      accelerator_for(i, nY, 1, {
+        ComplexF acc(0.0,0.0);
+        for(int j=0;j<nk;j++) acc += pp[(uint64_t)j*stride + i];
+        py[i] = acc;
+      });
+    }
+    double t4 = usecond();
+    acceleratorCopyFromDevice(&dY[0],&hY[0],nY*sizeof(ComplexF));
+    double t5 = usecond();
+    if (tprof) {
+      tprof[0] = devSum ? (t3-t2) : (t2-t1);   // allreduce
+      tprof[1] = devSum ? (t2-t1) : (t3-t2);   // H2D
+      tprof[2] = t4-t3;                        // gemm+reduce
+      tprof[3] = t5-t4;                        // D2H
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Single-RHS apply (also the VERIFY path => certifies device slab).
+  ////////////////////////////////////////////////////////////////////
+  virtual void operator()(const Field &src, Field &psi)
+  {
+    uint64_t nX = (uint64_t)N;
+    thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
+    for(int ss=0; ss<lsites; ss++){
+      sobj s;
+      peekLocalSite(s, src, myLcoor[ss]);
+      for(int b=0; b<nbasis; b++)
+        hX[ myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+    }
+    SlabApplyPacked(1, nullptr);
+    for(int ss=0; ss<lsites; ss++){
+      sobj s;
+      for(int b=0; b<nbasis; b++)
+        ((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]);
+      pokeLocalSite(s, psi, myLcoor[ss]);
+    }
+    if ( getenv("DENSE_CC_CHECK") ) {
+      Field tmp(grid);
+      _Op.M(psi, tmp);
+      tmp = tmp - src;
+      std::cout << GridLogMessage << "DenseCoarseMatrix: apply defect ||A x - b||/||b|| = "
+                << std::sqrt(norm2(tmp)/norm2(src)) << std::endl;
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // Batched (vector-of-fields) apply.
+  ////////////////////////////////////////////////////////////////////
+  void ApplyBatch(std::vector<Field> &src, std::vector<Field> &psi)
+  {
+    int nr = src.size();
+    GRID_ASSERT(nr <= MRHS_MAX);
+    double t0 = usecond();
+    uint64_t nX = (uint64_t)N*nr;
+    thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
+    for(int rr=0; rr<nr; rr++){
+      for(int ss=0; ss<lsites; ss++){
+        sobj s;
+        peekLocalSite(s, src[rr], myLcoor[ss]);
+        for(int b=0; b<nbasis; b++)
+          hX[ (uint64_t)rr*N + myGsite[ss]*nbasis + b ] = ComplexF(((ComplexD *)&s)[b]);
+      }
+    }
+    SlabApplyPacked(nr, nullptr);
+    for(int rr=0; rr<nr; rr++){
+      for(int ss=0; ss<lsites; ss++){
+        sobj s;
+        for(int b=0; b<nbasis; b++)
+          ((ComplexD *)&s)[b] = ComplexD(hY[(uint64_t)rr*nrows + (ss*nbasis+b)]);
+        pokeLocalSite(s, psi[rr], myLcoor[ss]);
+      }
+    }
+    double t1 = usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: batched apply " << nr << " rhs took "
+              << (t1-t0)/1000.0 << " ms  (" << (t1-t0)/1000.0/nr << " ms/rhs)" << std::endl;
+    if ( getenv("DENSE_CC_CHECK") ) {
+      Field tmp(grid);
+      for(int rr=0; rr<nr; rr++){
+        _Op.M(psi[rr], tmp);
+        tmp = tmp - src[rr];
+        std::cout << GridLogMessage << "DenseCoarseMatrix: batch defect["<<rr<<"] = "
+                  << std::sqrt(norm2(tmp)/norm2(src[rr])) << std::endl;
+      }
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // 6D mrhs apply: operates DIRECTLY on the packed 6D field (rhs = dim 0).
+  ////////////////////////////////////////////////////////////////////
+  void ApplyBatch6D(const Field &in6, Field &out6, int nr)
+  {
+    GRID_ASSERT(nr <= MRHS_MAX);
+    GRID_ASSERT(in6.Grid()->_ndimension == nd+1);   // {rhs, s, x,y,z,t}
+    double t0 = usecond();
+    Field &in = const_cast<Field &>(in6);
+    uint64_t nX = (uint64_t)N * nr;
+    thread_for(i, nX, { hX[i]=ComplexF(0.0,0.0); });
+    {
+      autoView(iv, in, CpuRead);
+      Coordinate c6(nd+1);
+      for(int ss=0; ss<lsites; ss++){
+        for(int d=0; d<nd; d++) c6[d+1] = myLcoor[ss][d];
+        for(int rr=0; rr<nr; rr++){
+          c6[0] = rr;
+          sobj s;
+          peekLocalSite(s, iv, c6);
+          for(int b=0; b<nbasis; b++)
+            hX[(uint64_t)rr*N + myGsite[ss]*nbasis + b] = ComplexF(((ComplexD *)&s)[b]);
+        }
+      }
+    }
+    double t1 = usecond();
+    double tprof[4];
+    SlabApplyPacked(nr, tprof);
+    double t5 = usecond();
+    {
+      autoView(ov, out6, CpuWrite);
+      Coordinate c6(nd+1);
+      for(int ss=0; ss<lsites; ss++){
+        for(int d=0; d<nd; d++) c6[d+1] = myLcoor[ss][d];
+        for(int rr=0; rr<nr; rr++){
+          c6[0] = rr;
+          sobj s;
+          for(int b=0; b<nbasis; b++)
+            ((ComplexD *)&s)[b] = ComplexD(hY[(uint64_t)rr*nrows + (ss*nbasis+b)]);  // Y col-major
+          pokeLocalSite(s, ov, c6);
+        }
+      }
+    }
+    double t6 = usecond();
+    std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D " << nr << " rhs took "
+              << (t6-t0)/1000.0 << " ms" << std::endl;
+    if ( getenv("DENSE_APPLY_PROFILE") ) {
+      std::cout << GridLogMessage << "DenseCoarseMatrix: apply6D profile:"
+                << " pack "        << (t1-t0)/1000.0
+                << "  allreduce "  << tprof[0]/1000.0
+                << "  H2D "        << tprof[1]/1000.0
+                << "  gemm+reduce "<< tprof[2]/1000.0
+                << "  D2H "        << tprof[3]/1000.0
+                << "  unpack "     << (t6-t5)/1000.0
+                << "  ms" << std::endl;
+    }
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/qcd/action/pseudofermion/TwoFlavourBosonPseudoFermion.h

@@ -0,0 +1,133 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/pseudofermion/TwoFlavourBosonPseudoFermion.h
+
+    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 */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Two flavour BOSON (wrong-sign) pseudofermion for any FermionOperator B:
+//
+//    S2 = chi^dag Bdag B chi = |B chi|^2
+//
+//    integral  ==>  det( Bdag B )^-1 = |det B|^-2
+//
+// A compensator monomial: supplies an INVERSE determinant with NO solve in
+// the force or the action -- both are matrix multiplies.  The only solve is
+// the heatbath  chi = B^-1 eta,  once per trajectory (for B = the
+// Pauli-Villars operator this is a mass-one solve, trivially cheap).
+//
+// Primary use: two instances with B = PV cancel the |det PV|^2 excess of
+// TwoFlavourPVdagMPseudoFermionAction down to the DWF quotient
+// |det M|^2/|det PV|^2 (two unsquared instances rather than one squared
+// kernel: first powers of PV in the force, milder).  Being generic in B it
+// also serves Hasenbusch-chain compensation at intermediate masses, or any
+// future inverse-det bookkeeping.  (Sibling of the domain-decomposed boson
+// in DomainDecomposedBoundaryTwoFlavourBosonPseudoFermion.h, without the
+// boundary machinery.)
+//
+// Heatbath exact by construction: S2 after refresh = |B B^-1 eta|^2 = |eta|^2.
+///////////////////////////////////////////////////////////////////////////////
+template<class Impl>
+class TwoFlavourBosonPseudoFermionAction : public Action<typename Impl::GaugeField> {
+public:
+  INHERIT_IMPL_TYPES(Impl);
+
+private:
+  FermionOperator<Impl> & BOp;  // the operator whose |det|^-2 is supplied
+
+  LinearFunction<FermionField> &HeatbathSolver; // b -> B^-1 b (heatbath only)
+
+  FermionField Chi; // the pseudo fermion field for this trajectory
+
+public:
+  TwoFlavourBosonPseudoFermionAction(FermionOperator<Impl>  &_BOp,
+				     LinearFunction<FermionField> & HS
+				     ) : BOp(_BOp),
+					 HeatbathSolver(HS),
+					 Chi(_BOp.FermionGrid())
+  {};
+
+  virtual std::string action_name(){return "TwoFlavourBosonPseudoFermionAction";}
+
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
+    return sstream.str();
+  }
+
+  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+    // P(chi) = e^{- chi^dag BdagB chi} ; chi = B^-1 eta ; P(eta) = e^{-eta^dag eta}
+    // e^{-x^2/2 sig^2} => sig^2 = 0.5 ; eta enters with width 1/sqrt(2).
+    RealD scale = std::sqrt(0.5);
+    FermionField eta(BOp.FermionGrid());
+    gaussian(pRNG,eta);
+    eta = eta * scale;
+    refresh(U,eta);
+  }
+
+  // Deterministic-noise variant (test hook):
+  // after this,  S(U) == norm2(eta)  exactly (to solver tolerance).
+  void refresh(const GaugeField &U, const FermionField &eta) {
+    BOp.ImportGauge(U);
+    Chi = Zero();
+    HeatbathSolver(eta,Chi);       // Chi = B^-1 eta : the ONLY solve
+    std::cout << GridLogMessage << action_name() << " refresh |Chi|^2 = "<< norm2(Chi)<<std::endl;
+  }
+
+  //////////////////////////////////////////////////////
+  // S2 = |B chi|^2  -- matrix multiply only
+  //////////////////////////////////////////////////////
+  virtual RealD S(const GaugeField &U) {
+    BOp.ImportGauge(U);
+
+    FermionField w(BOp.FermionGrid());
+    BOp.M(Chi,w);                  // w = B chi
+    RealD action = norm2(w);
+    return action;
+  }
+
+  //////////////////////////////////////////////////////
+  // dS2 = chi^dag dBdag w + w^dag dB chi ,   w = B chi
+  // NO solves.
+  //////////////////////////////////////////////////////
+  virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+    BOp.ImportGauge(U);
+
+    FermionField w(BOp.FermionGrid());
+    GaugeField force(BOp.GaugeGrid());
+
+    BOp.M(Chi,w);                  // w = B chi
+
+    BOp.MDeriv(force, Chi, w, DaggerYes); dSdU =      force;
+    BOp.MDeriv(force, w, Chi, DaggerNo ); dSdU = dSdU+force;
+
+    dSdU *= -1.0;  // Grid action sign convention (cf TwoFlavourRatio.h)
+  };
+};
+
+NAMESPACE_END(Grid);

Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h

@@ -0,0 +1,264 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h
+
+    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 */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Two flavour ratio with FOUR dimensional pseudofermion, UNpreconditioned
+// (full grid) solves.
+//
+// Companion to TwoFlavourRatioEO4DPseudoFermion.h but with the solver
+// plumbing exposed as LinearFunction<FermionField> objects that already
+// know their operator -- the natural interface for the non-Hermitian
+// multigrid GCR stack (PVdagM), which solves M and Mdag DIRECTLY rather
+// than through SchurRedBlack normal equations.
+//
+// Why: with 5D pseudofermions the squared-operator formulation hands
+// normal-equation solvers (MdagM)^-1 phi AND Mdag^-1 phi from ONE Krylov
+// space; a direct solver must solve twice, halving its per-solve gain.
+// The 4D pseudofermion action needs one M^-1 and one M^-dag solve per
+// force evaluation FOR BOTH solver families, so the direct-solver gain
+// carries through undiluted.  In addition phi4 is Ls-agnostic, so the
+// force can be evaluated with a reduced-Ls operator pair while the
+// accept/reject uses full Ls (inexact force, exact action).
+//
+// Solver slots (all full-grid 5D LinearFunctions, solution overwritten,
+// zero guess imposed internally):
+//   DerivMinvSolver     : x = M^-1 b        (DenOp)
+//   DerivMdagInvSolver  : x = M^-dag b      (DenOp).  For G5R5-hermitian
+//                         actions this may be implemented by the caller as
+//                         G5R5 . DerivMinvSolver . G5R5 -- no adjoint
+//                         multigrid needed.
+//   ActionMinvSolver    : x = M^-1 b        (DenOp, accept/reject tolerance)
+//   HeatbathVinvSolver  : x = V^-1 b        (NumOp)
+//
+// 4D <-> 5D wall maps: the action is S = | P (M^-1 V) Pdag phi4 |^2 where
+// (P,Pdag) MUST be a mutually adjoint pair for S and deriv to be
+// consistent.  Two candidate conventions, selected by solution_walls:
+//   true  : P = P_- psi(0) + P_+ psi(Ls-1)   (solution walls, matches
+//           ExportPhysicalFermionSolution) and Pdag its literal adjoint.
+//   false : P = P_+ psi(0) + P_- psi(Ls-1)   (source walls, Pdag matches
+//           ImportUnphysicalFermion).
+// The heatbath is exact iff [P M^-1 V Pdag][P V^-1 M Pdag] = 1 (the 4D
+// effective-operator composition identity); which convention satisfies it
+// is settled numerically by the refresh test  S == 0.5*|eta4|^2  exactly.
+///////////////////////////////////////////////////////////////////////////////
+template<class Impl>
+class TwoFlavourRatio4DPseudoFermionAction : public Action<typename Impl::GaugeField> {
+public:
+  INHERIT_IMPL_TYPES(Impl);
+
+private:
+  typedef FermionOperator<Impl> FermOp;
+  FermionOperator<Impl> & NumOp;// the basic operator (V)
+  FermionOperator<Impl> & DenOp;// the basic operator (M)
+
+  LinearFunction<FermionField> &DerivMinvSolver;
+  LinearFunction<FermionField> &DerivMdagInvSolver;
+  LinearFunction<FermionField> &ActionMinvSolver;
+  LinearFunction<FermionField> &HeatbathVinvSolver;
+
+  FermionField phi4; // the pseudo fermion field for this trajectory
+
+  int solution_walls; // wall convention for the (P,Pdag) pair; see header
+
+  ////////////////////////////////////////////////////////////////////
+  // The mutually adjoint 4D <-> 5D pair.
+  //   Wall4D    : q4  = P psi5      (extract)
+  //   Wall4DAdj : psi5 = Pdag q4    (insert; literal adjoint of Wall4D)
+  ////////////////////////////////////////////////////////////////////
+  void Wall4D(const FermionField &psi5, FermionField &q4)
+  {
+    int Ls = NumOp.FermionGrid()->_fdimensions[0];
+    FermionField tmp(NumOp.FermionGrid());
+    if ( solution_walls ) {
+      // q4 = P_- psi(0) + P_+ psi(Ls-1)
+      axpby_ssp_pminus(tmp, 0., psi5, 1., psi5, 0, 0);
+      axpby_ssp_pplus (tmp, 1., tmp , 1., psi5, 0, Ls-1);
+    } else {
+      // q4 = P_+ psi(0) + P_- psi(Ls-1)
+      axpby_ssp_pplus (tmp, 0., psi5, 1., psi5, 0, 0);
+      axpby_ssp_pminus(tmp, 1., tmp , 1., psi5, 0, Ls-1);
+    }
+    ExtractSlice(q4, tmp, 0, 0);
+  }
+  void Wall4DAdj(const FermionField &q4, FermionField &psi5)
+  {
+    int Ls = NumOp.FermionGrid()->_fdimensions[0];
+    FermionField tmp(NumOp.FermionGrid());
+    tmp = Zero();
+    InsertSlice(q4, tmp, 0   , 0);
+    InsertSlice(q4, tmp, Ls-1, 0);
+    if ( solution_walls ) {
+      // psi(0) = P_- q4 ; psi(Ls-1) = P_+ q4
+      axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, 0   , 0);
+      axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
+    } else {
+      // psi(0) = P_+ q4 ; psi(Ls-1) = P_- q4
+      axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, 0   , 0);
+      axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
+    }
+    psi5 = tmp;
+  }
+
+public:
+  TwoFlavourRatio4DPseudoFermionAction(FermionOperator<Impl>  &_NumOp,
+				       FermionOperator<Impl>  &_DenOp,
+				       LinearFunction<FermionField> & DMS,
+				       LinearFunction<FermionField> & DMDS,
+				       LinearFunction<FermionField> & AMS,
+				       LinearFunction<FermionField> & HVS,
+				       int _solution_walls = 1
+				       ) : NumOp(_NumOp),
+					   DenOp(_DenOp),
+					   DerivMinvSolver(DMS),
+					   DerivMdagInvSolver(DMDS),
+					   ActionMinvSolver(AMS),
+					   HeatbathVinvSolver(HVS),
+					   phi4(_NumOp.GaugeGrid()),
+					   solution_walls(_solution_walls)
+  {};
+
+  virtual std::string action_name(){return "TwoFlavourRatio4DPseudoFermionAction";}
+
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "["<<action_name()<<"] solution_walls " << solution_walls << std::endl;
+    return sstream.str();
+  }
+
+  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+
+    // P(phi4) = e^{- phi4^dag Beff^dag Beff phi4}   ;  Beff = P M^-1 V Pdag
+    //
+    // NumOp == V
+    // DenOp == M
+    //
+    // Take phi4 = P V^-1 M Pdag eta4   ( = Beff^-1 eta4 by the composition
+    // identity; verified numerically by S == 0.5 |eta4|^2 after refresh )
+    //
+    // P(eta) = e^{- eta^dag eta} ;  e^{-x^2/2 sig^2} => sig^2 = 0.5
+    // so eta enters with width 1/sqrt(2).
+    //
+    RealD scale = std::sqrt(0.5);
+
+    FermionField eta4(NumOp.GaugeGrid());
+    FermionField eta5(NumOp.FermionGrid());
+    FermionField tmp (NumOp.FermionGrid());
+    FermionField phi5(NumOp.FermionGrid());
+
+    gaussian(pRNG,eta4);
+
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    Wall4DAdj(eta4,eta5);            // eta5 = Pdag eta4
+    DenOp.M(eta5,tmp);               // tmp  = M eta5
+    phi5 = Zero();
+    HeatbathVinvSolver(tmp,phi5);    // phi5 = V^-1 M eta5
+    Wall4D(phi5,phi4);               // phi4 = P phi5
+    phi4 = phi4*scale;
+
+    std::cout << GridLogMessage << "4d pf (non-EO) refresh "<< norm2(phi4)<<"\n";
+  };
+
+  //////////////////////////////////////////////////////
+  // S = phi4^dag (Pdag^dag V^dag M^-dag P^dag) (P M^-1 V Pdag) phi4
+  //   = | P M^-1 V Pdag phi4 |^2
+  //////////////////////////////////////////////////////
+  virtual RealD S(const GaugeField &U) {
+
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    FermionField Y4  (NumOp.GaugeGrid());
+    FermionField phi5(NumOp.FermionGrid());
+    FermionField X   (NumOp.FermionGrid());
+    FermionField Y   (NumOp.FermionGrid());
+
+    Wall4DAdj(phi4,phi5);            // phi5 = Pdag phi4
+    NumOp.M(phi5,X);                 // X    = V phi5
+    Y = Zero();
+    ActionMinvSolver(X,Y);           // Y    = M^-1 V phi5
+    Wall4D(Y,Y4);                    // Y4   = P Y
+
+    RealD action = norm2(Y4);
+
+    return action;
+  };
+
+  //////////////////////////////////////////////////////
+  // dS/du = 2 Re [ (M^-dag Pdag w4)^dag dV Pdag phi4 ]
+  //       - 2 Re [ (M^-dag Pdag w4)^dag dM (M^-1 V Pdag phi4) ]
+  // with w4 = P M^-1 V Pdag phi4.
+  // Two first-power solves: one M^-1, one M^-dag.
+  //////////////////////////////////////////////////////
+  virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    FermionField  phi5           (NumOp.FermionGrid());
+    FermionField  Vphi           (NumOp.FermionGrid());
+    FermionField  MinvVphi       (NumOp.FermionGrid());
+    FermionField  w4             (NumOp.GaugeGrid());
+    FermionField  Y              (NumOp.FermionGrid());
+    FermionField  MdagInvPdagW   (NumOp.FermionGrid());
+
+    GaugeField   force(NumOp.GaugeGrid());
+
+    Wall4DAdj(phi4,phi5);              // phi5     = Pdag phi4
+    NumOp.M(phi5,Vphi);                // Vphi     = V phi5
+    MinvVphi = Zero();
+    DerivMinvSolver(Vphi,MinvVphi);    // MinvVphi = M^-1 V phi5
+    std::cout << GridLogMessage << "4d pf (non-EO) deriv solve "<< norm2(MinvVphi)<<"\n";
+
+    // Project onto the physical 4D subspace and back: Y = Pdag P MinvVphi.
+    // Pdag here MUST be the literal adjoint of the P used in S, else the
+    // force is inconsistent with the action.
+    Wall4D(MinvVphi,w4);               // w4 = P MinvVphi
+    Wall4DAdj(w4,Y);                   // Y  = Pdag w4
+
+    MdagInvPdagW = Zero();
+    DerivMdagInvSolver(Y,MdagInvPdagW);   // = M^-dag Pdag w4  (adjoint solve)
+    std::cout << GridLogMessage << "4d pf (non-EO) deriv solve dag "<< norm2(MdagInvPdagW)<<"\n";
+
+    // phi^dag (Pdag' Vdag Mdag^-1 P') (dV) Pdag phi   + h.c.
+    NumOp.MDeriv(force, MdagInvPdagW, phi5, DaggerNo );  dSdU=force;
+    NumOp.MDeriv(force, phi5, MdagInvPdagW, DaggerYes);  dSdU=dSdU+force;
+
+    // - phi^dag ( ... Mdag^-1 ) dM ( M^-1 V ... ) phi  + h.c.
+    DenOp.MDeriv(force, MdagInvPdagW, MinvVphi, DaggerNo );  dSdU=dSdU-force;
+    DenOp.MDeriv(force, MinvVphi, MdagInvPdagW, DaggerYes);  dSdU=dSdU-force;
+
+    dSdU *= -1.0;
+  };
+};
+
+NAMESPACE_END(Grid);

Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h

@@ -0,0 +1,206 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h
+
+    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 */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+///////////////////////////////////////////////////////////////////////////////
+// Two flavour ratio with LEFT-PRECONDITIONED solves.
+//
+// Same action content as TwoFlavourRatio.h:
+//
+//    S = phi^dag V (Mdag M)^-1 Vdag phi     ==>   det[ Mdag M / Vdag V ]
+//
+// (V = NumOp the heavier / Pauli-Villars operator, M = DenOp the lighter),
+// but organised around the composite
+//
+//    F = Vdag M
+//
+// which is the 2-hop-coarsenable operator the non-Hermitian multigrid
+// serves.  Solving M X = b as F X = Vdag b is LEFT PRECONDITIONING by
+// Vdag; the determinant/action layer is the standard quotient, and all
+// novelty is confined to the solver contract.
+//
+// TwoFlavourRatio.h is tied to a normal-equations solver: one (MdagM)^-1
+// solve, then Y = M X gives Mdag^-1 Vdag phi almost free.  The left-
+// preconditioned idiom is DIFFERENT: the chain
+//
+//    b = Vdag phi
+//    z : Fdag z = b          (adjoint  F solve)
+//    Y = V z                 (= Mdag^-1 Vdag phi -- harvested from solve 1)
+//    s = Vdag Y              (= Vdag V z)
+//    X : F X = s             (forward F solve;  X = (MdagM)^-1 Vdag phi)
+//
+// yields Y BEFORE X (so S(U) needs only the adjoint solve), with Y's
+// accuracy independent of the second solve.  Force terms are then the
+// standard four MDeriv insertions of TwoFlavourRatio.
+//
+// Solver slots are LinearFunctions with the F-SOLVE contract (solution
+// overwritten, zero guess imposed internally):
+//    ForwardSolver(b,x)  :  F x    = b
+//    AdjointSolver(b,z)  :  Fdag z = b
+// implemented in production by the multigrid-GCR stack (forward cycle and
+// adjoint cycle); in tests by CG on the composite normal equations.
+// HeatbathSolver(b,x) : x = (Vdag V)^-1 b -- heavy operator, plain CG.
+//
+// Heatbath is exact by operator algebra:  phi = V (VdagV)^-1 Mdag eta
+// ==>  S = | Mdag^-1 Vdag phi |^2 = |eta|^2  (to solver tolerance); the
+// deterministic refresh(U,eta) hook below is the test point.
+//
+// Hasenbusch: nothing requires V to have mass one; any (heavier,lighter)
+// pair works, F(V,M) = Vdag M coarsenable by the same machinery, rungs'
+// solves are F-family (mrhs-batchable, mass-shared coarse space).
+///////////////////////////////////////////////////////////////////////////////
+template<class Impl>
+class TwoFlavourRatioLeftPrecPseudoFermionAction : public Action<typename Impl::GaugeField> {
+public:
+  INHERIT_IMPL_TYPES(Impl);
+
+private:
+  FermionOperator<Impl> & NumOp;// V
+  FermionOperator<Impl> & DenOp;// M
+
+  LinearFunction<FermionField> &DerivForwardSolver;  // F x = b, MD tolerance
+  LinearFunction<FermionField> &DerivAdjointSolver;  // Fdag z = b, MD tolerance
+  LinearFunction<FermionField> &ActionAdjointSolver; // Fdag z = b, accept/reject tolerance
+  LinearFunction<FermionField> &HeatbathSolver;      // (VdagV)^-1 b, heavy op
+
+  FermionField Phi; // the pseudo fermion field for this trajectory
+
+public:
+  TwoFlavourRatioLeftPrecPseudoFermionAction(FermionOperator<Impl>  &_NumOp,
+					     FermionOperator<Impl>  &_DenOp,
+					     LinearFunction<FermionField> & DFS,
+					     LinearFunction<FermionField> & DAS,
+					     LinearFunction<FermionField> & AAS,
+					     LinearFunction<FermionField> & HS
+					     ) : NumOp(_NumOp),
+						 DenOp(_DenOp),
+						 DerivForwardSolver(DFS),
+						 DerivAdjointSolver(DAS),
+						 ActionAdjointSolver(AAS),
+						 HeatbathSolver(HS),
+						 Phi(_NumOp.FermionGrid())
+  {};
+
+  virtual std::string action_name(){return "TwoFlavourRatioLeftPrecPseudoFermionAction";}
+
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
+    return sstream.str();
+  }
+
+  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+    // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} ; phi = Vdag^-1 Mdag eta
+    // e^{-x^2/2 sig^2} => sig^2 = 0.5 ; eta enters with width 1/sqrt(2).
+    RealD scale = std::sqrt(0.5);
+    FermionField eta(NumOp.FermionGrid());
+    gaussian(pRNG,eta);
+    eta = eta * scale;
+    refresh(U,eta);
+  }
+
+  // Deterministic-noise variant (test hook):
+  // after this,  S(U) == norm2(eta)  exactly (to solver tolerance).
+  void refresh(const GaugeField &U, const FermionField &eta) {
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    FermionField tmp(NumOp.FermionGrid());
+    FermionField w  (NumOp.FermionGrid());
+
+    DenOp.Mdag(eta,tmp);           // tmp = Mdag eta
+    w = Zero();
+    HeatbathSolver(tmp,w);         // w   = (VdagV)^-1 Mdag eta
+    NumOp.M(w,Phi);                // Phi = V (VdagV)^-1 Mdag eta = Vdag^-1 Mdag eta
+    std::cout << GridLogMessage << action_name() << " refresh |Phi|^2 = "<< norm2(Phi)<<std::endl;
+  }
+
+  //////////////////////////////////////////////////////
+  // S = phi^dag V (MdagM)^-1 Vdag phi = | Mdag^-1 Vdag phi |^2
+  // ONE adjoint F solve:  Y = V Fdag^-1 Vdag phi = Mdag^-1 Vdag phi
+  //////////////////////////////////////////////////////
+  virtual RealD S(const GaugeField &U) {
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    FermionField b(NumOp.FermionGrid());
+    FermionField z(NumOp.FermionGrid());
+    FermionField Y(NumOp.FermionGrid());
+
+    NumOp.Mdag(Phi,b);             // b = Vdag phi
+    z = Zero();
+    ActionAdjointSolver(b,z);      // Fdag z = b
+    NumOp.M(z,Y);                  // Y = V z = Mdag^-1 Vdag phi
+
+    RealD action = norm2(Y);
+    return action;
+  }
+
+  //////////////////////////////////////////////////////
+  // dS/du = phi^dag dV (MdagM)^-1 Vdag phi
+  //       - phi^dag V (MdagM)^-1 [ Mdag dM + dMdag M ] (MdagM)^-1 Vdag phi
+  //       + phi^dag V (MdagM)^-1 dVdag phi
+  // Identical force insertions to TwoFlavourRatio.h; X and Y from the
+  // left-preconditioned chain (Y harvested from the adjoint solve).
+  //////////////////////////////////////////////////////
+  virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+    NumOp.ImportGauge(U);
+    DenOp.ImportGauge(U);
+
+    FermionField b(NumOp.FermionGrid());
+    FermionField z(NumOp.FermionGrid());
+    FermionField Y(NumOp.FermionGrid());
+    FermionField s(NumOp.FermionGrid());
+    FermionField X(NumOp.FermionGrid());
+
+    GaugeField force(NumOp.GaugeGrid());
+
+    NumOp.Mdag(Phi,b);             // b = Vdag phi
+    z = Zero();
+    DerivAdjointSolver(b,z);       // Fdag z = b
+    NumOp.M(z,Y);                  // Y = V z = Mdag^-1 Vdag phi   (solve-1 harvest)
+    NumOp.Mdag(Y,s);               // s = Vdag V z
+    X = Zero();
+    DerivForwardSolver(s,X);       // F X = s  ==>  X = (MdagM)^-1 Vdag phi
+
+    // phi^dag V (MdagM)^-1 dVdag phi
+    NumOp.MDeriv(force , X, Phi, DaggerYes);  dSdU  =      force;
+    // phi^dag dV (MdagM)^-1 Vdag phi
+    NumOp.MDeriv(force , Phi, X, DaggerNo );  dSdU  = dSdU+force;
+    //  - phi^dag V (MdagM)^-1 Mdag dM (MdagM)^-1 Vdag phi
+    //  - phi^dag V (MdagM)^-1 dMdag M (MdagM)^-1 Vdag phi
+    DenOp.MDeriv(force, Y, X, DaggerNo );     dSdU  = dSdU-force;
+    DenOp.MDeriv(force, X, Y, DaggerYes);     dSdU  = dSdU-force;
+
+    dSdU *= -1.0;
+  };
+};
+
+NAMESPACE_END(Grid);

examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc

@@ -0,0 +1,663 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.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.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+
+// MultiRHS (valence) THREE-level multigrid for PVdagM with a DENSE, EXACT,
+// non-iterative coarse-coarse bottom -- the LIBRARY-CLASS successor of
+// Example_pvdagm_mrhs_3level_dense.cc, which is FROZEN as the regression
+// baseline / champion-provenance artifact (21.7 s/RHS at BLOCK=2.2.3.3,
+// BLOCK2=8.4.2.4, nb60, CSO3/FSO6/CST0.04 on 36 Frontier nodes).
+//
+// The dense bottom is now Grid/algorithms/multigrid/DenseCoarseMatrix.h:
+//  - stencil -> dense DIRECT import (no probe assembly: rows are local data)
+//    + IMPORT CERTIFICATE  (DENSE_IMPORT_SIGN=-1 flips convention, no rebuild)
+//  - split-K apply via GridBLAS.gemmBatched with explicit leading dimensions
+//    (DENSE_SPLITK chunks, default 32) -- the fig-11 software split-K
+//  - deviceVector / GridBLAS throughout the apply: platform-agnostic
+//
+// INTERCHANGE: same SLAB_FILE per-rank format as the frozen example (stem MUST
+// encode cfg/mass/blocking/nbasis; the header guards only N/nrows/nbasis) and
+// the same env-var set, so existing sbatch scripts drive either binary.
+//
+// A/B acceptance (old binary = control):
+//   slab-cached  : outer counts match EXACTLY (identical apply data; split-K
+//                  changes only fp32 reduction order); wall delta = split-K gain.
+//   fresh setup  : outer equal-or-+-1 (import vs probe = rounding); VERIFY
+//                  ~7e-4 both; setup delta = import gain (~93 s probe retired).
+//
+// Level structure, solvers, and tuning knobs are UNCHANGED from the frozen
+// example.  Env: MASS SUBSPACE_FILE NRHS BLOCK BLOCK2 FineSmootherShift/Order
+// CoarseSmootherShift/Nstep CoarseSolverTol/Order DENSE_CC DENSE_CC_CHECK
+// DENSE_SPLITK DENSE_DEVICE_SUM DENSE_IMPORT_SIGN DENSE_APPLY_PROFILE
+// L3_TOL L3_MAXIT L3_NSTEP OuterMmax OuterNstep OuterTol
+
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
+#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
+
+#include <memory>
+
+using namespace std;
+using namespace Grid;
+
+RealD FineSmootherShift    = 0.1;
+int   FineSmootherOrder    = 16;
+RealD CoarseSmootherShift  = 0.1;
+int   CoarseSmootherNstep  = 4;
+RealD CoarseSolverTol      = 0.03;
+int   CoarseSolverOrder    = 200;
+RealD L3Tol                = 2.5e-1;
+int   L3MaxIt              = 50;
+int   L3Nstep              = 50;
+RealD OuterTol             = 1.0e-8;
+int   OuterMmax            = 8;
+int   OuterNstep           = 8;
+int   Nrhs                 = 12;
+int   UseDenseCC           = 1;
+RealD mass                 = 0.00078;
+
+void ParseEnvironment(void)
+{
+  if(getenv("MASS"))               mass               = atof(getenv("MASS"));
+  if(getenv("FineSmootherShift"))  FineSmootherShift  = atof(getenv("FineSmootherShift"));
+  if(getenv("FineSmootherOrder"))  FineSmootherOrder  = atoi(getenv("FineSmootherOrder"));
+  if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
+  if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
+  if(getenv("CoarseSolverTol"))    CoarseSolverTol    = atof(getenv("CoarseSolverTol"));
+  if(getenv("CoarseSolverOrder"))  CoarseSolverOrder  = atoi(getenv("CoarseSolverOrder"));
+  if(getenv("L3_TOL"))             L3Tol              = atof(getenv("L3_TOL"));
+  if(getenv("L3_MAXIT"))           L3MaxIt            = atoi(getenv("L3_MAXIT"));
+  if(getenv("L3_NSTEP"))           L3Nstep            = atoi(getenv("L3_NSTEP"));
+  if(getenv("OuterTol"))           OuterTol           = atof(getenv("OuterTol"));
+  if(getenv("OuterMmax"))          OuterMmax          = atoi(getenv("OuterMmax"));
+  if(getenv("OuterNstep"))         OuterNstep         = atoi(getenv("OuterNstep"));
+  if(getenv("NRHS"))               Nrhs               = atoi(getenv("NRHS"));
+  if(getenv("DENSE_CC"))           UseDenseCC         = atoi(getenv("DENSE_CC"));
+
+  std::cout << GridLogMessage << "PARAM: MASS               " << mass               << std::endl;
+  std::cout << GridLogMessage << "PARAM: NRHS               " << Nrhs               << std::endl;
+  std::cout << GridLogMessage << "PARAM: DENSE_CC           " << UseDenseCC         << std::endl;
+  std::cout << GridLogMessage << "PARAM: FineSmootherShift  " << FineSmootherShift  << std::endl;
+  std::cout << GridLogMessage << "PARAM: FineSmootherOrder  " << FineSmootherOrder  << std::endl;
+  std::cout << GridLogMessage << "PARAM: CoarseSmootherShift" << CoarseSmootherShift<< std::endl;
+  std::cout << GridLogMessage << "PARAM: CoarseSmootherNstep" << CoarseSmootherNstep<< std::endl;
+  std::cout << GridLogMessage << "PARAM: CoarseSolverTol    " << CoarseSolverTol    << std::endl;
+  std::cout << GridLogMessage << "PARAM: CoarseSolverOrder  " << CoarseSolverOrder  << std::endl;
+  std::cout << GridLogMessage << "PARAM: L3_TOL             " << L3Tol              << std::endl;
+  std::cout << GridLogMessage << "PARAM: OuterMmax          " << OuterMmax          << std::endl;
+  std::cout << GridLogMessage << "PARAM: OuterNstep         " << OuterNstep         << std::endl;
+}
+
+template <class Field>
+void saveSubspace(std::vector<Field> &subspace, std::string const fname){
+#ifdef HAVE_LIME
+  Grid::emptyUserRecord record;
+  Grid::ScidacWriter SW(subspace[0].Grid()->IsBoss());
+  SW.open(fname);
+  for (int k = 0; k < (int)subspace.size(); k++) SW.writeScidacFieldRecord(subspace[k], record);
+  SW.close();
+#endif
+}
+template <class Field>
+void loadSubspace(std::vector<Field> &subspace, std::string const fname){
+#ifdef HAVE_LIME
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader SR;
+  SR.open(fname);
+  for (int k = 0; k < (int)subspace.size(); k++) SR.readScidacFieldRecord(subspace[k], record);
+  SR.close();
+#endif
+}
+
+//////////////////////////////////////////////////////////////////////
+// A = PV^dag M (non-Hermitian), and shifted variant for smoothers.
+//////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class PVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+template<class Matrix,class Field>
+class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat; Matrix &_PV;
+public:
+  RealD shift;
+  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
+  void OpDiag (const Field &in, Field &out) { assert(0); }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); };
+  void Op     (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; }
+  void AdjOp  (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+// Generic shift wrapper (for the coarse-level smoother on the 6D mrhs coarse operator).
+template<class Field>
+class ShiftedLinearOperator : public LinearOperatorBase<Field> {
+  LinearOperatorBase<Field> &_Op; RealD shift;
+public:
+  ShiftedLinearOperator(RealD _shift, LinearOperatorBase<Field> &Op) : _Op(Op), shift(_shift) {}
+  void OpDiag  (const Field &in, Field &out) { assert(0); }
+  void OpDir   (const Field &in, Field &out,int dir,int disp) { assert(0); }
+  void OpDirAll (const Field &in, std::vector<Field> &out) { assert(0); }
+  void Op      (const Field &in, Field &out) { _Op.Op(in,out);    out = out + shift*in; }
+  void AdjOp   (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
+  void HermOp  (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
+};
+
+//////////////////////////////////////////////////////////////////////
+// Dense CC solve on the PACKED 6D mrhs coarse-coarse field: drop-in
+// for the L3 PGCR, delegating to the library DenseCoarseMatrix.
+//////////////////////////////////////////////////////////////////////
+template<class DenseType, class CoarseCoarseField>
+class MrhsDenseCCSolve : public LinearFunction<CoarseCoarseField> {
+public:
+  DenseType &_Dense;
+  GridBase *_CoarseCoarse5d;
+  int _nrhs;
+  MrhsDenseCCSolve(DenseType &D, GridBase *cc5d, int nrhs)
+    : _Dense(D), _CoarseCoarse5d(cc5d), _nrhs(nrhs) {}
+  using LinearFunction<CoarseCoarseField>::operator();
+  virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){
+    if ( getenv("DENSE_CC_CHECK") ) {
+      // Audit path: per-rhs 5D unpack so ApplyBatch can run the _Op defect
+      // check per rhs.  ~50ms/call of slice/split overhead -- audit only.
+      CoarseCoarseField tmp(in.Grid());
+      tmp = in;
+      std::vector<CoarseCoarseField> split_in (_nrhs,_CoarseCoarse5d);
+      std::vector<CoarseCoarseField> split_out(_nrhs,_CoarseCoarse5d);
+      for(int r=0;r<_nrhs;r++) ExtractSliceFast(split_in[r], tmp, r, 0);
+      _Dense.ApplyBatch(split_in, split_out);
+      for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0);
+    } else {
+      _Dense.ApplyBatch6D(in, out, _nrhs);
+    }
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from the
+// frozen Example_pvdagm_mrhs_3level_dense.cc)
+//////////////////////////////////////////////////////////////////////
+template<class Field>
+class MrhsLinearFunction {
+public:
+  virtual void operator()(std::vector<Field> &in, std::vector<Field> &out) = 0;
+};
+
+template<class Field>
+class MrhsPGCRNonHermitian {
+public:
+  RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level;
+  int ZeroGuess = 0; int FirstCycle = 0;   // caller contract: zero guess => first-cycle r0 = src
+  std::string name = "Level 1";
+  LinearOperatorBase<Field> &Linop;
+  MrhsLinearFunction<Field> &Preconditioner;
+  void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; }
+  void Name(std::string n){ name = n; }
+  void SetZeroGuess(int z){ ZeroGuess=z; }
+  MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,MrhsLinearFunction<Field> &Prec,int _mmax,int _nstep)
+    : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; }
+  static RealD vnorm2(std::vector<Field> &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; }
+  static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; }
+  static void vaxpy(std::vector<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); }
+  void vOp(std::vector<Field> &in,std::vector<Field> &out){ for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
+  void operator()(std::vector<Field> &src,std::vector<Field> &psi){
+    RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
+    ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq;
+    std::vector<Field> r(nrhs,grid);
+    GridStopWatch T; T.Start(); steps=0; FirstCycle=1;
+    for(int k=0;k<MaxIterations;k++){
+      cp=GCRnStep(src,psi,rsq);
+      std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR("<<mmax<<","<<nstep<<") "<<steps<<" steps cp = "<<cp<<" target "<<rsq<<std::endl;
+      if(cp<rsq){
+        T.Stop(); vOp(psi,r); for(int rr=0;rr<nrhs;rr++) axpy(r[rr],-1.0,src[rr],r[rr]);
+        RealD tr=vnorm2(r);
+        std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR: Converged on iteration "<<steps
+                 <<" computed residual "<<std::sqrt(cp/ssq)<<" true residual "<<std::sqrt(tr/ssq)<<" target "<<Tolerance<<std::endl;
+        std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR Time elapsed: Total "<<T.Elapsed()<<std::endl;
+        for(int rr=0;rr<nrhs;rr++){ RealD rn=std::sqrt(norm2(r[rr])/norm2(src[rr])); std::cout<<GridLogMessage<<"MrhsPGCR per-rhs true residual["<<rr<<"] = "<<rn<<std::endl; }
+        return;
+      }
+    }
+    std::cout<<GridLogMessage<<"MrhsPGCR: did not converge"<<std::endl;
+  }
+  RealD GCRnStep(std::vector<Field> &src,std::vector<Field> &psi,RealD rsq){
+    RealD cp; ComplexD a,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid();
+    std::vector<Field> r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid);
+    std::vector< std::vector<Field> > q(mmax,std::vector<Field>(nrhs,grid));
+    std::vector< std::vector<Field> > p(mmax,std::vector<Field>(nrhs,grid));
+    std::vector<RealD> qq(mmax);
+    std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR nStep("<<nstep<<")"<<std::endl;
+    if (ZeroGuess && FirstCycle) { for(int rr=0;rr<nrhs;rr++){ psi[rr]=Zero(); r[rr]=src[rr]; } }
+    else                         { vOp(psi,Az); for(int rr=0;rr<nrhs;rr++) r[rr]=src[rr]-Az[rr]; }
+    FirstCycle=0;
+    Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
+    p[0]=z; q[0]=Az; qq[0]=zAAz; cp=vnorm2(r);
+    for(int k=0;k<nstep;k++){
+      steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
+      rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
+      vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
+      std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"]  resid "<<cp<<" target "<<rsq<<std::endl;
+      if((k==nstep-1)||(cp<rsq)) return cp;
+      Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
+      q[peri_kp]=Az; p[peri_kp]=z;
+      int northog=((kp)>(mmax-1))?(mmax-1):(kp);
+      for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0);
+        b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back];
+        vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); }
+      qq[peri_kp]=vnorm2(q[peri_kp]);
+    }
+    GRID_ASSERT(0); return cp;
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// L2->L3 mrhs V-cycle: LinearFunction on the 6D mrhs COARSE field.
+// The coarse-coarse solve slot takes EITHER the dense mrhs solve
+// (DENSE_CC=1) or the L3 PGCR (DENSE_CC=0).
+//////////////////////////////////////////////////////////////////////
+template<class CoarseField, class CoarseCoarseField>
+class MrhsCoarseThreeLevelPrec : public LinearFunction<CoarseField> {
+public:
+  LinearOperatorBase<CoarseField>          &_CoarseOp;          // mrhs coarse op (6D)
+  LinearFunction<CoarseField>              &_CoarseSmoother;    // shifted 6D coarse smoother
+  MultiRHSBlockProject<CoarseField>        &_Projector;         // L2->L3 (vector-based)
+  LinearFunction<CoarseCoarseField>        &_CoarseCoarseSolve; // L3 solve (6D cc)
+  GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs;
+  int _nrhs;
+
+  MrhsCoarseThreeLevelPrec(LinearOperatorBase<CoarseField> &CoarseOp,
+                           LinearFunction<CoarseField> &CoarseSmoother,
+                           MultiRHSBlockProject<CoarseField> &Projector,
+                           LinearFunction<CoarseCoarseField> &CoarseCoarseSolve,
+                           GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs)
+    : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector),
+      _CoarseCoarseSolve(CoarseCoarseSolve),
+      _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {}
+
+  using LinearFunction<CoarseField>::operator();
+  virtual void operator()(const CoarseField &in, CoarseField &out) {
+    int nrhs=_nrhs; double t;
+    CoarseField vec1(in.Grid());
+    CoarseField vec2(in.Grid());
+
+    // trivial pre-smoother
+    out = in;
+
+    // residual (6D coarse)
+    _CoarseOp.Op(out,vec1);  sub(vec1,in,vec1);
+
+    // restrict: unpack 6D coarse -> vector<Coarse> -> blockProject -> vector<CoarseCoarse> -> pack 6D cc
+    std::vector<CoarseField>       csplit(nrhs,_Coarse5d);
+    std::vector<CoarseCoarseField> ccsplit(nrhs,_CoarseCoarse5d);
+    CoarseCoarseField CCsrc(_CoarseCoarseMrhs);
+    CoarseCoarseField CCsol(_CoarseCoarseMrhs);
+
+    t=-usecond();
+    for(int r=0;r<nrhs;r++) ExtractSliceFast(csplit[r],vec1,r,0);
+    _Projector.blockProject(csplit,ccsplit);
+    for(int r=0;r<nrhs;r++) InsertSliceFast(ccsplit[r],CCsrc,r,0);
+    t+=usecond();
+    std::cout<<GridLogMessage<<"L2->L3 restrict took "<<t/1000.0<<"ms"<<std::endl;
+
+    // L3 solve (dense mrhs GEMM, or PGCR)
+    t=-usecond();
+    CCsol=Zero();
+    _CoarseCoarseSolve(CCsrc,CCsol);
+    t+=usecond();
+    std::cout<<GridLogMessage<<"L3 coarse-coarse solve took "<<t/1000.0<<"ms  ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
+
+    // prolong: unpack 6D cc -> blockPromote -> pack 6D coarse; add correction
+    t=-usecond();
+    for(int r=0;r<nrhs;r++) ExtractSliceFast(ccsplit[r],CCsol,r,0);
+    _Projector.blockPromote(csplit,ccsplit);
+    for(int r=0;r<nrhs;r++) InsertSliceFast(csplit[r],vec1,r,0);
+    add(out,out,vec1);
+    t+=usecond();
+    std::cout<<GridLogMessage<<"L2->L3 prolong took "<<t/1000.0<<"ms"<<std::endl;
+
+    // residual + coarse smoother (6D coarse)
+    _CoarseOp.Op(out,vec1);  sub(vec1,in,vec1);
+    vec2=Zero();
+    _CoarseSmoother(vec1,vec2);
+    add(out,out,vec2);
+  }
+};
+
+//////////////////////////////////////////////////////////////////////
+// L1->L2 mrhs V-cycle (verbatim from the frozen example)
+//////////////////////////////////////////////////////////////////////
+template<class FineField, class MrhsCoarseVector, class FineSmoother>
+class MrhsTwoLevelMG : public MrhsLinearFunction<FineField> {
+public:
+  typedef MrhsCoarseVector CoarseVector;
+  LinearOperatorBase<FineField>   &_FineOperator;
+  FineSmoother                    &_PostSmoother;
+  MultiRHSBlockProject<FineField> &_Projector;
+  LinearFunction<CoarseVector>    &_CoarseSolve;
+  GridBase *_CoarseGrid, *_CoarseGridMrhs;
+  MrhsTwoLevelMG(LinearOperatorBase<FineField> &FineOp, FineSmoother &Post,
+                 MultiRHSBlockProject<FineField> &Projector, LinearFunction<CoarseVector> &CoarseSolve,
+                 GridBase *CoarseGrid, GridBase *CoarseGridMrhs)
+    : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve),
+      _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){}
+  virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){
+    int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t;
+    std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid);
+    for(int r=0;r<nrhs;r++) out[r]=in[r];
+    for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    std::vector<CoarseVector> Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid);
+    CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs);
+    t=-usecond();
+    _Projector.blockProject(vec1,Csrc_split);
+    for(int r=0;r<nrhs;r++) InsertSliceFast(Csrc_split[r],CsrcMrhs,r,0);
+    t+=usecond(); std::cout<<GridLogMessage<<"Mrhs project+pack took "<<t/1000.0<<"ms"<<std::endl;
+    t=-usecond(); CsolMrhs=Zero(); _CoarseSolve(CsrcMrhs,CsolMrhs); t+=usecond();
+    std::cout<<GridLogMessage<<"Mrhs coarse solve took "<<t/1000.0<<"ms  ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
+    t=-usecond();
+    for(int r=0;r<nrhs;r++) ExtractSliceFast(Csol_split[r],CsolMrhs,r,0);
+    _Projector.blockPromote(vec1,Csol_split);
+    for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
+    t+=usecond(); std::cout<<GridLogMessage<<"Mrhs unpack+promote took "<<t/1000.0<<"ms"<<std::endl;
+    for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
+    t=-usecond();
+    for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); }
+    t+=usecond(); std::cout<<GridLogMessage<<"Mrhs post-smooth took "<<t/1000.0<<"ms  ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
+  }
+};
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+  ParseEnvironment();
+
+  const int Ls=24; RealD M5=1.8, b=1.5, c=0.5;
+  const int nbasis=60; const int nrhs=Nrhs;
+  GRID_ASSERT(nrhs % vComplex::Nsimd() == 0);
+
+  std::vector<int> lat_size {48,48,48,96};
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(lat_size, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  // Level 1 blocking (default 2^4)
+  Coordinate clatt = lat_size;
+  Coordinate Block({2,2,2,2});
+  if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; }
+  std::cout << GridLogMessage << "Block  " << Block  << "  coarse lattice       " << clatt << std::endl;
+
+  // Level 2 blocking: SUPERCOARSE default 8,4,3,6 -> CC [3,6,8,8], the dense floor.
+  Coordinate cclatt = clatt;
+  Coordinate Block2({8,4,3,6});
+  if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); }
+  for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; }
+  std::cout << GridLogMessage << "Block2 " << Block2 << "  coarse-coarse lattice " << cclatt << std::endl;
+
+  GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+  GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
+
+  // 6D mrhs grids: rhs is dim 0, SIMD across rhs
+  Coordinate mpi=GridDefaultMpi();
+  Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  Coordinate rhSimd({vComplex::Nsimd(),1,1,1,1,1});
+  Coordinate rhLatt ({nrhs,1,clatt[0], clatt[1], clatt[2], clatt[3]});
+  Coordinate rhLatt2({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
+  GridCartesian *CoarseMrhs       = new GridCartesian(rhLatt, rhSimd,rhMpi);
+  GridCartesian *CoarseCoarseMrhs = new GridCartesian(rhLatt2,rhSimd,rhMpi);
+
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  LatticeGaugeField Umu(UGrid);
+  std::cout << GridLogMessage << "Reading gauge field" << std::endl;
+  FieldMetaData header;
+  std::string file("/ccs/home/poare/ckpoint_lat.1000");
+  NerscIO::readConfiguration(Umu,header,file);
+
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+  MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
+
+  typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD>        PVdagM_t;
+  typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
+
+  // Level 1 tensor types
+  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>       LittleDiracOperator;
+  typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>  MrhsLittleDiracOperator;
+  typedef LittleDiracOperator::CoarseVector                                CoarseVector;
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis>                  Subspace;
+
+  // Level 2 tensor types (coarsening deepens the nest by one iScalar)
+  typedef CoarseVector::vector_object                                      CoarseSiteObj;
+  typedef iScalar<vTComplex>                                              vTTComplex;
+  typedef GeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis>          LittleDiracOperatorL2;
+  typedef MultiGeneralCoarsenedMatrix<CoarseSiteObj,vTTComplex,nbasis>     MrhsLittleDiracOperatorL2;
+  typedef LittleDiracOperatorL2::CoarseVector                              CoarseCoarseVector;
+  typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis>                     SubspaceL2;
+
+  // The library dense bottom over the L2 coarse operator
+  typedef DenseCoarseMatrix<CoarseSiteObj,vTTComplex,nbasis>               DenseCC_t;
+
+  PVdagM_t        PVdagM(Ddwf,Dpv);
+  ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
+
+  NextToNearestStencilGeometry5D geom (Coarse5d);
+  NextToNearestStencilGeometry5D geom2(CoarseCoarse5d);
+
+  //////////////////////////////////////////////////////////////////////
+  // Subspace: load RAW (no Orthogonalise!), or generate.
+  //////////////////////////////////////////////////////////////////////
+  std::string subspace_file = "/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/subspace_nb"
+                            + std::to_string(nbasis) + ".scidac";
+  if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE"));
+  uint64_t file_exists=0;
+  if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; }
+  UGrid->GlobalSum(file_exists);
+
+  const int cb=0;
+  Subspace AggregatesGCR(Coarse5d,FGrid,cb);
+  if ( file_exists ){
+    std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl;
+    loadSubspace(AggregatesGCR.subspace, subspace_file);
+  } else {
+    std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl;
+    AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis);
+    saveSubspace(AggregatesGCR.subspace, subspace_file);
+  }
+
+  // RAW copy of the fine null vectors BEFORE CoarsenOperator block-orthonormalises in place.
+  std::vector<LatticeFermionD> rawNull(nbasis,FGrid);
+  for(int k=0;k<nbasis;k++) rawNull[k]=AggregatesGCR.subspace[k];
+
+  //////////////////////////////////////////////////////////////////////
+  // Coarsen L1->L2 and L2->L3 with SINGLE-RHS machinery; import to mrhs via
+  // CopyMatrix.  The L2 (coarse-coarse) single-RHS operator is HOISTED to
+  // main scope: DenseCoarseMatrix imports its stencil and uses its M for
+  // certificates, so it must stay alive for the whole run.
+  //////////////////////////////////////////////////////////////////////
+  MrhsLittleDiracOperator   mrhsLittleDiracOpPV(geom,  CoarseMrhs);
+  MrhsLittleDiracOperatorL2 mrhsLittleDiracOpL2(geom2, CoarseCoarseMrhs);
+  MultiRHSBlockProject<LatticeFermionD> MrhsProjector;
+  MultiRHSBlockProject<CoarseVector>    MrhsProjectorL2;
+
+  LittleDiracOperatorL2 LittleDiracOpL2(geom2,Coarse5d,CoarseCoarse5d);
+  NonHermitianLinearOperator<LittleDiracOperatorL2,CoarseCoarseVector> LinOpCC5d(LittleDiracOpL2);
+
+  {
+    // --- L1->L2 single-RHS coarse operator (scoped: its padded _A is the memory peak) ---
+    LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
+    LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesGCR);   // orthonormalises AggregatesGCR.subspace in place
+    mrhsLittleDiracOpPV.CopyMatrix(LittleDiracOpPV);
+    MrhsProjector.Allocate(nbasis,FGrid,Coarse5d);
+    MrhsProjector.ImportBasis(AggregatesGCR.subspace);
+    NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
+
+    // --- psi_coarse = P^dag (RAW fine null) -> Galerkin images, NOT e_k ---
+    std::vector<CoarseVector> psi_coarse(nbasis,Coarse5d);
+    for(int k=0;k<nbasis;k++) AggregatesGCR.ProjectToSubspace(psi_coarse[k], rawNull[k]);
+    rawNull.clear(); rawNull.shrink_to_fit();
+    {
+      RealD s2=0.0;
+      for(int i=0;i<nbasis;i++) for(int j=0;j<nbasis;j++){
+        ComplexD sij=TensorRemove(innerProduct(psi_coarse[i],psi_coarse[j]));
+        ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0)); s2+=real(d)*real(d)+imag(d)*imag(d);
+      }
+      std::cout<<GridLogMessage<<"GUARD: ||<psi_coarse|psi_coarse> - I||_F = "<<std::sqrt(s2)
+               <<"   (~0.23 good; ~sqrt(N_coarse)="<<std::sqrt((double)Coarse5d->gSites())<<" = e_k leak)"<<std::endl;
+    }
+
+    // --- L2->L3 single-RHS coarsening ---
+    SubspaceL2 AggregatesL2(CoarseCoarse5d,Coarse5d,cb);
+    for(int k=0;k<nbasis;k++) AggregatesL2.subspace[k]=psi_coarse[k];
+    LittleDiracOpL2.CoarsenOperator(LinOpCoarse, AggregatesL2);
+    mrhsLittleDiracOpL2.CopyMatrix(LittleDiracOpL2);
+    MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d);
+    MrhsProjectorL2.ImportBasis(AggregatesL2.subspace);
+
+    // --- guard psi_cc ---
+    {
+      std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d);
+      for(int k=0;k<nbasis;k++) AggregatesL2.ProjectToSubspace(psi_cc[k], psi_coarse[k]);
+      RealD s2=0.0;
+      for(int i=0;i<nbasis;i++) for(int j=0;j<nbasis;j++){
+        ComplexD sij=TensorRemove(innerProduct(psi_cc[i],psi_cc[j]));
+        ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0)); s2+=real(d)*real(d)+imag(d)*imag(d);
+      }
+      std::cout<<GridLogMessage<<"GUARD: ||<psi_cc|psi_cc> - I||_F = "<<std::sqrt(s2)
+               <<"   (~0.23 good; ~sqrt(N_cc)="<<std::sqrt((double)CoarseCoarse5d->gSites())<<" = e_k leak)"<<std::endl;
+    }
+  } // single-RHS FINE op + padded _A + AggregatesL2 + psi_coarse freed here
+
+  NonHermitianLinearOperator<MrhsLittleDiracOperator,CoarseVector>         mrhsLinOpCoarse(mrhsLittleDiracOpPV);
+  NonHermitianLinearOperator<MrhsLittleDiracOperatorL2,CoarseCoarseVector> mrhsLinOpCC(mrhsLittleDiracOpL2);
+
+  //////////////////////////////////////////////////////////////////////
+  // DENSE coarse-coarse bottom: the LIBRARY class, constructed AFTER the
+  // fine coarsening frees its memory peak.  Imports the stencil of the
+  // hoisted single-RHS LittleDiracOpL2 directly (no probing).
+  //////////////////////////////////////////////////////////////////////
+  std::unique_ptr<DenseCC_t> DenseCC;
+  std::unique_ptr<MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector>> MrhsDenseCC;
+  if (UseDenseCC) {
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+    std::cout << GridLogMessage << " Dense CC inverse setup (library DenseCoarseMatrix)" << std::endl;
+    std::cout << GridLogMessage << "**********************************************" << std::endl;
+    DenseCC.reset(new DenseCC_t(LittleDiracOpL2, CoarseCoarse5d));
+    MrhsDenseCC.reset(new MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector>(*DenseCC, CoarseCoarse5d, nrhs));
+  }
+
+  //////////////////////////////////////////////////////////////////////
+  // Solvers, innermost first.
+  //////////////////////////////////////////////////////////////////////
+  TrivialPrecon<CoarseVector>       simpleC;
+  TrivialPrecon<CoarseCoarseVector> simpleCC;
+  TrivialPrecon<LatticeFermionD>    simple_fine;
+
+  // L3 (coarse-coarse) iterative solve: PGCR on the 6D cc operator (DENSE_CC=0 branch)
+  PrecGeneralisedConjugateResidualNonHermitian<CoarseCoarseVector>
+    L3PGCR(L3Tol,L3MaxIt,mrhsLinOpCC,simpleCC,L3Nstep,L3Nstep);
+  L3PGCR.Level(3);
+  L3PGCR.Name("CCouter");
+
+  LinearFunction<CoarseCoarseVector> *ccSolve;
+  if (UseDenseCC) ccSolve = MrhsDenseCC.get();
+  else            ccSolve = &L3PGCR;
+
+  // L2 coarse smoother: shifted 6D coarse op, fixed nstep
+  ShiftedLinearOperator<CoarseVector> ShiftedMrhsCoarse(CoarseSmootherShift, mrhsLinOpCoarse);
+  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+    CoarseSmootherGCR(0.01,1,ShiftedMrhsCoarse,simpleC,CoarseSmootherNstep,CoarseSmootherNstep);
+  CoarseSmootherGCR.Level(2);
+  CoarseSmootherGCR.Name("Csmoother");
+  CoarseSmootherGCR.SetZeroGuess(1);   // caller zeroes vec2: skip r0 apply every L2 iteration
+
+  // L2->L3 V-cycle preconditioner (operates on 6D coarse field)
+  MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
+    L2to3Precon(mrhsLinOpCoarse, CoarseSmootherGCR, MrhsProjectorL2, *ccSolve,
+                Coarse5d, CoarseCoarse5d, CoarseCoarseMrhs, nrhs);
+
+  // L2 coarse solve: PGCR on 6D coarse op, preconditioned by the L2->L3 V-cycle
+  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
+    L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, mrhsLinOpCoarse, L2to3Precon, 16, 16);
+  L2PGCR.Level(2);
+  L2PGCR.Name("Couter");
+  L2PGCR.SetZeroGuess(1);              // caller zeroes CsolMrhs; restarts still recompute r
+
+  // Fine smoother (per-rhs, looped in the L1->L2 V-cycle)
+  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD>
+    SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherOrder,FineSmootherOrder);
+  SmootherGCR.Level(1);
+  SmootherGCR.Name("Fsmoother");
+  SmootherGCR.SetZeroGuess(1);         // caller zeroes vec2[r]: saves 12 fine mults/outer
+
+  // L1->L2 V-cycle (fine); its coarse solve is the three-level L2PGCR
+  typedef PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> FineSmoother_t;
+  MrhsTwoLevelMG<LatticeFermionD,CoarseVector,FineSmoother_t>
+    ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CoarseMrhs);
+
+  // Outer mrhs solve
+  MrhsPGCRNonHermitian<LatticeFermionD>
+    L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
+  L1PGCR.Level(1);
+  L1PGCR.Name("Fouter");
+  L1PGCR.SetZeroGuess(1);              // sol[r]=Zero() below; restarts recompute r as always
+
+  //////////////////////////////////////////////////////////////////////
+  // Sources and solve
+  //////////////////////////////////////////////////////////////////////
+  std::vector<LatticeFermionD> src(nrhs,FGrid), sol(nrhs,FGrid);
+  for(int r=0;r<nrhs;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
+
+  std::cout << GridLogMessage << "**********************************************" << std::endl;
+  std::cout << GridLogMessage << " MultiRHS THREE-level solve (DenseCoarseMatrix bottom): " << nrhs << " RHS " << std::endl;
+  std::cout << GridLogMessage << "**********************************************" << std::endl;
+
+  GridStopWatch w; w.Start();
+  L1PGCR(src,sol);
+  w.Stop();
+  std::cout << GridLogMessage << "MultiRHS 3-level dense solve total " << w.Elapsed()
+            << "  (per RHS: " << w.useconds()/1.0e6/nrhs << " s)" << std::endl;
+
+  { LatticeFermionD Ax(FGrid); RealD worst=0.0;
+    for(int r=0;r<nrhs;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
+      RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
+      std::cout << GridLogMessage << "FINAL: rhs["<<r<<"] true residual = " << rn << std::endl;
+      worst=std::max(worst,rn); }
+    std::cout << GridLogMessage << "FINAL: worst-case residual = " << worst << std::endl;
+  }
+
+  std::cout << GridLogMessage << "Done" << std::endl;
+  Grid_finalize();
+  return 0;
+}

tests/debug/Test_batched_blas.cc

@@ -0,0 +1,256 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_batched_blas.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 */
+
+//
+// Unit tests for the blas-layer dense primitives:
+//
+//  T1 : GridBLASInverse::inverseBatched (ComplexF and ComplexD) --
+//       || A A^{-1} - I ||_max over a random well-conditioned batch.
+//       On a CPU build this exercises the Eigen reference (the semantic
+//       oracle); the SAME binary exercises rocSOLVER/cuBLAS/oneMKL on a
+//       device build.
+//  T2 : explicit-leading-dimension gemmBatched == SOFTWARE SPLIT-K in
+//       miniature.  Y = op(A,T).X computed (a) in one compact batch=1 call
+//       and (b) as NK K-chunks by POINTER OFFSET with lda = full K, partials
+//       reduced on the host.  (a)==(b) validates the strided overload AND
+//       the production dense-slab apply pattern (arXiv:2409.03904 fig 11).
+//
+// Hard asserts throughout: regression gate for Grid/algorithms/blas.
+//
+#include <Grid/Grid.h>
+#include <Grid/algorithms/blas/BatchedInverse.h>
+
+#include <random>
+
+using namespace std;
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  GridBLAS        blas;
+  GridBLASInverse inverse;
+
+  std::mt19937 rng(12345);
+  std::uniform_real_distribution<double> dist(-1.0,1.0);
+
+  ////////////////////////////////////////////////////////////////
+  // T1a : batched inversion, ComplexF
+  ////////////////////////////////////////////////////////////////
+  {
+    const int64_t N = 64;
+    const int  batch = 4;
+    const uint64_t elems = (uint64_t)batch*N*N;
+
+    // Random diagonally-dominant batch: A = N*I + R, |R_ij| <= 1
+    std::vector<ComplexF> Ahost(elems);
+    for(uint64_t i=0;i<elems;i++) Ahost[i] = ComplexF(dist(rng),dist(rng));
+    for(int b=0;b<batch;b++)
+      for(int64_t d=0;d<N;d++)
+	Ahost[(uint64_t)b*N*N + d*N + d] += ComplexF((RealF)N,0.0);
+
+    deviceVector<ComplexF> Adev(elems);   // gets inverted in place
+    deviceVector<ComplexF> Aorig(elems);  // untouched copy for the residual
+    deviceVector<ComplexF> Cdev(elems);
+    acceleratorCopyToDevice(&Ahost[0],&Adev[0], elems*sizeof(ComplexF));
+    acceleratorCopyToDevice(&Ahost[0],&Aorig[0],elems*sizeof(ComplexF));
+
+    deviceVector<ComplexF*> Ap(batch);
+    deviceVector<ComplexF*> Op(batch);
+    deviceVector<ComplexF*> Cp(batch);
+    std::vector<ComplexF*> ptr_h(batch);
+    for(int b=0;b<batch;b++) ptr_h[b] = &Adev[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Ap[0],batch*sizeof(ComplexF*));
+    for(int b=0;b<batch;b++) ptr_h[b] = &Aorig[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Op[0],batch*sizeof(ComplexF*));
+    for(int b=0;b<batch;b++) ptr_h[b] = &Cdev[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Cp[0],batch*sizeof(ComplexF*));
+
+    inverse.inverseBatched(N,Ap);                      // A <- A^{-1}
+
+    ComplexF one (1.0,0.0);
+    ComplexF zero(0.0,0.0);
+    blas.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
+		     (int)N,(int)N,(int)N,
+		     one, Op, Ap, zero, Cp);           // C = A A^{-1}
+    blas.synchronise();
+
+    std::vector<ComplexF> Chost(elems);
+    acceleratorCopyFromDevice(&Cdev[0],&Chost[0],elems*sizeof(ComplexF));
+
+    double maxdev = 0.0;
+    for(int b=0;b<batch;b++)
+      for(int64_t j=0;j<N;j++)
+	for(int64_t i=0;i<N;i++){
+	  ComplexF expect = (i==j) ? ComplexF(1.0,0.0) : ComplexF(0.0,0.0);
+	  ComplexF got = Chost[(uint64_t)b*N*N + j*N + i];
+	  maxdev = std::max(maxdev,(double)abs(got-expect));
+	}
+    std::cout << GridLogMessage << "T1a inverseBatched ComplexF  ||A Ainv - I||_max = "
+	      << maxdev << ( maxdev < 1.0e-4 ? "   PASS" : "   FAIL" ) << std::endl;
+    GRID_ASSERT(maxdev < 1.0e-4);
+  }
+
+  ////////////////////////////////////////////////////////////////
+  // T1b : batched inversion, ComplexD
+  ////////////////////////////////////////////////////////////////
+  {
+    const int64_t N = 48;
+    const int  batch = 3;
+    const uint64_t elems = (uint64_t)batch*N*N;
+
+    std::vector<ComplexD> Ahost(elems);
+    for(uint64_t i=0;i<elems;i++) Ahost[i] = ComplexD(dist(rng),dist(rng));
+    for(int b=0;b<batch;b++)
+      for(int64_t d=0;d<N;d++)
+	Ahost[(uint64_t)b*N*N + d*N + d] += ComplexD((RealD)N,0.0);
+
+    deviceVector<ComplexD> Adev(elems);
+    deviceVector<ComplexD> Aorig(elems);
+    deviceVector<ComplexD> Cdev(elems);
+    acceleratorCopyToDevice(&Ahost[0],&Adev[0], elems*sizeof(ComplexD));
+    acceleratorCopyToDevice(&Ahost[0],&Aorig[0],elems*sizeof(ComplexD));
+
+    deviceVector<ComplexD*> Ap(batch);
+    deviceVector<ComplexD*> Op(batch);
+    deviceVector<ComplexD*> Cp(batch);
+    std::vector<ComplexD*> ptr_h(batch);
+    for(int b=0;b<batch;b++) ptr_h[b] = &Adev[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Ap[0],batch*sizeof(ComplexD*));
+    for(int b=0;b<batch;b++) ptr_h[b] = &Aorig[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Op[0],batch*sizeof(ComplexD*));
+    for(int b=0;b<batch;b++) ptr_h[b] = &Cdev[(uint64_t)b*N*N];
+    acceleratorCopyToDevice(&ptr_h[0],&Cp[0],batch*sizeof(ComplexD*));
+
+    inverse.inverseBatched(N,Ap);
+
+    ComplexD one (1.0,0.0);
+    ComplexD zero(0.0,0.0);
+    blas.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
+		     (int)N,(int)N,(int)N,
+		     one, Op, Ap, zero, Cp);
+    blas.synchronise();
+
+    std::vector<ComplexD> Chost(elems);
+    acceleratorCopyFromDevice(&Cdev[0],&Chost[0],elems*sizeof(ComplexD));
+
+    double maxdev = 0.0;
+    for(int b=0;b<batch;b++)
+      for(int64_t j=0;j<N;j++)
+	for(int64_t i=0;i<N;i++){
+	  ComplexD expect = (i==j) ? ComplexD(1.0,0.0) : ComplexD(0.0,0.0);
+	  ComplexD got = Chost[(uint64_t)b*N*N + j*N + i];
+	  maxdev = std::max(maxdev,(double)abs(got-expect));
+	}
+    std::cout << GridLogMessage << "T1b inverseBatched ComplexD  ||A Ainv - I||_max = "
+	      << maxdev << ( maxdev < 1.0e-10 ? "   PASS" : "   FAIL" ) << std::endl;
+    GRID_ASSERT(maxdev < 1.0e-10);
+  }
+
+  ////////////////////////////////////////////////////////////////
+  // T2 : explicit-ld gemmBatched == software split-K miniature.
+  // Slab layout as in the dense coarse-coarse apply: A is K x nrows
+  // column major (lda=K); Y = op(A,T).X with X K x nrhs (ldb=K).
+  ////////////////////////////////////////////////////////////////
+  {
+    const int nrows = 8;
+    const int nrhs  = 4;
+    const int K     = 256;
+    const int NK    = 8;        // split-K chunks
+    const int Kc    = K/NK;
+
+    std::vector<ComplexF> Ahost((uint64_t)K*nrows);
+    std::vector<ComplexF> Xhost((uint64_t)K*nrhs);
+    for(auto &z : Ahost) z = ComplexF(dist(rng),dist(rng));
+    for(auto &z : Xhost) z = ComplexF(dist(rng),dist(rng));
+
+    deviceVector<ComplexF> Adev(Ahost.size());
+    deviceVector<ComplexF> Xdev(Xhost.size());
+    deviceVector<ComplexF> Yref((uint64_t)nrows*nrhs);
+    deviceVector<ComplexF> Ypart((uint64_t)NK*nrows*nrhs);
+    acceleratorCopyToDevice(&Ahost[0],&Adev[0],Ahost.size()*sizeof(ComplexF));
+    acceleratorCopyToDevice(&Xhost[0],&Xdev[0],Xhost.size()*sizeof(ComplexF));
+
+    ComplexF one (1.0,0.0);
+    ComplexF zero(0.0,0.0);
+
+    // (a) reference: one compact batch=1 call (compact lda == K for OP_T)
+    {
+      deviceVector<ComplexF*> Ap(1), Xp(1), Yp(1);
+      std::vector<ComplexF*> h(1);
+      h[0]=&Adev[0]; acceleratorCopyToDevice(&h[0],&Ap[0],sizeof(ComplexF*));
+      h[0]=&Xdev[0]; acceleratorCopyToDevice(&h[0],&Xp[0],sizeof(ComplexF*));
+      h[0]=&Yref[0]; acceleratorCopyToDevice(&h[0],&Yp[0],sizeof(ComplexF*));
+      blas.gemmBatched(GridBLAS_OP_T,GridBLAS_OP_N,
+		       nrows,nrhs,K,
+		       one, Ap, Xp, zero, Yp);
+      blas.synchronise();
+    }
+
+    // (b) split-K: NK chunk-pointers into the SAME allocations, lda/ldb = K
+    {
+      deviceVector<ComplexF*> Ap(NK), Xp(NK), Yp(NK);
+      std::vector<ComplexF*> h(NK);
+      for(int j=0;j<NK;j++) h[j] = &Adev[(uint64_t)j*Kc];        // K-offset slice
+      acceleratorCopyToDevice(&h[0],&Ap[0],NK*sizeof(ComplexF*));
+      for(int j=0;j<NK;j++) h[j] = &Xdev[(uint64_t)j*Kc];
+      acceleratorCopyToDevice(&h[0],&Xp[0],NK*sizeof(ComplexF*));
+      for(int j=0;j<NK;j++) h[j] = &Ypart[(uint64_t)j*nrows*nrhs];
+      acceleratorCopyToDevice(&h[0],&Yp[0],NK*sizeof(ComplexF*));
+      blas.gemmBatched(GridBLAS_OP_T,GridBLAS_OP_N,
+		       nrows,nrhs,Kc,
+		       one, Ap, /*lda*/ K,
+		            Xp, /*ldb*/ K,
+		       zero,Yp, /*ldc*/ nrows);
+      blas.synchronise();
+    }
+
+    std::vector<ComplexF> Yref_h((uint64_t)nrows*nrhs);
+    std::vector<ComplexF> Ypart_h((uint64_t)NK*nrows*nrhs);
+    acceleratorCopyFromDevice(&Yref[0], &Yref_h[0], Yref_h.size()*sizeof(ComplexF));
+    acceleratorCopyFromDevice(&Ypart[0],&Ypart_h[0],Ypart_h.size()*sizeof(ComplexF));
+
+    double maxdev = 0.0;
+    double maxval = 0.0;
+    for(int i=0;i<nrows*nrhs;i++){
+      ComplexF sum(0.0,0.0);
+      for(int j=0;j<NK;j++) sum = sum + Ypart_h[(uint64_t)j*nrows*nrhs + i];
+      maxdev = std::max(maxdev,(double)abs(sum-Yref_h[i]));
+      maxval = std::max(maxval,(double)abs(Yref_h[i]));
+    }
+    double rel = maxdev/maxval;
+    std::cout << GridLogMessage << "T2  split-K strided gemmBatched  max rel dev vs compact = "
+	      << rel << ( rel < 1.0e-4 ? "   PASS" : "   FAIL" ) << std::endl;
+    GRID_ASSERT(rel < 1.0e-4);
+  }
+
+  std::cout << GridLogMessage << "All batched-blas tests PASSED" << std::endl;
+
+  Grid_finalize();
+}

tests/forces/Test_dwf_ratio_4dpf_force.cc

@@ -0,0 +1,279 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/forces/Test_dwf_ratio_4dpf_force.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 */
+
+//
+// Validation of TwoFlavourRatio4DPseudoFermionAction (non-EO, LinearFunction
+// solver slots).  Three tests, run for BOTH wall conventions:
+//
+//  T1 HeatbathIdentityTest : refresh then S; PASS iff S == 0.5*|eta4|^2 to
+//     solver tolerance.  This adjudicates the 4D effective-operator
+//     composition identity [P M^-1 V Pdag][P V^-1 M Pdag] = 1 for the chosen
+//     (P,Pdag) wall pair.  NO PREDICTION is made about which convention
+//     passes -- that is what the test decides.
+//  T2 ForceTest (idiom from Test_double_ratio.cc) : midpoint-derivative
+//     check of deriv against S.  Should PASS for BOTH conventions (S and
+//     deriv use the same literal-adjoint pair by construction).
+//  T3 Trivial-ratio control (V == M) : T1 with NumOp = DenOp.  The solve
+//     cancels against the multiply, so S = 0.5|eta4|^2 requires only
+//     P Pdag = 1_4d.  Should PASS for BOTH conventions; isolates plumbing
+//     from the composition identity.
+//
+// Solvers here are plain CG on the normal equations (CGNR), tolerance 1e-12,
+// so every defect above ~1e-10 is structural, not solver noise.  Run small,
+// e.g.:  ./Test_dwf_ratio_4dpf_force --grid 8.8.8.8
+//
+#include <Grid/Grid.h>
+#include <Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h>
+
+using namespace std;
+using namespace Grid;
+
+////////////////////////////////////////////////////////////////////
+// LinearFunction wrappers: direct M^-1 and M^-dag via CG on the
+// normal equations.  These stand in for the MG-GCR stack in this
+// test; the action class sees only LinearFunction.
+////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class CGNRLinearFunction : public LinearFunction<Field> {   // out = M^-1 in
+public:
+  using LinearFunction<Field>::operator();
+  Matrix &_Mat; RealD tol; Integer maxit;
+  CGNRLinearFunction(Matrix &Mat,RealD _tol,Integer _maxit) : _Mat(Mat), tol(_tol), maxit(_maxit) {};
+  void operator()(const Field &in, Field &out) {
+    MdagMLinearOperator<Matrix,Field> MdagM(_Mat);
+    Field src(in.Grid());
+    _Mat.Mdag(in,src);                 // src = Mdag in
+    ConjugateGradient<Field> CG(tol,maxit);
+    out = Zero();
+    CG(MdagM,src,out);                 // out = (MdagM)^-1 Mdag in = M^-1 in
+  }
+};
+template<class Matrix,class Field>
+class CGNRDagLinearFunction : public LinearFunction<Field> { // out = M^-dag in
+public:
+  using LinearFunction<Field>::operator();
+  Matrix &_Mat; RealD tol; Integer maxit;
+  CGNRDagLinearFunction(Matrix &Mat,RealD _tol,Integer _maxit) : _Mat(Mat), tol(_tol), maxit(_maxit) {};
+  void operator()(const Field &in, Field &out) {
+    MdagMLinearOperator<Matrix,Field> MdagM(_Mat);
+    Field tmp(in.Grid());
+    tmp = Zero();
+    ConjugateGradient<Field> CG(tol,maxit);
+    CG(MdagM,in,tmp);                  // tmp = (MdagM)^-1 in
+    _Mat.M(tmp,out);                   // out = M (MdagM)^-1 in = M^-dag in
+  }
+};
+
+////////////////////////////////////////////////////////////////////
+// T1 / T3 : heatbath composition-identity test.
+// Twin-seeded RNG reproduces the eta4 drawn inside refresh.
+////////////////////////////////////////////////////////////////////
+template<class Impl>
+RealD HeatbathIdentityTest(TwoFlavourRatio4DPseudoFermionAction<Impl> &action,
+			   LatticeGaugeField &U,
+			   GridCartesian *UGrid,
+			   const std::string &tag)
+{
+  typedef typename Impl::FermionField FermionField;
+
+  std::vector<int> seeds({9,11,13,17});
+  GridSerialRNG   sRNG;                 sRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG RNG4(UGrid);          RNG4.SeedFixedIntegers(seeds);
+  GridParallelRNG RNG4check(UGrid);     RNG4check.SeedFixedIntegers(seeds);
+
+  FermionField eta4check(UGrid);
+  gaussian(RNG4check,eta4check);        // identical to the draw inside refresh
+
+  action.refresh(U,sRNG,RNG4);
+  RealD S       = action.S(U);
+  RealD Sexpect = 0.5*norm2(eta4check);
+  RealD defect  = std::abs(S-Sexpect)/Sexpect;
+
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+  std::cout << GridLogMessage << " HeatbathIdentityTest ["<<tag<<"]" << std::endl;
+  std::cout << GridLogMessage << "   S              = " << S << std::endl;
+  std::cout << GridLogMessage << "   0.5|eta4|^2    = " << Sexpect << std::endl;
+  std::cout << GridLogMessage << "   relative defect = " << defect
+	    << ( defect < 1.0e-8 ? "   PASS" : "   FAIL" ) << std::endl;
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+  return defect;
+}
+
+////////////////////////////////////////////////////////////////////
+// T2 : ForceTest idiom from Test_double_ratio.cc (midpoint derivative)
+////////////////////////////////////////////////////////////////////
+template<class Gimpl>
+void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
+{
+  GridBase *UGrid = U.Grid();
+
+  std::vector<int> seeds({1,2,3,5});
+  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
+
+  LatticeColourMatrix Pmu(UGrid);
+  LatticeGaugeField P(UGrid);
+  LatticeGaugeField UdSdU(UGrid);
+
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+
+  RealD eps=0.005;
+
+  Gimpl::generate_momenta(P,sRNG,RNG4);
+  Filter.applyFilter(P);
+
+  action.refresh(U,sRNG,RNG4);
+
+  RealD S1 = action.S(U);
+
+  Gimpl::update_field(P,U,eps);
+
+  action.deriv(U,UdSdU);
+  UdSdU = Ta(UdSdU);
+  Filter.applyFilter(UdSdU);
+
+  DumpSliceNorm("Force",UdSdU,Nd-1);
+
+  Gimpl::update_field(P,U,eps);
+
+  RealD S2 = action.S(U);
+
+  // Use the derivative
+  LatticeComplex dS(UGrid); dS = Zero();
+  for(int mu=0;mu<Nd;mu++){
+    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
+    Pmu= PeekIndex<LorentzIndex>(P,mu);
+    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
+  }
+  ComplexD dSpred    = sum(dS);
+  RealD diff =  S2-S1-dSpred.real();
+
+  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
+  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
+  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
+  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
+  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
+  std::cout<< GridLogMessage << "diff/dS : "<< diff/(S2-S1)<<std::endl;
+  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
+  //  GRID_ASSERT(diff<1.0);
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::cout << std::setprecision(14);
+
+  const int Ls=8;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  GridParallelRNG RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField U(UGrid);
+  SU<Nc>::HotConfiguration(RNG4,U);
+
+  ////////////////////////////////////////////////////////////////
+  // Operators: Mobius, campaign-like b,c; heavyish masses so CGNR
+  // is fast and well-conditioned even on a hot configuration.
+  ////////////////////////////////////////////////////////////////
+  RealD mden = 0.2;
+  RealD mnum = 0.5;
+  RealD M5   = 1.8;
+  RealD b    = 1.5;
+  RealD c    = 0.5;
+
+  WilsonImplParams p;
+  p.boundary_phases[0] = 1.0;
+  p.boundary_phases[1] = 1.0;
+  p.boundary_phases[2] = 1.0;
+  p.boundary_phases[3] = -1.0;
+
+  MobiusFermionD DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mden,M5,b,c,p);
+  MobiusFermionD NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mnum,M5,b,c,p);
+
+  RealD tol   = 1.0e-12;
+  Integer maxit = 20000;
+
+  typedef WilsonImplD::FermionField FermionField;
+  CGNRLinearFunction<MobiusFermionD,FermionField>    MinvSolver   (DenOp,tol,maxit);
+  CGNRDagLinearFunction<MobiusFermionD,FermionField> MdagInvSolver(DenOp,tol,maxit);
+  CGNRLinearFunction<MobiusFermionD,FermionField>    VinvSolver   (NumOp,tol,maxit);
+
+  ////////////////////////////////////////////////////////////////
+  // Actions: both wall conventions, plus the V==M trivial control
+  ////////////////////////////////////////////////////////////////
+  TwoFlavourRatio4DPseudoFermionAction<WilsonImplD> ActSol(NumOp,DenOp,MinvSolver,MdagInvSolver,MinvSolver,VinvSolver,1);
+  TwoFlavourRatio4DPseudoFermionAction<WilsonImplD> ActSrc(NumOp,DenOp,MinvSolver,MdagInvSolver,MinvSolver,VinvSolver,0);
+  TwoFlavourRatio4DPseudoFermionAction<WilsonImplD> ActTrivSol(DenOp,DenOp,MinvSolver,MdagInvSolver,MinvSolver,MinvSolver,1);
+  TwoFlavourRatio4DPseudoFermionAction<WilsonImplD> ActTrivSrc(DenOp,DenOp,MinvSolver,MdagInvSolver,MinvSolver,MinvSolver,0);
+
+  ////////////////////////////////////////////////////////////////
+  // T3 controls first (must both pass; isolates plumbing)
+  ////////////////////////////////////////////////////////////////
+  RealD d3s = HeatbathIdentityTest(ActTrivSol,U,UGrid,"T3 trivial V==M, solution walls");
+  RealD d3q = HeatbathIdentityTest(ActTrivSrc,U,UGrid,"T3 trivial V==M, source walls");
+
+  ////////////////////////////////////////////////////////////////
+  // T1 : the composition-identity adjudication
+  ////////////////////////////////////////////////////////////////
+  RealD d1s = HeatbathIdentityTest(ActSol,U,UGrid,"T1 ratio, solution walls");
+  RealD d1q = HeatbathIdentityTest(ActSrc,U,UGrid,"T1 ratio, source walls");
+
+  ////////////////////////////////////////////////////////////////
+  // T2 : force consistency (expected PASS for both conventions)
+  ////////////////////////////////////////////////////////////////
+  MomentumFilterNone<LatticeGaugeField> FilterNone;
+  ForceTest<GimplTypesR>(ActSol,U,FilterNone);
+  ForceTest<GimplTypesR>(ActSrc,U,FilterNone);
+
+  ////////////////////////////////////////////////////////////////
+  // Summary
+  ////////////////////////////////////////////////////////////////
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+  std::cout << GridLogMessage << " SUMMARY (relative heatbath defects)" << std::endl;
+  std::cout << GridLogMessage << "   T3 trivial  solution walls : " << d3s << std::endl;
+  std::cout << GridLogMessage << "   T3 trivial  source   walls : " << d3q << std::endl;
+  std::cout << GridLogMessage << "   T1 ratio    solution walls : " << d1s << std::endl;
+  std::cout << GridLogMessage << "   T1 ratio    source   walls : " << d1q << std::endl;
+  std::cout << GridLogMessage << " T3 must pass for both; T1 selects the wall convention." << std::endl;
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+
+  GRID_ASSERT(d3s < 1.0e-8);
+  GRID_ASSERT(d3q < 1.0e-8);
+
+  Grid_finalize();
+}

tests/forces/Test_dwf_ratio_leftprec.cc

@@ -0,0 +1,341 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/forces/Test_dwf_ratio_leftprec.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 */
+
+//
+// Correctness of TwoFlavourRatioLeftPrecPseudoFermionAction against the
+// decades-proven TwoFlavourRatioPseudoFermionAction.  Both classes compute
+// the SAME action S = phi^dag V (MdagM)^-1 Vdag phi through different solve
+// chains (normal-equations vs left-preconditioned F = Vdag M), so with
+// twin-seeded refreshes and 1e-12 solvers they must agree to solver
+// tolerance.  Tests:
+//
+//   E0a/E0b : heatbath identity, S == 0.5|eta|^2 after RNG refresh, for
+//             BOTH classes (E0a also validates the twin-eta capture).
+//   E1      : S_classic == S_leftprec       (relative, ~1e-8)
+//   E2      : deriv_classic == deriv_leftprec (pointwise field norm, ~1e-8)
+//   F1      : ForceTest (Test_double_ratio.cc idiom) on the LeftPrec class.
+//
+// All asserts are hard: this is the regression gate for the new class.
+// Run small, e.g.:  ./Test_dwf_ratio_leftprec --grid 8.8.8.8
+//
+#include <Grid/Grid.h>
+#include <Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h>
+
+using namespace std;
+using namespace Grid;
+
+////////////////////////////////////////////////////////////////////
+// Minimal LinearOperator for the composite F = Vdag M, exposing the
+// Hermitian normal operator FdagF for CG.  Stencil entries assert.
+////////////////////////////////////////////////////////////////////
+template<class Impl>
+class VdagMNormalOperator : public LinearOperatorBase<typename Impl::FermionField> {
+public:
+  typedef typename Impl::FermionField Field;
+  FermionOperator<Impl> &VOp;
+  FermionOperator<Impl> &MOp;
+  VdagMNormalOperator(FermionOperator<Impl> &V,FermionOperator<Impl> &M) : VOp(V), MOp(M) {};
+
+  void Fapply(const Field &in, Field &out) {          // out = Vdag M in
+    Field tmp(in.Grid());
+    MOp.M(in,tmp);
+    VOp.Mdag(tmp,out);
+  }
+  void FdagApply(const Field &in, Field &out) {       // out = Mdag V in
+    Field tmp(in.Grid());
+    VOp.M(in,tmp);
+    MOp.Mdag(tmp,out);
+  }
+  virtual void Op     (const Field &in, Field &out) { Fapply(in,out); }
+  virtual void AdjOp  (const Field &in, Field &out) { FdagApply(in,out); }
+  virtual void HermOp (const Field &in, Field &out) {
+    Field tmp(in.Grid());
+    Fapply(in,tmp);
+    FdagApply(tmp,out);
+  }
+  virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) {
+    HermOp(in,out);
+    ComplexD dot = innerProduct(in,out);
+    n1 = real(dot);
+    n2 = norm2(out);
+  }
+  virtual void OpDiag (const Field &in, Field &out)                  { GRID_ASSERT(0); }
+  virtual void OpDir  (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); }
+  virtual void OpDirAll(const Field &in, std::vector<Field> &out)    { GRID_ASSERT(0); }
+};
+
+////////////////////////////////////////////////////////////////////
+// F-contract LinearFunctions for the test, both via CG on FdagF:
+//   forward :  F x = b     ==>  x = (FdagF)^-1 Fdag b
+//   adjoint :  Fdag z = b  ==>  z = F (FdagF)^-1 b
+////////////////////////////////////////////////////////////////////
+template<class Impl>
+class ForwardFSolve : public LinearFunction<typename Impl::FermionField> {
+public:
+  typedef typename Impl::FermionField Field;
+  using LinearFunction<Field>::operator();
+  VdagMNormalOperator<Impl> &FdagF; RealD tol; Integer maxit;
+  ForwardFSolve(VdagMNormalOperator<Impl> &Op,RealD _tol,Integer _maxit) : FdagF(Op), tol(_tol), maxit(_maxit) {};
+  void operator()(const Field &in, Field &out) {
+    Field src(in.Grid());
+    FdagF.FdagApply(in,src);
+    ConjugateGradient<Field> CG(tol,maxit);
+    out = Zero();
+    CG(FdagF,src,out);
+  }
+};
+template<class Impl>
+class AdjointFSolve : public LinearFunction<typename Impl::FermionField> {
+public:
+  typedef typename Impl::FermionField Field;
+  using LinearFunction<Field>::operator();
+  VdagMNormalOperator<Impl> &FdagF; RealD tol; Integer maxit;
+  AdjointFSolve(VdagMNormalOperator<Impl> &Op,RealD _tol,Integer _maxit) : FdagF(Op), tol(_tol), maxit(_maxit) {};
+  void operator()(const Field &in, Field &out) {
+    Field y(in.Grid());
+    y = Zero();
+    ConjugateGradient<Field> CG(tol,maxit);
+    CG(FdagF,in,y);
+    FdagF.Fapply(y,out);
+  }
+};
+template<class Matrix,class Field>
+class NormalEqSolve : public LinearFunction<Field> {   // out = (MdagM)^-1 in
+public:
+  using LinearFunction<Field>::operator();
+  Matrix &_Mat; RealD tol; Integer maxit;
+  NormalEqSolve(Matrix &Mat,RealD _tol,Integer _maxit) : _Mat(Mat), tol(_tol), maxit(_maxit) {};
+  void operator()(const Field &in, Field &out) {
+    MdagMLinearOperator<Matrix,Field> MdagM(_Mat);
+    ConjugateGradient<Field> CG(tol,maxit);
+    out = Zero();
+    CG(MdagM,in,out);
+  }
+};
+
+////////////////////////////////////////////////////////////////////
+// ForceTest idiom from Test_double_ratio.cc (midpoint derivative)
+////////////////////////////////////////////////////////////////////
+template<class Gimpl>
+void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
+{
+  GridBase *UGrid = U.Grid();
+
+  std::vector<int> seeds({1,2,3,5});
+  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
+
+  LatticeColourMatrix Pmu(UGrid);
+  LatticeGaugeField P(UGrid);
+  LatticeGaugeField UdSdU(UGrid);
+
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+
+  RealD eps=0.005;
+
+  Gimpl::generate_momenta(P,sRNG,RNG4);
+  Filter.applyFilter(P);
+
+  action.refresh(U,sRNG,RNG4);
+
+  RealD S1 = action.S(U);
+
+  Gimpl::update_field(P,U,eps);
+
+  action.deriv(U,UdSdU);
+  UdSdU = Ta(UdSdU);
+  Filter.applyFilter(UdSdU);
+
+  DumpSliceNorm("Force",UdSdU,Nd-1);
+
+  Gimpl::update_field(P,U,eps);
+
+  RealD S2 = action.S(U);
+
+  LatticeComplex dS(UGrid); dS = Zero();
+  for(int mu=0;mu<Nd;mu++){
+    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
+    Pmu= PeekIndex<LorentzIndex>(P,mu);
+    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
+  }
+  ComplexD dSpred    = sum(dS);
+  RealD diff =  S2-S1-dSpred.real();
+
+  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
+  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
+  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
+  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
+  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
+  std::cout<< GridLogMessage << "diff/dS : "<< diff/(S2-S1)<<std::endl;
+  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::cout << std::setprecision(14);
+
+  const int Ls=8;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  GridParallelRNG RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField U(UGrid);
+  SU<Nc>::HotConfiguration(RNG4,U);
+
+  ////////////////////////////////////////////////////////////////
+  // Operators: quotient pair (V = PV mass 1, M light-ish), Mobius,
+  // campaign b,c.  Heavyish M so CG is quick on a hot configuration.
+  ////////////////////////////////////////////////////////////////
+  RealD mden = 0.2;
+  RealD mnum = 1.0;
+  RealD M5   = 1.8;
+  RealD b    = 1.5;
+  RealD c    = 0.5;
+
+  WilsonImplParams p;
+  p.boundary_phases[0] = 1.0;
+  p.boundary_phases[1] = 1.0;
+  p.boundary_phases[2] = 1.0;
+  p.boundary_phases[3] = -1.0;
+
+  MobiusFermionD DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mden,M5,b,c,p);
+  MobiusFermionD NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mnum,M5,b,c,p);
+
+  RealD tol     = 1.0e-12;
+  Integer maxit = 30000;
+
+  typedef WilsonImplD::FermionField FermionField;
+
+  ////////////////////////////////////////////////////////////////
+  // Solvers.  Classic: CG as OperatorFunction on the MdagM linop the
+  // action supplies.  LeftPrec: F-contract solves via CG on FdagF.
+  ////////////////////////////////////////////////////////////////
+  ConjugateGradient<FermionField> CG(tol,maxit);
+
+  VdagMNormalOperator<WilsonImplD> FdagF(NumOp,DenOp);
+  ForwardFSolve<WilsonImplD>       Ffwd (FdagF,tol,maxit);
+  AdjointFSolve<WilsonImplD>       Fadj (FdagF,tol,maxit);
+  NormalEqSolve<MobiusFermionD,FermionField> VdagVinv(NumOp,tol,maxit);
+
+  TwoFlavourRatioPseudoFermionAction<WilsonImplD>         Classic (NumOp,DenOp,CG,CG);
+  TwoFlavourRatioLeftPrecPseudoFermionAction<WilsonImplD> LeftPrec(NumOp,DenOp,Ffwd,Fadj,Fadj,VdagVinv);
+
+  ////////////////////////////////////////////////////////////////
+  // Twin-seeded refreshes: identical eta into both classes.
+  ////////////////////////////////////////////////////////////////
+  std::vector<int> seedsR({9,11,13,17});
+  GridSerialRNG   sRNGa;         sRNGa.SeedFixedIntegers(seedsR);
+  GridSerialRNG   sRNGb;         sRNGb.SeedFixedIntegers(seedsR);
+  GridParallelRNG RNG5a(FGrid);  RNG5a.SeedFixedIntegers(seedsR);
+  GridParallelRNG RNG5b(FGrid);  RNG5b.SeedFixedIntegers(seedsR);
+  GridParallelRNG RNG5c(FGrid);  RNG5c.SeedFixedIntegers(seedsR);
+
+  FermionField etaTwin(FGrid);
+  gaussian(RNG5c,etaTwin);                    // identical to both refresh draws
+
+  Classic.refresh (U,sRNGa,RNG5a);
+  LeftPrec.refresh(U,sRNGb,RNG5b);
+
+  ////////////////////////////////////////////////////////////////
+  // E0 : heatbath identity for both classes
+  ////////////////////////////////////////////////////////////////
+  RealD Sexpect = 0.5*norm2(etaTwin);
+  RealD Sc = Classic.S(U);
+  RealD Sl = LeftPrec.S(U);
+
+  RealD e0a = std::abs(Sc-Sexpect)/Sexpect;
+  RealD e0b = std::abs(Sl-Sexpect)/Sexpect;
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+  std::cout << GridLogMessage << " E0 heatbath identity: 0.5|eta|^2 = " << Sexpect << std::endl;
+  std::cout << GridLogMessage << "   classic  S = " << Sc << "  rel defect " << e0a
+	    << ( e0a < 1.0e-8 ? "   PASS" : "   FAIL" ) << std::endl;
+  std::cout << GridLogMessage << "   leftprec S = " << Sl << "  rel defect " << e0b
+	    << ( e0b < 1.0e-8 ? "   PASS" : "   FAIL" ) << std::endl;
+
+  ////////////////////////////////////////////////////////////////
+  // E1 : action equivalence
+  ////////////////////////////////////////////////////////////////
+  RealD e1 = std::abs(Sc-Sl)/std::abs(Sc);
+  std::cout << GridLogMessage << " E1 action equivalence: rel diff = " << e1
+	    << ( e1 < 1.0e-8 ? "   PASS" : "   FAIL" ) << std::endl;
+
+  ////////////////////////////////////////////////////////////////
+  // E2 : derivative equivalence (pointwise field comparison)
+  ////////////////////////////////////////////////////////////////
+  LatticeGaugeField dSdUc(UGrid);
+  LatticeGaugeField dSdUl(UGrid);
+  LatticeGaugeField dDiff(UGrid);
+
+  Classic.deriv (U,dSdUc);
+  LeftPrec.deriv(U,dSdUl);
+  dDiff = dSdUc - dSdUl;
+
+  RealD e2 = std::sqrt( norm2(dDiff) / norm2(dSdUc) );
+  std::cout << GridLogMessage << " E2 deriv equivalence: |diff|/|classic| = " << e2
+	    << ( e2 < 1.0e-8 ? "   PASS" : "   FAIL" ) << std::endl;
+  std::cout << GridLogMessage << "   |dSdU classic |^2 = " << norm2(dSdUc) << std::endl;
+  std::cout << GridLogMessage << "   |dSdU leftprec|^2 = " << norm2(dSdUl) << std::endl;
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+
+  ////////////////////////////////////////////////////////////////
+  // F1 : standalone force test on the LeftPrec class
+  ////////////////////////////////////////////////////////////////
+  MomentumFilterNone<LatticeGaugeField> FilterNone;
+  ForceTest<GimplTypesR>(LeftPrec,U,FilterNone);
+
+  ////////////////////////////////////////////////////////////////
+  // Summary + hard asserts (this is the regression gate)
+  ////////////////////////////////////////////////////////////////
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+  std::cout << GridLogMessage << " SUMMARY" << std::endl;
+  std::cout << GridLogMessage << "   E0a classic  heatbath defect : " << e0a << std::endl;
+  std::cout << GridLogMessage << "   E0b leftprec heatbath defect : " << e0b << std::endl;
+  std::cout << GridLogMessage << "   E1  action   equivalence     : " << e1  << std::endl;
+  std::cout << GridLogMessage << "   E2  deriv    equivalence     : " << e2  << std::endl;
+  std::cout << GridLogMessage << "=========================================================" << std::endl;
+
+  GRID_ASSERT(e0a < 1.0e-8);
+  GRID_ASSERT(e0b < 1.0e-8);
+  GRID_ASSERT(e1  < 1.0e-8);
+  GRID_ASSERT(e2  < 1.0e-8);
+
+  std::cout << GridLogMessage << "All equivalence tests PASSED" << std::endl;
+
+  Grid_finalize();
+}