Updates

tests/debug/Test_blockcyclic.cc

@@ -20,7 +20,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 //////////////////////////////////////////////////////////////////////////////
 // Regression gate for BlockCyclicLayout -- stage 1 of the 2D distributed
-// dense inverse (documentation/DistributedDenseInverse2D.tex).
+// dense inverse.
 //
 // The layout is pure index arithmetic, so this test is EXHAUSTIVE rather
 // than statistical: every stage sweeps a battery of (N, nb, Pr, Pc)

tests/debug/Test_coarse.cc

@@ -0,0 +1,261 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_coarse.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 */
+
+//
+// MultiGeneralCoarsenedOperator against the existing mrhs coarse operator.
+//
+// The reference (the deprecated operator) is constructed on the D+1 grid; the
+// operator under test on the D dimensional grid, with SetGrid() adopting the
+// caller owned D+1 grid and building its padded cell and neighbour table from
+// the D dimensional stencil, with the Nrhs factor multiplied in.
+//
+// Both are given identical matrix elements, so any difference in the apply is
+// the restructured neighbour table. The same geometry object is passed to
+// both: the reference adds one to skip for the rhs direction, the operator
+// under test uses it as is on the D dimensional grid, so both describe the
+// same stencil over the D dimensions.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const int nbasis = 8;
+
+typedef vSpinColourVector FineObj;
+typedef sTComplexD        CComplexT;   // unvectorised coarse space
+
+typedef DeprecatedMultiGeneralCoarsenedMatrix    <FineObj,CComplexT,nbasis> RefOperator;
+typedef MultiGeneralCoarsenedOperator<FineObj,CComplexT,nbasis> TestOperator;
+
+////////////////////////////////////////////////////////////////////////
+// Identical random matrix elements into both operators
+////////////////////////////////////////////////////////////////////////
+template<class OpA,class OpB>
+void SeedMatrixElements(OpA &A,OpB &B,int npoint,GridSerialRNG &sRNG)
+{
+  typedef typename OpA::calcMatrix calcMatrix;
+
+  deviceVector<calcMatrix> bufA,bufB;
+  for(int p=0;p<npoint;p++){
+
+    // MatrixPointOut/In carry one stencil point in lSite order whichever
+    // internal layout the operator uses.
+    A.MatrixPointOut(p,bufA);
+    B.MatrixPointOut(p,bufB);
+    GRID_ASSERT(bufA.size() == bufB.size());
+
+    int64_t sites = bufA.size();
+    std::vector<calcMatrix> host(sites);
+
+    ComplexD *w = (ComplexD *)&host[0];
+    int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      RealD re,im;
+      random(sRNG,re);
+      random(sRNG,im);
+      w[i] = ComplexD(re-0.5,im-0.5);
+    }
+
+    acceleratorCopyToDevice(&host[0],&bufA[0],sites*sizeof(calcMatrix));
+    A.MatrixPointIn(p,bufA);
+    B.MatrixPointIn(p,bufA);
+  }
+}
+
+template<class Op>
+RealD MatrixChecksum(Op &O,int npoint)
+{
+  typedef typename Op::calcMatrix calcMatrix;
+  RealD sum=0.0;
+  deviceVector<calcMatrix> buf;
+  for(int p=0;p<npoint;p++){
+    O.MatrixPointOut(p,buf);
+    int64_t sites = buf.size();
+    std::vector<calcMatrix> host(sites);
+    acceleratorCopyFromDevice(&buf[0],&host[0],sites*sizeof(calcMatrix));
+    ComplexD *w = (ComplexD *)&host[0];
+    int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++) sum += real(w[i])*real(w[i]) + imag(w[i])*imag(w[i]);
+  }
+  return sum;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int nrhs = 4;
+
+  Coordinate clatt = GridDefaultLatt();
+  Coordinate csimd(Nd,1);                 // the coarse space is unvectorised
+  Coordinate cmpi  = GridDefaultMpi();
+
+  ////////////////////////////////////////////////
+  // D dimensional coarse grid, and D+1 for the reference
+  ////////////////////////////////////////////////
+  GridCartesian *CoarseD = new GridCartesian(clatt,csimd,cmpi);
+
+  std::cout << GridLogMessage << "coarse D   grid "; for(int d=0;d<Nd;d++) std::cout<<clatt[d]<<" ";
+  std::cout << "  Nsimd " << CoarseD->Nsimd() << std::endl;
+
+  ////////////////////////////////////////////////
+  // One geometry object for both, and one D+1 grid
+  // owned here and shared by both operators: fields
+  // conform only across a shared grid object.
+  ////////////////////////////////////////////////
+  NextToNearestStencilGeometry4D geom(CoarseD);
+
+  Coordinate mlatt(1,nrhs), msimd(1,1), mmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    mlatt.push_back(clatt[d]);
+    msimd.push_back(csimd[d]);
+    mmpi .push_back(cmpi[d]);
+  }
+  GridCartesian *CoarseMulti = new GridCartesian(mlatt,msimd,mmpi);
+
+  TestOperator OpTest(geom,CoarseD);
+  RefOperator OpRef(geom,CoarseMulti);
+  OpTest.SetGrid(CoarseMulti);
+
+  std::cout << GridLogMessage << "coarse D+1 grid nrhs " << nrhs
+	    << "  Nsimd " << CoarseMulti->Nsimd() << std::endl;
+
+  std::cout << GridLogMessage << "npoint ref " << OpRef.geom.npoint
+	    << "   npoint test " << OpTest.geom.npoint << std::endl;
+  GRID_ASSERT(OpRef.geom.npoint == OpTest.geom.npoint);
+
+  int npoint = OpRef.geom.npoint;
+
+  ////////////////////////////////////////////////
+  // Identical matrix elements
+  ////////////////////////////////////////////////
+  GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector<int>({7,8,9,10}));
+  SeedMatrixElements(OpRef,OpTest,npoint,sRNG);
+
+  RealD ckRef = MatrixChecksum(OpRef,npoint);
+  RealD ckTest = MatrixChecksum(OpTest,npoint);
+  std::cout << GridLogMessage << "matrix element checksum ref " << ckRef
+	    << "  test " << ckTest << std::endl;
+  GRID_ASSERT( ckRef == ckTest );
+
+  ////////////////////////////////////////////////
+  // Same input, compare the applies
+  ////////////////////////////////////////////////
+  typedef RefOperator::CoarseVector CoarseVector;
+  // RNG on the D dimensional grid fills any D+1 field: the rhs direction is
+  // undistributed and divides cleanly. One RNG serves every Nrhs.
+  GridParallelRNG pRNG(CoarseD); pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  CoarseVector in  (CoarseMulti); random(pRNG,in);
+  CoarseVector out1(CoarseMulti);
+  CoarseVector out2(CoarseMulti);
+  CoarseVector err (CoarseMulti);
+
+  OpRef.M(in,out1);
+  OpTest.M(in,out2);
+
+  err = out1 - out2;
+  std::cout << GridLogMessage << "|ref out|^2 = " << norm2(out1)
+	    << "   |test out|^2 = " << norm2(out2) << std::endl;
+  std::cout << GridLogMessage << "|ref - test|^2 = " << norm2(err) << std::endl;
+  GRID_ASSERT( norm2(out1) > 0.0 );
+  GRID_ASSERT( norm2(err) == 0.0 );
+
+  ////////////////////////////////////////////////
+  // SetGrid is idempotent on pointer identity
+  ////////////////////////////////////////////////
+  OpTest.SetGrid(CoarseMulti);
+  GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );
+  OpTest.M(in,out2);
+  err = out1 - out2;
+  GRID_ASSERT( norm2(err) == 0.0 );
+  std::cout << GridLogMessage << "SetGrid idempotent on identity" << std::endl;
+
+  ////////////////////////////////////////////////
+  // Move to a different Nrhs and back. The matrix
+  // elements are Nrhs independent and must survive
+  // both the release and the rebuild.
+  ////////////////////////////////////////////////
+  // Nrhs 1 is the single RHS case through the multiRHS path, and each slice
+  // of the Nrhs 4 apply must come back unchanged.
+  OpTest.M(in,out2);
+  for(int nr=2;nr>=1;nr--){
+    Coordinate latt2(1,nr), simd2(1,1), mpi2(1,1);
+    for(int d=0;d<Nd;d++){
+      latt2.push_back(clatt[d]);
+      simd2.push_back(csimd[d]);
+      mpi2 .push_back(cmpi[d]);
+    }
+    GridCartesian *CoarseMulti2 = new GridCartesian(latt2,simd2,mpi2);
+
+    OpTest.SetGrid(CoarseMulti2);
+    GRID_ASSERT( OpTest.Nrhs() == nr );
+    GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );
+
+    CoarseVector in2 (CoarseMulti2);
+    CoarseVector out(CoarseMulti2);
+    for(int r=0;r<nr;r++){
+      CoarseVector slice(CoarseD);
+      ExtractSliceFast(slice,in,r,0);
+      InsertSliceFast(slice,in2,r,0);
+    }
+    OpTest.M(in2,out);
+
+    RealD sdiff=0.0;
+    for(int r=0;r<nr;r++){
+      CoarseVector a(CoarseD),b(CoarseD),e(CoarseD);
+      ExtractSliceFast(a,out ,r,0);
+      ExtractSliceFast(b,out2,r,0);
+      e = a-b;
+      sdiff += norm2(e);
+    }
+    // Not bit exact: a different Nrhs is a different GEMM shape
+    std::cout << GridLogMessage << "Nrhs " << nr << " slices agree with Nrhs "
+	      << nrhs << " : |diff|^2/|out|^2 = " << sdiff/norm2(out) << std::endl;
+    GRID_ASSERT( norm2(out) > 0.0 );
+    GRID_ASSERT( sdiff/norm2(out) < 1.0e-20 );
+
+    OpTest.ReleaseGrid();
+    GRID_ASSERT( MatrixChecksum(OpTest,npoint) == ckTest );  // survives release
+
+    delete CoarseMulti2;
+  }
+
+  OpTest.SetGrid(CoarseMulti);
+  GRID_ASSERT( OpTest.Nrhs() == nrhs );
+
+  OpTest.M(in,out2);
+  err = out1 - out2;
+  std::cout << GridLogMessage << "after Nrhs 4 -> 2 -> 1 -> release -> 4, |ref - test|^2 = "
+	    << norm2(err) << std::endl;
+  GRID_ASSERT( norm2(err) == 0.0 );
+
+  std::cout << GridLogMessage << "Test_coarse: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/debug/Test_coarse_coarsen.cc

@@ -0,0 +1,279 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_coarse_coarsen.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 */
+
+//
+// Coarsen the same fine operator two ways and compare the matrix elements.
+//
+//   ref  : DeprecatedMultiGeneralCoarsenedMatrix::CoarsenOperator, vectorised
+//          coarse space matched to the fine SIMD layout, single RHS fine
+//          applications
+//   test : MultiGeneralCoarsenedOperator::CoarsenOperator on the D+1 grid,
+//          unvectorised sComplexD coarse space, the batch of phased basis
+//          vectors carried in the rhs direction and applied through
+//          MrhsPromotedOperator
+//
+// BLAS_A is written by GridtoBLAS in lSite order, which does not depend on
+// the SIMD layout, so the two are directly comparable.
+//
+#include <Grid/Grid.h>
+
+using namespace Grid;
+
+const int nbasis = 8;
+const int batch  = 9;
+
+typedef vSpinColourVector FineObj;
+typedef vTComplex         CComplexV;   // vectorised coarse space, reference
+typedef sTComplexD        CComplexS;   // unvectorised coarse space, under test
+
+typedef DeprecatedMultiGeneralCoarsenedMatrix    <FineObj,CComplexV,nbasis> RefOperator;
+typedef MultiGeneralCoarsenedOperator<FineObj,CComplexS,nbasis> TestOperator;
+
+template<class Op>
+void ReadMatrix(Op &O,int npoint,std::vector<std::vector<typename Op::calcMatrix> > &host)
+{
+  typedef typename Op::calcMatrix calcMatrix;
+  host.resize(npoint);
+  // MatrixPointOut carries one stencil point in lSite order whichever internal
+  // layout the operator uses.
+  deviceVector<calcMatrix> buf;
+  for(int p=0;p<npoint;p++){
+    O.MatrixPointOut(p,buf);
+    int64_t sites = buf.size();
+    host[p].resize(sites);
+    acceleratorCopyFromDevice(&buf[0],&host[p][0],sites*sizeof(calcMatrix));
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate flatt = GridDefaultLatt();
+  Coordinate fsimd = GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate fmpi  = GridDefaultMpi();
+
+  Coordinate block({2,2,2,2});
+  Coordinate clatt(Nd);
+  for(int d=0;d<Nd;d++){
+    GRID_ASSERT(flatt[d]%block[d]==0);
+    clatt[d] = flatt[d]/block[d];
+  }
+
+  GridCartesian         *FineGrid   = new GridCartesian(flatt,fsimd,fmpi);
+  GridRedBlackCartesian *FrbGrid    = SpaceTimeGrid::makeFourDimRedBlackGrid(FineGrid);
+
+  ////////////////////////////////////////////////
+  // Reference coarse space: SIMD layout matched to the fine
+  ////////////////////////////////////////////////
+  Coordinate cvsimd = GridDefaultSimd(Nd,CComplexV::Nsimd());
+  GridCartesian *CoarseV = new GridCartesian(clatt,cvsimd,fmpi);
+
+  // The reference puts all of the SIMD in the rhs direction. CoarsenOperator
+  // does not use the multiRHS grid; it only sizes BLAS_A, so one lane of rhs suffices.
+  int nrhs_ref = CComplexV::Nsimd();
+  Coordinate reflatt(1,nrhs_ref),refsimd(1,CComplexV::Nsimd()),refmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    reflatt.push_back(clatt[d]);
+    refsimd.push_back(1);
+    refmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *CoarseVMulti = new GridCartesian(reflatt,refsimd,refmpi);
+
+  ////////////////////////////////////////////////
+  // Coarse space under test: unvectorised
+  ////////////////////////////////////////////////
+  Coordinate cssimd(Nd,1);
+  GridCartesian *CoarseS = new GridCartesian(clatt,cssimd,fmpi);
+
+  Coordinate cmlatt(1,batch),cmsimd(1,1),cmmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    cmlatt.push_back(clatt[d]);
+    cmsimd.push_back(1);
+    cmmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *CoarseSMulti = new GridCartesian(cmlatt,cmsimd,cmmpi);
+
+  // D+1 fine grid carrying the batch
+  Coordinate fmlatt(1,batch), fmsimd(1,1), fmmpi(1,1);
+  for(int d=0;d<Nd;d++){
+    fmlatt.push_back(flatt[d]);
+    fmsimd.push_back(fsimd[d]);
+    fmmpi .push_back(fmpi[d]);
+  }
+  GridCartesian *FineGridMulti = new GridCartesian(fmlatt,fmsimd,fmmpi);
+
+  std::cout << GridLogMessage << "fine "<<flatt<<"   coarse "<<clatt<<"   batch "<<batch<<std::endl;
+  std::cout << GridLogMessage << "Nsimd   fine "<<FineGrid->Nsimd()
+	    << "   coarse ref "<<CoarseV->Nsimd()
+	    << "   coarse test "<<CoarseS->Nsimd()<<std::endl;
+
+  ////////////////////////////////////////////////
+  // Fine operator
+  ////////////////////////////////////////////////
+  GridParallelRNG pRNG(FineGrid); pRNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+  LatticeGaugeFieldD Umu(FineGrid); SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  RealD mass = 0.1;
+  WilsonFermionD Dw(Umu,*FineGrid,*FrbGrid,mass);
+  MdagMLinearOperator<WilsonFermionD,LatticeFermionD> HermOp(Dw);
+
+  ////////////////////////////////////////////////
+  // One random subspace, two copies: CoarsenOperator orthogonalises in place
+  ////////////////////////////////////////////////
+  Aggregation<FineObj,CComplexV,nbasis> Subspace(CoarseV,FineGrid,0);
+  std::vector<LatticeFermionD> subspace(nbasis,FineGrid);
+  for(int i=0;i<nbasis;i++){
+    random(pRNG,Subspace.subspace[i]);
+    subspace[i] = Subspace.subspace[i];
+  }
+
+  NextToNearestStencilGeometry4D geomV(CoarseV);
+  NextToNearestStencilGeometry4D geomS(CoarseS);
+
+  ////////////////////////////////////////////////
+  // Reference coarsening, matched layouts
+  ////////////////////////////////////////////////
+  RefOperator OpRef(geomV,CoarseVMulti);
+  std::cout << GridLogMessage << "reference CoarsenOperator" << std::endl;
+  OpRef.CoarsenOperator(HermOp,Subspace,CoarseV);
+
+  ////////////////////////////////////////////////
+  // Coarsening under test, D+1 fine applications, unvectorised coarse
+  ////////////////////////////////////////////////
+  TestOperator OpTest(geomS,CoarseS);
+  OpTest.SetGrid(CoarseSMulti);
+
+  MrhsPromotedOperator<LatticeFermionD> MrhsHermOp(HermOp,FineGrid,batch);
+
+  std::cout << GridLogMessage << "test CoarsenOperator (D+1)" << std::endl;
+  OpTest.CoarsenOperator(MrhsHermOp,FineGridMulti,subspace,CoarseS);
+
+  ////////////////////////////////////////////////
+  // Compare matrix elements
+  ////////////////////////////////////////////////
+  int npoint = OpRef.geom.npoint;
+  GRID_ASSERT(npoint == OpTest.geom.npoint);
+  typedef RefOperator::calcMatrix calcMatrix;
+  std::vector<std::vector<calcMatrix> > A1,A2;
+  ReadMatrix(OpRef,npoint,A1);
+  ReadMatrix(OpTest,npoint,A2);
+
+  RealD num=0.0, den=0.0;
+  for(int p=0;p<npoint;p++){
+    GRID_ASSERT(A1[p].size()==A2[p].size());
+    ComplexD *w1 = (ComplexD *)&A1[p][0];
+    ComplexD *w2 = (ComplexD *)&A2[p][0];
+    int64_t words = A1[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      ComplexD d = w1[i]-w2[i];
+      num += real(d)*real(d)+imag(d)*imag(d);
+      den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
+    }
+  }
+  std::cout << GridLogMessage << "|A_ref|^2 = " << den << std::endl;
+  std::cout << GridLogMessage << "|A_ref - A_test|^2 / |A_ref|^2 = " << num/den << std::endl;
+  GRID_ASSERT( den > 0.0 );
+  GRID_ASSERT( num/den < 1.0e-20 );
+
+  ////////////////////////////////////////////////
+  // Same coarsening through the single RHS variant: no mrhs packing, the
+  // batch assembled on the coarse side by the mixed blockProject. Block
+  // Gram-Schmidt is idempotent so the subspace may be reused in place.
+  ////////////////////////////////////////////////
+  TestOperator OpTestSrhs(geomS,CoarseS);
+  OpTestSrhs.SetGrid(CoarseSMulti);
+
+  std::cout << GridLogMessage << "test CoarsenOperator (single RHS fine op)" << std::endl;
+  OpTestSrhs.CoarsenOperator(HermOp,subspace,CoarseS,batch);
+
+  std::vector<std::vector<calcMatrix> > A3;
+  ReadMatrix(OpTestSrhs,npoint,A3);
+
+  RealD nums=0.0;
+  for(int p=0;p<npoint;p++){
+    GRID_ASSERT(A1[p].size()==A3[p].size());
+    ComplexD *w1 = (ComplexD *)&A1[p][0];
+    ComplexD *w3 = (ComplexD *)&A3[p][0];
+    int64_t words = A1[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+    for(int64_t i=0;i<words;i++){
+      ComplexD d = w1[i]-w3[i];
+      nums += real(d)*real(d)+imag(d)*imag(d);
+    }
+  }
+  std::cout << GridLogMessage << "|A_ref - A_test_srhs|^2 / |A_ref|^2 = " << nums/den << std::endl;
+  GRID_ASSERT( nums/den < 1.0e-20 );
+
+  ////////////////////////////////////////////////
+  // GetMatrix/SetMatrix round trip through the BLAS layout array. This is
+  // how a bilingual DenseCoarseMatrix will retrieve the elements, and it
+  // needs no SetGrid since the matrix elements are Nrhs independent.
+  ////////////////////////////////////////////////
+  {
+    typedef TestOperator::CoarseMatrix CoarseMatrixS;
+    std::vector<CoarseMatrixS> Aget(npoint,CoarseS);
+    for(int p=0;p<npoint;p++) OpTest.GetMatrix(p,Aget);
+
+    TestOperator OpTestRoundTrip(geomS,CoarseS);
+    for(int p=0;p<npoint;p++) OpTestRoundTrip.SetMatrix(p,Aget);
+
+    std::vector<std::vector<calcMatrix> > A4;
+    ReadMatrix(OpTestRoundTrip,npoint,A4);
+
+    RealD numrt=0.0;
+    for(int p=0;p<npoint;p++){
+      GRID_ASSERT(A2[p].size()==A4[p].size());
+      ComplexD *w2 = (ComplexD *)&A2[p][0];
+      ComplexD *w4 = (ComplexD *)&A4[p][0];
+      int64_t words = A2[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
+      for(int64_t i=0;i<words;i++){
+        ComplexD d = w2[i]-w4[i];
+        numrt += real(d)*real(d)+imag(d)*imag(d);
+      }
+    }
+    std::cout << GridLogMessage << "GetMatrix/SetMatrix round trip |diff|^2 = " << numrt << std::endl;
+    GRID_ASSERT( numrt == 0.0 );
+  }
+
+  ////////////////////////////////////////////////
+  // and the operator under test applies the matrix it just built
+  ////////////////////////////////////////////////
+  typedef TestOperator::CoarseVector CoarseVectorS;
+  GridParallelRNG cRNG(CoarseS); cRNG.SeedFixedIntegers(std::vector<int>({5,6,7,8}));
+  CoarseVectorS in(CoarseSMulti);  random(cRNG,in);
+  CoarseVectorS out(CoarseSMulti);
+
+  OpTest.M(in,out);
+  std::cout << GridLogMessage << "|in|^2 = " << norm2(in)
+	    << "   |M_test in|^2 = " << norm2(out) << std::endl;
+  GRID_ASSERT( norm2(out) > 0.0 );
+
+  std::cout << GridLogMessage << "Test_coarse_coarsen: ALL PASS" << std::endl;
+
+  Grid_finalize();
+}

tests/debug/Test_mrhs_deflation.cc

@@ -0,0 +1,89 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_mrhs_deflation.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 */
+#include <Grid/Grid.h>
+
+// MultiRHSDeflation: the D+1 multiRHS interface (one permutation pass each
+// way) must reproduce the vector-of-D-fields interface bit for bit up to
+// the summation order of the GEMMs.  Random "eigenvectors" and values: the
+// deflation formula G = E (E^dag R)/lambda does not care that they are not
+// eigenpairs of anything.
+
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int nbasis = 8;
+  const int nev    = 12;
+  const int nrhs   = 5;
+
+  typedef iVector<sTComplexD,nbasis> siteVector;
+  typedef Lattice<siteVector>        CoarseVector;
+
+  // Unvectorised D and D+1 coarse grids, as MGCoarseGrids builds them
+  Coordinate latt({1,4,4,4,4});
+  Coordinate simd({1,1,1,1,1});
+  Coordinate mpi = GridDefaultMpi();
+  Coordinate mpi5({1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian Coarse5d(latt,simd,mpi5);
+
+  Coordinate latt6({nrhs,1,4,4,4,4});
+  Coordinate simd6({1,1,1,1,1,1});
+  Coordinate mpi6({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  GridCartesian CoarseMrhs(latt6,simd6,mpi6);
+
+  GridParallelRNG RNG(&Coarse5d); RNG.SeedFixedIntegers(std::vector<int>({1,2,3,4}));
+
+  std::vector<CoarseVector> evec(nev,&Coarse5d);
+  std::vector<RealD>        eval(nev);
+  for(int e=0;e<nev;e++){ random(RNG,evec[e]); eval[e] = 0.5 + 0.1*e; }
+
+  std::vector<CoarseVector> src(nrhs,&Coarse5d), guess(nrhs,&Coarse5d);
+  for(int r=0;r<nrhs;r++) random(RNG,src[r]);
+
+  MultiRHSDeflation<CoarseVector> Deflator;
+  Deflator.ImportEigenBasis(evec,eval);
+
+  // Reference: the vector interface
+  Deflator.DeflateSources(src,guess);
+
+  // The D+1 interface on the same sources
+  CoarseVector src_mrhs(&CoarseMrhs), guess_mrhs(&CoarseMrhs);
+  for(int r=0;r<nrhs;r++) InsertSliceFast(src[r],src_mrhs,r,0);
+  guess_mrhs = Zero();
+  Deflator.DeflateSources(src_mrhs,guess_mrhs);
+
+  RealD worst = 0.0;
+  CoarseVector g(&Coarse5d), d(&Coarse5d);
+  for(int r=0;r<nrhs;r++){
+    ExtractSliceFast(g,guess_mrhs,r,0);
+    d = g - guess[r];
+    RealD rel = std::sqrt(norm2(d)/norm2(guess[r]));
+    std::cout << GridLogMessage << "rhs " << r << " ||guess_mrhs - guess||/||guess|| = " << rel << std::endl;
+    worst = std::max(worst,rel);
+  }
+  std::cout << GridLogMessage << "MultiRHSDeflation D+1 vs vector interface: worst " << worst << std::endl;
+  GRID_ASSERT( worst < 1.0e-12 );
+  std::cout << GridLogMessage << "PASS" << std::endl;
+
+  Grid_finalize();
+  return 0;
+}

tests/debug/Test_sloppy_dagger.cc

@@ -0,0 +1,207 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_sloppy_dagger.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 */
+
+//
+// Isolate a sloppy-comms dagger defect: on a CPU/GEN multi-rank build the
+// sloppy dagger halo can be wrong by ~5% of norm2, deterministically, which
+// Benchmark_dwf's sloppy pass catches as a failed dagger Cshift check.
+//
+// The sloppy compressed-buffer pool (StencilBuffer::DeviceCommBuf) is a
+// SHARED STATIC across all stencil objects, so scenarios contaminate each
+// other inside one process: this program runs exactly ONE scenario per
+// invocation, selected with --seq <name>:
+//
+//   fresh-dag       sloppy dagger is the operator's first call
+//   nodag-dag       sloppy nodag x3 (one source), then sloppy dagger (same source)
+//   nodag-dag-fresh nodag x3 (fresh source each), then dagger (fresh source)
+//   nodag-dag-newsrc nodag x3 (one source), then dagger (different source)
+//   dag-nodag-dag   dagger, nodag x3, dagger  (dagger-first history)
+//
+// Every call is checked against a separate never-sloppy reference
+// operator (exact Dhop == the Cshift construction at 1e-31, certified by
+// Benchmark_dwf's exact pass).  The last dagger's error is fingerprinted
+// by slice along each MPI-decomposed direction.
+//
+//   OMP_NUM_THREADS=1 mpirun -n 2 ./Test_sloppy_dagger --grid 8.8.8.16 \
+//        --mpi 1.1.1.2 --seq nodag-dag
+//
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::string seq("nodag-dag");
+  if( GridCmdOptionExists(argv,argv+argc,"--seq") )
+    seq = GridCmdOptionPayload(argv,argv+argc,"--seq");
+
+  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);
+
+  GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4});
+  GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8});
+
+  LatticeGaugeField Umu(UGrid);
+  SU<Nc>::HotConfiguration(RNG4,Umu);
+
+  RealD mass=0.1, M5=1.8;
+  DomainWallFermionD Dref(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);  // never sloppy
+  DomainWallFermionD Dw  (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);  // under test
+  Dw.SloppyComms(1);
+
+  Coordinate mpi = GridDefaultMpi();
+  LatticeFermionD src(FGrid), exact(FGrid), sloppy(FGrid), err(FGrid);
+
+  // One call: fresh or reused source, dag or not; always checked.
+  LatticeFermionD srcA(FGrid); gaussian(RNG5,srcA);
+  auto call = [&](int dag, int fresh, const char *tag){
+    if ( fresh ) gaussian(RNG5,src); else src = srcA;
+    Dref.Dhop(src,exact,dag);
+    Dw.Dhop(src,sloppy,dag);
+    err = sloppy - exact;
+    std::cout << GridLogMessage << "SEQ[" << seq << "] " << tag
+              << (dag?" DAG   ":" NODAG ") << " norm2(err) " << norm2(err) << std::endl;
+  };
+  // First bad site: coordinate + exact-vs-sloppy hex words (PB: look for
+  // expected/actual mismatching in the high-order half word).
+  auto firstbad = [&](LatticeFermionD &ex, LatticeFermionD &sl){
+    typedef LatticeFermionD::scalar_object sobj;
+    Coordinate gdims = FGrid->GlobalDimensions();
+    int64_t gsites = FGrid->gSites();
+    for(int64_t g=0; g<gsites; g++){
+      Coordinate gcoor(5);
+      Lexicographic::CoorFromIndex(gcoor,g,gdims);
+      sobj se, ss_;
+      peekSite(se,ex,gcoor);
+      peekSite(ss_,sl,gcoor);
+      uint64_t *we=(uint64_t *)&se;
+      uint64_t *ws=(uint64_t *)&ss_;
+      int nw = sizeof(sobj)/8;
+      for(int w=0;w<nw;w++){
+        double de=((double *)we)[w], ds=((double *)ws)[w];
+        if ( fabs(de-ds) > 1.0e-4 ) {
+          std::cout << GridLogMessage << "FIRSTBAD site " << gcoor
+                    << " word " << w << " (spin "<<(w/6)%4<<" col "<<(w/2)%3<<" reim "<<(w%2)<<")"
+                    << "  exact " << std::hex << we[w]
+                    << "  sloppy " << ws[w] << std::dec << std::endl;
+          for(int w2=0;w2<8&&w2<nw;w2++)
+            std::cout << GridLogMessage << "   word["<<w2<<"] exact "<<std::hex<<we[w2]
+                      <<" sloppy "<<ws[w2]<<std::dec<<std::endl;
+          return;
+        }
+      }
+    }
+    std::cout << GridLogMessage << "FIRSTBAD: none above 1e-4" << std::endl;
+  };
+  auto fingerprint = [&](void){
+    for(int mu=0;mu<Nd;mu++){
+      if ( mpi[mu] == 1 ) continue;
+      std::vector<RealD> sn;
+      sliceNorm(sn,err,mu+1);
+      std::cout << GridLogMessage << "SEQ[" << seq << "] err by slice of dim " << mu << ":";
+      for(int t=0;t<(int)sn.size();t++) std::cout << " " << sn[t];
+      std::cout << std::endl;
+    }
+  };
+
+  if        ( seq == "fresh-dag" ) {
+    call(1,0,"call0");
+  } else if ( seq == "nodag-dag" ) {
+    call(0,0,"call0"); call(0,0,"call1"); call(0,0,"call2");
+    call(1,0,"call3");
+  } else if ( seq == "nodag-dag-fresh" ) {
+    call(0,1,"call0"); call(0,1,"call1"); call(0,1,"call2");
+    call(1,1,"call3");
+  } else if ( seq == "nodag-dag-newsrc" ) {
+    call(0,0,"call0"); call(0,0,"call1"); call(0,0,"call2");
+    call(1,1,"call3");
+  } else if ( seq == "dag-nodag-dag" ) {
+    call(1,0,"call0");
+    call(0,0,"call1"); call(0,0,"call2"); call(0,0,"call3");
+    call(1,0,"call4");
+  } else if ( seq == "noleave" ) {
+    // NO exact-operator call between sloppy calls: references precomputed.
+    LatticeFermionD refnodag(FGrid), refdag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);
+    Dref.Dhop(srcA,refdag,DaggerYes);
+    for(int i=0;i<2;i++){
+      Dw.Dhop(srcA,sloppy,DaggerNo);
+      err = sloppy - refnodag;
+      std::cout << GridLogMessage << "SEQ[noleave] call" << i << " NODAG  norm2(err) " << norm2(err) << std::endl;
+      if ( i==0 ) firstbad(refnodag,sloppy);
+    }
+    Dw.Dhop(srcA,sloppy,DaggerYes);
+    err = sloppy - refdag;
+    std::cout << GridLogMessage << "SEQ[noleave] call3 DAG    norm2(err) " << norm2(err) << std::endl;
+  } else if ( seq == "twoop" ) {
+    // A SECOND sloppy operator shares the static pool; sequence runs on it
+    // with interleaved exact references (as the honest runs had).
+    DomainWallFermionD Dw2(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    Dw2.SloppyComms(1);
+    Dw.Dhop(srcA,sloppy,DaggerNo);   // prime the FIRST sloppy op's state
+    for(int i=0;i<3;i++){
+      Dref.Dhop(srcA,exact,DaggerNo);
+      Dw2.Dhop(srcA,sloppy,DaggerNo);
+      err = sloppy - exact;
+      std::cout << GridLogMessage << "SEQ[twoop] call" << i << " NODAG  norm2(err) " << norm2(err) << std::endl;
+    }
+    Dref.Dhop(srcA,exact,DaggerYes);
+    Dw2.Dhop(srcA,sloppy,DaggerYes);
+    err = sloppy - exact;
+    std::cout << GridLogMessage << "SEQ[twoop] call3 DAG    norm2(err) " << norm2(err) << std::endl;
+  } else if ( seq == "self-exact" ) {
+    // Exact call on the SAME object (Dw), opposite sense, then sloppy:
+    // distinguishes per-object from shared state.
+    LatticeFermionD refnodag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);              // reference only, early
+    Dw.SloppyComms(0);
+    Dw.Dhop(srcA,sloppy,DaggerYes);                 // exact DAG on Dw itself
+    Dw.SloppyComms(1);
+    Dw.Dhop(srcA,sloppy,DaggerNo);                  // sloppy NODAG: sense mismatch
+    err = sloppy - refnodag;
+    std::cout << GridLogMessage << "SEQ[self-exact] sloppy NODAG after own exact DAG: " << norm2(err) << std::endl;
+  } else if ( seq == "cold" ) {
+    // NO exact Dhop anywhere before the sloppy calls.
+    LatticeFermionD outn(FGrid), outd(FGrid);
+    Dw.Dhop(srcA,outn,DaggerNo);
+    Dw.Dhop(srcA,outd,DaggerYes);
+    LatticeFermionD refnodag(FGrid), refdag(FGrid);
+    Dref.Dhop(srcA,refnodag,DaggerNo);
+    Dref.Dhop(srcA,refdag,DaggerYes);
+    err = outn - refnodag;
+    std::cout << GridLogMessage << "SEQ[cold] first-ever sloppy NODAG: " << norm2(err) << std::endl;
+    err = outd - refdag;
+    std::cout << GridLogMessage << "SEQ[cold] then      sloppy DAG:   " << norm2(err) << std::endl;
+  } else {
+    std::cout << GridLogMessage << "unknown --seq " << seq << std::endl;
+  }
+  fingerprint();
+
+  Grid_finalize();
+}