Coarse space is expensive, but gives a speed up in fine matrix multiplies now.

Down to optimisation

Grid/algorithms/CoarsenedMatrix.h

@@ -170,7 +170,8 @@ public:
       subspace[b] = noise;
     }
   }
-  virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
+  virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis)
+  {
 
     RealD scale;
 
@@ -229,6 +230,11 @@ public:
     scale = std::pow(norm2(noise),-0.5); 
     noise=noise*scale;
 
+    std::cout << GridLogMessage<<" Chebyshev subspace pass-1 : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
+    std::cout << GridLogMessage<<" Chebyshev subspace pass-2 : nbasis"<<nn<<" min "
+	      <<ordermin<<" step "<<orderstep
+	      <<" lo"<<filterlo<<std::endl;
+
     // Initial matrix element
     hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
 
@@ -302,6 +308,43 @@ public:
     }
     assert(b==nn);
   }
+  virtual void CreateSubspaceChebyshev(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,
+				       int nn,
+				       double hi,
+				       double lo,
+				       int orderfilter
+				       ) {
+
+    RealD scale;
+
+    FineField noise(FineGrid);
+    FineField Mn(FineGrid);
+    FineField tmp(FineGrid);
+
+    // New normalised noise
+    std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
+    std::cout << GridLogMessage<<" Chebyshev subspace pure noise  : nbasis"<<nn<<std::endl;
+
+
+    for(int b =0;b<nbasis;b++)
+    {
+      gaussian(RNG,noise);
+      scale = std::pow(norm2(noise),-0.5); 
+      noise=noise*scale;
+
+      // Initial matrix element
+      hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
+      // Filter
+      Chebyshev<FineField> Cheb(lo,hi,orderfilter);
+      Cheb(hermop,noise,Mn);
+      // normalise
+      scale = std::pow(norm2(Mn),-0.5); 	Mn=Mn*scale;
+      subspace[b]   = Mn;
+      hermop.Op(Mn,tmp); 
+      std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
+    }
+
+  }
 
 };
 

Grid/algorithms/GeneralCoarsenedMatrix.h

@@ -34,6 +34,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 
 NAMESPACE_BEGIN(Grid);
 
+// Fixme need coalesced read gpermute
 template<class vobj> void gpermute(vobj & inout,int perm){
   vobj tmp=inout;
   if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;}
@@ -50,175 +51,139 @@ template<class vobj> void gpermute(vobj & inout,int perm){
 class NonLocalStencilGeometry {
 public:
   int depth;
+  int hops;
   int npoint;
   std::vector<Coordinate> shifts;
-  virtual void BuildShifts(void) { assert(0); } ;
+  Coordinate stencil_size;
-  int Depth(void){return depth;};
+  Coordinate stencil_lo;
-  NonLocalStencilGeometry(int _depth) : depth(_depth)
+  Coordinate stencil_hi;
-  {
+  GridCartesian *grid;
-  };
+  GridCartesian *Grid() {return grid;};
+  int Depth(void){return 1;};   // Ghost zone depth
+  int Hops(void){return hops;}; // # of hops=> level of corner fill in in stencil
+
+  virtual int DimSkip(void) =0;
+
   virtual ~NonLocalStencilGeometry() {};
-};
+
-// Need to worry about red-black now
+  int  Reverse(int point)
-class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
-public:
-  NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(2)
   {
-      this->BuildShifts();
+    int Nd = Grid()->Nd();
-  };
+    Coordinate shft = shifts[point];
-  virtual ~NextToNearestStencilGeometry4D() {};
+    Coordinate rev(Nd);
-  virtual void BuildShifts(void)
+    for(int mu=0;mu<Nd;mu++) rev[mu]= -shft[mu];
+    for(int p=0;p<npoint;p++){
+      if(rev==shifts[p]){
+	return p;
+      }
+    }
+    assert(0);
+    return -1;
+  }
+  void BuildShifts(void)
   {
     this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
+    int Nd = this->grid->Nd();
-    this->shifts.push_back(Coordinate({0,0,0,0}));
+
-    // +-x, +-y, +-z, +-t  : 8
+    int dd = this->DimSkip();
-    for(int s=-1;s<=1;s+=2){
+    for(int s0=this->stencil_lo[dd+0];s0<=this->stencil_hi[dd+0];s0++){
-      this->shifts.push_back(Coordinate({s,0,0,0}));
+    for(int s1=this->stencil_lo[dd+1];s1<=this->stencil_hi[dd+1];s1++){
-      this->shifts.push_back(Coordinate({0,s,0,0}));
+    for(int s2=this->stencil_lo[dd+2];s2<=this->stencil_hi[dd+2];s2++){
-      this->shifts.push_back(Coordinate({0,0,s,0}));
+    for(int s3=this->stencil_lo[dd+3];s3<=this->stencil_hi[dd+3];s3++){
-      this->shifts.push_back(Coordinate({0,0,0,s}));
+      Coordinate sft(Nd,0);
-    }
+      sft[dd+0] = s0;
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
+      sft[dd+1] = s1;
-    for(int s1=-1;s1<=1;s1+=2){
+      sft[dd+2] = s2;
-    for(int s2=-1;s2<=1;s2+=2){
+      sft[dd+3] = s3;
-      this->shifts.push_back(Coordinate({s1,s2,0,0}));
+      int nhops = abs(s0)+abs(s1)+abs(s2)+abs(s3);
-      this->shifts.push_back(Coordinate({s1,0,s2,0}));
+      if(nhops<=this->hops) this->shifts.push_back(sft);
-      this->shifts.push_back(Coordinate({s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2}));
-    }}
-    this->npoint = this->shifts.size();
-  }
-};
-// Need to worry about red-black now
-class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
-public:
-  NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(4)
-  {
-    this->BuildShifts();
-  };
-  virtual ~NextToNextToNextToNearestStencilGeometry4D() {}
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2}));
-    }}
-    // +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-      this->shifts.push_back(Coordinate({s1,s2,s3,0})); // 8x4 = 32
-      this->shifts.push_back(Coordinate({s1,s2,0,s3}));
-      this->shifts.push_back(Coordinate({s1,0,s2,s3}));
-      this->shifts.push_back(Coordinate({0,s1,s2,s3}));
-    }}}
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-    for(int s4=-1;s4<=1;s4+=2){
-      this->shifts.push_back(Coordinate({s1,s2,s3,s4})); // 16 
     }}}}
     this->npoint = this->shifts.size();
+    std::cout << GridLogMessage << "NonLocalStencilGeometry has "<< this->npoint << " terms in stencil "<<std::endl;
   }
-};
+  
-class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
+  NonLocalStencilGeometry(GridCartesian *_coarse_grid,int _hops) : grid(_coarse_grid), hops(_hops)
-public:
-  NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(2)
   {
-      this->BuildShifts();
+    Coordinate latt = grid->GlobalDimensions();
+    stencil_size.resize(grid->Nd());
+    stencil_lo.resize(grid->Nd());
+    stencil_hi.resize(grid->Nd());
+    for(int d=0;d<grid->Nd();d++){
+     if ( latt[d] == 1 ) {
+      stencil_lo[d] = 0;
+      stencil_hi[d] = 0;
+      stencil_size[d]= 1;
+     } else if ( latt[d] == 2 ) {
+      stencil_lo[d] = -1;
+      stencil_hi[d] = 0;
+      stencil_size[d]= 2;
+     } else if ( latt[d] > 2 ) {
+       stencil_lo[d] = -1;
+       stencil_hi[d] =  1;
+       stencil_size[d]= 3;
+     }
+    }
   };
-  virtual ~NextToNearestStencilGeometry5D() {};
+
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({0,s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
-    }}
-    this->npoint = this->shifts.size();
-  }
 };
+
 // Need to worry about red-black now
-class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
+class NonLocalStencilGeometry4D : public NonLocalStencilGeometry {
 public:
-  NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(4)
+  virtual int DimSkip(void) { return 0;};
+  NonLocalStencilGeometry4D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops) { };
+  virtual ~NonLocalStencilGeometry4D() {};
+};
+class NonLocalStencilGeometry5D : public NonLocalStencilGeometry {
+public:
+  virtual int DimSkip(void) { return 1; }; 
+  NonLocalStencilGeometry5D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops)  { };
+  virtual ~NonLocalStencilGeometry5D() {};
+};
+/*
+ * Bunch of different options classes
+ */
+class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry4D {
+public:
+  NextToNextToNextToNearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,4)
   {
     this->BuildShifts();
   };
-  virtual ~NextToNextToNextToNearestStencilGeometry5D() {}
+};
-  virtual void BuildShifts(void)
+class NextToNextToNextToNearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
+public:
+  NextToNextToNextToNearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,4)
   {
-    this->shifts.resize(0);
+    this->BuildShifts();
-    // Like HDCG:  81 point stencil including self connection
+  };
-    this->shifts.push_back(Coordinate({0,0,0,0,0}));
+};
-    // +-x, +-y, +-z, +-t  : 8
+class NextToNearestStencilGeometry4D : public  NonLocalStencilGeometry4D {
-    for(int s=-1;s<=1;s+=2){
+public:
-      this->shifts.push_back(Coordinate({0,s,0,0,0}));
+  NextToNearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,2)
-      this->shifts.push_back(Coordinate({0,0,s,0,0}));
+  {
-      this->shifts.push_back(Coordinate({0,0,0,s,0}));
+    this->BuildShifts();
-      this->shifts.push_back(Coordinate({0,0,0,0,s}));
+  };
-    }
+};
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
+class NextToNearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
-    for(int s1=-1;s1<=1;s1+=2){
+public:
-    for(int s2=-1;s2<=1;s2+=2){
+  NextToNearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,2)
-      this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
+  {
-      this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
+    this->BuildShifts();
-      this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
+  };
-      this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
+};
-      this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
+class NearestStencilGeometry4D : public  NonLocalStencilGeometry4D {
-      this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
+public:
-    }}
+  NearestStencilGeometry4D(GridCartesian *Coarse) :  NonLocalStencilGeometry4D(Coarse,1)
-    // +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
+  {
-    for(int s1=-1;s1<=1;s1+=2){
+    this->BuildShifts();
-    for(int s2=-1;s2<=1;s2+=2){
+  };
-    for(int s3=-1;s3<=1;s3+=2){
+};
-      this->shifts.push_back(Coordinate({0,s1,s2,s3,0})); // 8x4 = 32
+class NearestStencilGeometry5D : public  NonLocalStencilGeometry5D {
-      this->shifts.push_back(Coordinate({0,s1,s2,0,s3}));
+public:
-      this->shifts.push_back(Coordinate({0,s1,0,s2,s3}));
+  NearestStencilGeometry5D(GridCartesian *Coarse) :  NonLocalStencilGeometry5D(Coarse,1)
-      this->shifts.push_back(Coordinate({0,0,s1,s2,s3}));
+  {
-    }}}
+    this->BuildShifts();
-    for(int s1=-1;s1<=1;s1+=2){
+  };
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-    for(int s4=-1;s4<=1;s4+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,s3,s4})); // 16 
-    }}}}
-    this->npoint = this->shifts.size();
-  }
 };
 
 // Fine Object == (per site) type of fine field
@@ -228,7 +193,7 @@ class GeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,
 public:
     
   typedef iVector<CComplex,nbasis >           siteVector;
-  typedef Lattice<CComplex >                  CoarseComplexField;
+  typedef Lattice<iScalar<CComplex> >         CoarseComplexField;
   typedef Lattice<siteVector>                 CoarseVector;
   typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
   typedef iMatrix<CComplex,nbasis >  Cobj;
@@ -239,7 +204,7 @@ public:
   // Data members
   ////////////////////
   int hermitian;
-  GridCartesian *       _FineGrid; 
+  GridBase      *       _FineGrid; 
   GridCartesian *       _CoarseGrid; 
   NonLocalStencilGeometry &geom;
   PaddedCell Cell;
@@ -251,11 +216,11 @@ public:
   ///////////////////////
   // Interface
   ///////////////////////
-  GridCartesian * Grid(void)           { return _FineGrid; };   // this is all the linalg routines need to know
+  GridBase      * Grid(void)           { return _FineGrid; };   // this is all the linalg routines need to know
-  GridCartesian * FineGrid(void)       { return _FineGrid; };   // this is all the linalg routines need to know
+  GridBase      * FineGrid(void)       { return _FineGrid; };   // this is all the linalg routines need to know
   GridCartesian * CoarseGrid(void)     { return _CoarseGrid; };   // this is all the linalg routines need to know
 
-  GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *FineGrid, GridCartesian * CoarseGrid)
+  GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid)
     : geom(_geom),
       _FineGrid(FineGrid),
       _CoarseGrid(CoarseGrid),
@@ -263,6 +228,20 @@ public:
       Cell(_geom.Depth(),_CoarseGrid),
       Stencil(Cell.grids.back(),geom.shifts)
   {
+    {
+      int npoint = _geom.npoint;
+      StencilEntry *SE;
+      autoView( Stencil_v  , Stencil, AcceleratorRead);
+      int osites=Stencil.Grid()->oSites();
+      for(int ss=0;ss<osites;ss++){
+	for(int point=0;point<npoint;point++){
+	  auto SE = Stencil_v.GetEntry(point,ss);
+	  int o = SE->_offset;
+	  assert( o< osites);
+	}
+      }    
+    }
+
     _A.resize(geom.npoint,CoarseGrid);
     _Adag.resize(geom.npoint,CoarseGrid);
   }
@@ -280,9 +259,9 @@ public:
     conformable(in.Grid(),out.Grid());
     out.Checkerboard() = in.Checkerboard();
     CoarseVector tin=in;
-    std::cout << "Calling Exchange"<<std::endl;
+
     CoarseVector pin  = Cell.Exchange(tin);
-    //    std::cout << "Called Exchange"<<std::endl;
+
     CoarseVector pout(pin.Grid());
     
     autoView( in_v , pin, AcceleratorRead);
@@ -304,10 +283,9 @@ public:
     int osites=pin.Grid()->oSites();
 
     for(int point=0;point<npoint;point++){
-      conformable(_A[point],pin);
+      conformable(A[point],pin);
     }
     
-    // Should also exchange "A" and "Adag"
     accelerator_for(sss, osites*nbasis, 1, {
       int ss = sss/nbasis;
       int b  = sss%nbasis;
@@ -318,17 +296,15 @@ public:
       int ptype;
       StencilEntry *SE;
       
-      // FIXME -- exchange the A and the A dag
-      
       for(int point=0;point<npoint;point++){
 
 	auto SE = Stencil_v.GetEntry(point,ss);
 	
 	int o = SE->_offset;
 
+	assert( o < osites);
 	// gpermute etc..
 	nbr = in_v[o];
-	assert( o< osites);
 	gpermute(nbr,SE->_permute);
 
 	for(int bb=0;bb<nbasis;bb++) {
@@ -343,44 +319,45 @@ public:
     out = Cell.Extract(pout);
   };
 
-  void Test(LinearOperatorBase<Lattice<Fobj> > &linop,
+  void PopulateAdag(void)
-	    Aggregation<Fobj,CComplex,nbasis> & Subspace)
   {
-    // Create a random
+    for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
-    GridCartesian *grid = FineGrid();
+      Coordinate bcoor;
-    FineField MbV(grid);
+      CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
-    FineField tmp(grid);
+      
-    FineField f_src(grid);
+      for(int p=0;p<geom.npoint;p++){
-    FineField f_res(grid);
+	Coordinate scoor = bcoor;
-    FineField f_ref(grid);
+	for(int mu=0;mu<bcoor.size();mu++){
-    CoarseVector c_src(CoarseGrid());
+	  int L = CoarseGrid()->GlobalDimensions()[mu];
-    CoarseVector c_res(CoarseGrid());
+	  scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
-    CoarseVector coarseInner(CoarseGrid());
+	}
-    GridParallelRNG RNG(CoarseGrid());  RNG.SeedUniqueString(std::string("Coarse RNG"));
+	// Flip to poke/peekLocalSite and not too bad
-    random(RNG,c_src);
+	auto link = peekSite(_A[p],scoor);
-    blockPromote(c_src,f_src,Subspace.subspace);
+	int pp = geom.Reverse(p);
-    linop.op(f_src,f_ref);
+	pokeSite(adj(link),_Adag[pp],bcoor);
-    this->Mult (_A,c_src,c_res);
+      }
-    blockPromote(c_res,f_res,Subspace.subspace);
+    }
-    std::cout << " GeneralCoarsenedMatrix comparison res  "<<norm2(f_res)<<std::endl;
-    std::cout << " GeneralCoarsenedMatrix comparison ref  "<<norm2(f_ref)<<std::endl;
-    f_res = f_res - f_ref;
-    std::cout << " GeneralCoarsenedMatrix comparison diff "<<norm2(f_res)<<std::endl;
   }
   void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
 		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
   {
+    RealD tproj=0.0;
+    RealD tpick=0.0;
+    RealD tmat=0.0;
+    RealD tpeek=0.0;
     std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl;
-    GridCartesian *grid = FineGrid();
+    GridBase *grid = FineGrid();
+
+    ////////////////////////////////////////////////
     // Orthogonalise the subblocks over the basis
+    ////////////////////////////////////////////////
     CoarseScalar InnerProd(CoarseGrid()); 
-    for(int b=0;b<nbasis;b++){
-      std::cout << "subspace["<<b<<"] " <<norm2(Subspace.subspace[b])<<std::endl;
-    }
     blockOrthogonalise(InnerProd,Subspace.subspace);
 
+    ////////////////////////////////////////////////
     // Now compute the matrix elements of linop between this orthonormal
     // set of vectors.
+    ////////////////////////////////////////////////
     FineField bV(grid);
     FineField MbV(grid);
     FineField tmp(grid);
@@ -393,39 +370,240 @@ public:
     for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
       Coordinate bcoor;
       CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
-      std::cout << GridLogMessage<< "CoarsenMatrix block "<< bcoor << std::endl;
+
       for(int b=0;b<nbasis;b++){
+	tpick-=usecond();
 	blockPick(CoarseGrid(),Subspace.subspace[b],bV,bcoor);
-	linop.HermOp(bV,MbV);
+	tpick+=usecond();
+	tmat-=usecond();
+	linop.Op(bV,MbV);
+	tmat+=usecond();
+	tproj-=usecond();
 	blockProject(coarseInner,MbV,Subspace.subspace);
+	tproj+=usecond();
+
+	tpeek-=usecond();
 	for(int p=0;p<geom.npoint;p++){
 	  Coordinate scoor = bcoor;
 	  for(int mu=0;mu<bcoor.size();mu++){
 	    int L = CoarseGrid()->GlobalDimensions()[mu];
 	    scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
 	  }
+	  // Flip to peekLocalSite
+	  // Flip to pokeLocalSite
 	  auto ip    = peekSite(coarseInner,scoor);
-	  std::cout << "A["<<b<<"]["<<p<<"]"<<scoor<<" "<<" shift "<<geom.shifts[p]<<" "<< ip <<std::endl;
 	  auto Ab    = peekSite(_A[p],scoor);
-	  auto Adagb = peekSite(_Adag[p],bcoor);
+	  int pp = geom.Reverse(p);
+	  auto Adagb = peekSite(_Adag[pp],bcoor);
 	  for(int bb=0;bb<nbasis;bb++){
 	    Ab(bb,b)    = ip(bb);
 	    Adagb(b,bb) = conjugate(ip(bb));
 	  }
 	  pokeSite(Ab,_A[p],scoor);
-	  pokeSite(Adagb,_Adag[p],bcoor);
+	  pokeSite(Adagb,_Adag[pp],bcoor);
+	}
+	tpeek+=usecond();
       }
     }
+    for(int p=0;p<geom.npoint;p++){
+      Coordinate coor({0,0,0,0,0});
+      auto sval = peekSite(_A[p],coor);
     }
-    std::cout << " Exchanging _A " <<std::endl;
     for(int p=0;p<geom.npoint;p++){
       _A[p] = Cell.Exchange(_A[p]);
-      _Adag[p] = Cell.Exchange(_Adag[p]);
+      _Adag[p]= Cell.Exchange(_Adag[p]);
+    }
+    std::cout << GridLogMessage<<"CoarsenOperator pick "<<tpick<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator mat  "<<tmat <<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator projection "<<tproj<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator peek/poke  "<<tpeek<<" us"<<std::endl;
+  }
+
+  /////////////////////////////////////////////////////////////
+  // 
+  // A) Only reduced flops option is to use a padded cell of depth 4
+  // and apply MpcDagMpc in the padded cell.
+  //
+  // Makes for ONE application of MpcDagMpc per vector instead of 30 or 80.
+  // With the effective cell size around (B+8)^4 perhaps 12^4/4^4 ratio
+  // Cost is 81x more, same as stencil size.
+  //
+  // But: can eliminate comms and do as local dirichlet.
+  //
+  // Local exchange gauge field once.
+  // Apply to all vectors, local only computation.
+  // Must exchange ghost subcells in reverse process of PaddedCell to take inner products
+  //
+  // B) Can reduce cost: pad by 1, apply Deo      (4^4+6^4+8^4+8^4 )/ (4x 4^4)
+  //                     pad by 2, apply Doe
+  //                     pad by 3, apply Deo
+  //                     then break out 8x directions; cost is ~10x MpcDagMpc per vector
+  //
+  // => almost factor of 10 in setup cost, excluding data rearrangement
+  //
+  // Intermediates -- ignore the corner terms, leave approximate and force Hermitian
+  // Intermediates -- pad by 2 and apply 1+8+24 = 33 times.
+  /////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////
+    // BFM HDCG style approach: Solve a system of equations to get Aij
+    //////////////////////////////////////////////////////////
+    /*
+     *     Here, k,l index which possible shift within the 3^Nd "ball" connected by MdagM.
+     *
+     *     conj(phases[block]) proj[k][ block*Nvec+j ] =  \sum_ball  e^{i q_k . delta} < phi_{block,j} | MdagM | phi_{(block+delta),i} > 
+     *                                                 =  \sum_ball e^{iqk.delta} A_ji
+     *
+     *     Must invert matrix M_k,l = e^[i q_k . delta_l]
+     *
+     *     Where q_k = delta_k . (2*M_PI/global_nb[mu])
+     */
+  void CoarsenOperatorColoured(LinearOperatorBase<Lattice<Fobj> > &linop,
+			       Aggregation<Fobj,CComplex,nbasis> & Subspace)
+  {
+    std::cout << GridLogMessage<< "CoarsenMatrixColoured "<< std::endl;
+    GridBase *grid = FineGrid();
+
+    RealD tproj=0.0;
+    RealD teigen=0.0;
+    RealD tmat=0.0;
+    RealD tphase=0.0;
+    RealD tinv=0.0;
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid()); 
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    const int npoint = geom.npoint;
+      
+    Coordinate clatt = CoarseGrid()->GlobalDimensions();
+    int Nd = CoarseGrid()->Nd();
+
+      /*
+       *     Here, k,l index which possible momentum/shift within the N-points connected by MdagM.
+       *     Matrix index i is mapped to this shift via 
+       *               geom.shifts[i]
+       *
+       *     conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] 
+       *       =  \sum_{l in ball}  e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > 
+       *       =  \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l}
+       *       = M_{kl} A_ji^{b.b+l}
+       *
+       *     Must assemble and invert matrix M_k,l = e^[i q_k . delta_l]
+       *  
+       *     Where q_k = delta_k . (2*M_PI/global_nb[mu])
+       *
+       *     Then A{ji}^{b,b+l} = M^{-1}_{lm} ComputeProj_{m,b,i,j}
+       */
+    teigen-=usecond();
+    ComplexD ci(0.0,1.0);
+    Eigen::MatrixXcd Mkl    = Eigen::MatrixXcd::Zero(npoint,npoint);
+    Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint);
+    for(int k=0;k<npoint;k++){ // Loop over momenta
+
+      for(int l=0;l<npoint;l++){ // Loop over nbr relative
+	std::complex<double> phase(0.0,0.0);
+	for(int mu=0;mu<Nd;mu++){
+	  RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	  phase=phase+TwoPiL*geom.shifts[k][mu]*geom.shifts[l][mu];
+	}
+	phase=exp(phase*ci);
+	Mkl(k,l) = phase;
       }
     }
+    invMkl = Mkl.inverse();
+    teigen+=usecond();
+
+    ///////////////////////////////////////////////////////////////////////
+    // Now compute the matrix elements of linop between the orthonormal
+    // set of vectors.
+    ///////////////////////////////////////////////////////////////////////
+    FineField phaV(grid); // Phased block basis vector
+    FineField MphaV(grid);// Matrix applied
+    CoarseVector coarseInner(CoarseGrid());
+
+    std::vector<CoarseVector> ComputeProj(npoint,CoarseGrid());
+    std::vector<CoarseVector>          FT(npoint,CoarseGrid());
+    for(int i=0;i<nbasis;i++){// Loop over basis vectors
+      std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<"/"<<nbasis<< std::endl;
+      for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
+	/////////////////////////////////////////////////////
+	// Stick a phase on every block
+	/////////////////////////////////////////////////////
+	tphase-=usecond();
+	CoarseComplexField coor(CoarseGrid());
+	CoarseComplexField pha(CoarseGrid());	pha=Zero();
+	for(int mu=0;mu<Nd;mu++){
+	  LatticeCoordinate(coor,mu);
+	  RealD TwoPiL =  M_PI * 2.0/ clatt[mu];
+	  pha = pha + (TwoPiL * geom.shifts[p][mu]) * coor;
+	}
+	pha  =exp(pha*ci);
+	phaV=Zero();
+	blockZAXPY(phaV,pha,Subspace.subspace[i],phaV);
+	tphase+=usecond();
+
+	/////////////////////////////////////////////////////////////////////
+	// Multiple phased subspace vector by matrix and project to subspace
+	// Remove local bulk phase to leave relative phases
+	/////////////////////////////////////////////////////////////////////
+	tmat-=usecond();
+	linop.Op(phaV,MphaV);
+	tmat+=usecond();
+
+	tproj-=usecond();
+	blockProject(coarseInner,MphaV,Subspace.subspace);
+	coarseInner = conjugate(pha) * coarseInner;
+
+	ComputeProj[p] = coarseInner;
+	tproj+=usecond();
+
+      }
+
+      tinv-=usecond();
+      for(int k=0;k<npoint;k++){
+	FT[k] = Zero();
+	for(int l=0;l<npoint;l++){
+	  FT[k]= FT[k]+ invMkl(l,k)*ComputeProj[l];
+	}
+      
+	int osites=CoarseGrid()->oSites();
+	autoView( A_v  , _A[k], AcceleratorWrite);
+	autoView( FT_v  , FT[k], AcceleratorRead);
+	accelerator_for(sss, osites, 1, {
+	    for(int j=0;j<nbasis;j++){
+	      A_v[sss](j,i) = FT_v[sss](j);
+	    }
+        });
+      }
+      tinv+=usecond();
+    }
+
+    for(int p=0;p<geom.npoint;p++){
+      Coordinate coor({0,0,0,0,0});
+      auto sval = peekSite(_A[p],coor);
+    }
+
+    PopulateAdag();
+
+    // Need to write something to populate Adag from A
+    for(int p=0;p<geom.npoint;p++){
+      _A[p] = Cell.Exchange(_A[p]);
+      _Adag[p]= Cell.Exchange(_Adag[p]);
+    }
+    std::cout << GridLogMessage<<"CoarsenOperator eigen  "<<teigen<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator phase  "<<tphase<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator mat    "<<tmat <<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator proj   "<<tproj<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"CoarsenOperator inv    "<<tinv<<" us"<<std::endl;
+  }
+
   virtual  void Mdiag    (const Field &in, Field &out){ assert(0);};
   virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);};
   virtual  void MdirAll  (const Field &in, std::vector<Field> &out){assert(0);};
 };
 
+
 NAMESPACE_END(Grid);

Grid/algorithms/iterative/AdefGeneric.h

@@ -33,15 +33,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
    * Script A = SolverMatrix 
    * Script P = Preconditioner
    *
-   * Deflation methods considered
-   *      -- Solve P A x = P b        [ like Luscher ]
-   * DEF-1        M P A x = M P b     [i.e. left precon]
-   * DEF-2        P^T M A x = P^T M b
-   * ADEF-1       Preconditioner = M P + Q      [ Q + M + M A Q]
-   * ADEF-2       Preconditioner = P^T M + Q
-   * BNN          Preconditioner = P^T M P + Q
-   * BNN2         Preconditioner = M P + P^TM +Q - M P A M 
-   * 
    * Implement ADEF-2
    *
    * Vstart = P^Tx + Qb
@@ -49,45 +40,157 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
    * M2=M3=1
    * Vout = x
    */
+NAMESPACE_BEGIN(Grid);
 
-// abstract base
+template<class Field, class CoarseField, class Aggregation>
-template<class Field, class CoarseField>
 class TwoLevelFlexiblePcg : public LinearFunction<Field>
 {
  public:
-  int verbose;
   RealD   Tolerance;
   Integer MaxIterations;
-  const int mmax = 5;
+  const int mmax = 1;
   GridBase *grid;
   GridBase *coarsegrid;
 
-  LinearOperatorBase<Field>   *_Linop
+  // Fine operator, Smoother, CoarseSolver
-  OperatorFunction<Field>     *_Smoother,
+  LinearOperatorBase<Field>   &_FineLinop;
-  LinearFunction<CoarseField> *_CoarseSolver;
+  LinearFunction<Field>   &_Smoother;
+  LinearFunction<CoarseField> &_CoarseSolver;
+  LinearFunction<CoarseField> &_CoarseSolverPrecise;
 
-  // Need somthing that knows how to get from Coarse to fine and back again
+  // Need something that knows how to get from Coarse to fine and back again
+  //  void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
+  //  void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
+  Aggregation &_Aggregates;                    
   
   // more most opertor functions
   TwoLevelFlexiblePcg(RealD tol,
 		      Integer maxit,
-		     LinearOperatorBase<Field> *Linop,
+		      LinearOperatorBase<Field>   &FineLinop,
-		     LinearOperatorBase<Field> *SmootherLinop,
+		      LinearFunction<Field>   &Smoother,
-		     OperatorFunction<Field>   *Smoother,
+		      LinearFunction<CoarseField>  &CoarseSolver,
-		     OperatorFunction<CoarseField>  CoarseLinop
+		      LinearFunction<CoarseField>  &CoarseSolverPrecise,
+		      Aggregation &Aggregates
 		     ) : 
       Tolerance(tol), 
       MaxIterations(maxit),
-      _Linop(Linop),
+      _FineLinop(FineLinop),
-      _PreconditionerLinop(PrecLinop),
+      _Smoother(Smoother),
-      _Preconditioner(Preconditioner)
+      _CoarseSolver(CoarseSolver),
+      _CoarseSolverPrecise(CoarseSolverPrecise),
+      _Aggregates(Aggregates)
   {
-    verbose=0;
+    coarsegrid = Aggregates.CoarseGrid;
+    grid       = Aggregates.FineGrid;
   };
 
+  void Inflexible(Field &src,Field &psi)
+  {
+    Field resid(grid);
+    RealD f;
+    RealD rtzp,rtz,a,d,b;
+    RealD rptzp;
+    
+    Field x(grid); 
+    Field p(grid);
+    Field z(grid);
+    Field tmp(grid);
+    Field mmp(grid);
+    Field r  (grid);
+    Field mu (grid);
+    Field rp (grid);
+
+    //Initial residual computation & set up
+    RealD guess = norm2(psi);
+    double tn;
+
+    //////////////////////////
+    // x0 = Vstart -- possibly modify guess
+    //////////////////////////
+    x=Zero();
+    Vstart(x,src);
+
+    // r0 = b -A x0
+    _FineLinop.HermOp(x,mmp);
+
+    axpy(r, -1.0, mmp, src);    // Recomputes r=src-x0
+    rp=r;
+
+    //////////////////////////////////
+    // Compute z = M1 x
+    //////////////////////////////////
+    PcgM1(r,z);
+    rtzp =real(innerProduct(r,z));
+
+    ///////////////////////////////////////
+    // Except Def2, M2 is trivial
+    ///////////////////////////////////////
+    p=z;
+
+    RealD ssq =  norm2(src);
+    RealD rsq =  ssq*Tolerance*Tolerance;
+
+    std::cout<<GridLogMessage<<"HDCG: k=0 residual "<<rtzp<<" target rsq "<<rsq<<" ssq "<<ssq<<std::endl;
+    
+    for (int k=1;k<=MaxIterations;k++){
+
+      rtz=rtzp;
+      d= PcgM3(p,mmp);
+      a = rtz/d;
+
+      axpy(x,a,p,x);
+      RealD rn = axpy_norm(r,-a,mmp,r);
+
+      PcgM1(r,z);
+
+      rtzp =real(innerProduct(r,z));
+
+      int ipcg=1; // almost free inexact preconditioned CG
+      if (ipcg) {
+	rptzp =real(innerProduct(rp,z));
+      } else {
+	rptzp =0;
+      }
+      b = (rtzp-rptzp)/rtz;
+
+      PcgM2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
+
+      axpy(p,b,p,mu);  // mu = A r
+
+      RealD rrn=sqrt(rn/ssq);
+      RealD rtn=sqrt(rtz/ssq);
+      std::cout<<GridLogMessage<<"HDCG: Pcg k= "<<k<<" residual = "<<rrn<<std::endl;
+
+      if ( ipcg ) {
+	axpy(rp,0.0,r,r);
+      }
+
+      // Stopping condition
+      if ( rn <= rsq ) { 
+
+	std::cout<<GridLogMessage<<"HDCG: Pcg converged in "<<k<<" iterations"<<std::endl;;
+
+	_FineLinop.HermOp(x,mmp);			  
+	axpy(tmp,-1.0,src,mmp);
+
+	RealD  mmpnorm = sqrt(norm2(mmp));
+	RealD  psinorm = sqrt(norm2(x));
+	RealD  srcnorm = sqrt(norm2(src));
+	RealD  tmpnorm = sqrt(norm2(tmp));
+	RealD  true_residual = tmpnorm/srcnorm;
+	std::cout<<GridLogMessage<<"HDCG: true residual is "<<true_residual
+		 <<" solution "<<psinorm<<" source "<<srcnorm<<std::endl;
+
+	return;
+      }
+
+    }
+    std::cout << "HDCG: Pcg not converged"<<std::endl;
+    return ;
+  }
+  
   // The Pcg routine is common to all, but the various matrices differ from derived 
   // implementation to derived implmentation
-  void operator() (const Field &src, Field &psi){
   void operator() (const Field &src, Field &psi){
 
     psi.Checkerboard() = src.Checkerboard();
@@ -108,7 +211,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     std::vector<Field> mmp(mmax,grid);
     std::vector<RealD> pAp(mmax);
 
-    Field x  (grid); x = psi;
+    Field x  (grid);
     Field z  (grid);
     Field tmp(grid);
     Field r  (grid);
@@ -117,25 +220,23 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     //////////////////////////
     // x0 = Vstart -- possibly modify guess
     //////////////////////////
-    x=src;
+    x=Zero();
     Vstart(x,src);
 
     // r0 = b -A x0
-    HermOp(x,mmp); // Shouldn't this be something else?
+    _FineLinop.HermOp(x,mmp[0]); // Fine operator
     axpy (r, -1.0,mmp[0], src);    // Recomputes r=src-Ax0
 
     //////////////////////////////////
-    // Compute z = M1 x
+    // Compute z = M1 r
     //////////////////////////////////
-    M1(r,z,tmp,mp,SmootherMirs);
+    PcgM1(r,z);
     rtzp =real(innerProduct(r,z));
 
     ///////////////////////////////////////
     // Solve for Mss mu = P A z and set p = z-mu
-    // Def2: p = 1 - Q Az = Pright z 
-    // Other algos M2 is trivial
     ///////////////////////////////////////
-    M2(z,p[0]);
+    PcgM2(z,p[0]);
 
     for (int k=0;k<=MaxIterations;k++){
     
@@ -143,26 +244,38 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
       int peri_kp = (k+1) % mmax;
 
       rtz=rtzp;
-      d= M3(p[peri_k],mp,mmp[peri_k],tmp);
+      d= PcgM3(p[peri_k],mmp[peri_k]);
       a = rtz/d;
     
       // Memorise this
       pAp[peri_k] = d;
+      std::cout << GridLogMessage << " pCG d "<< d<<std::endl;
 
       axpy(x,a,p[peri_k],x);
+      //      std::cout << GridLogMessage << " pCG x "<< norm2(x)<<std::endl;
       RealD rn = axpy_norm(r,-a,mmp[peri_k],r);
 
+      std::cout << GridLogMessage << " pCG rn "<< rn<<std::endl;
+
       // Compute z = M x
-      M1(r,z,tmp,mp);
+      PcgM1(r,z);
+      //      std::cout << GridLogMessage << " pCG z "<< norm2(z)<<std::endl;
 
       rtzp =real(innerProduct(r,z));
+      std::cout << GridLogMessage << " pCG rtzp "<<rtzp<<std::endl;
+      //      std::cout << GridLogMessage << " pCG r "<<norm2(r)<<std::endl;
 
-      M2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
+      PcgM2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
 
-      p[peri_kp]=p[peri_k];
+      //      std::cout << GridLogMessage << " pCG mu "<<norm2(mu)<<std::endl;
       
-      // Standard search direction  p -> z + b p    ; b = 
+      p[peri_kp]=mu;
+
+      //      std::cout << GridLogMessage << " pCG p[peri_kp] "<<norm2(p[peri_kp])<<std::endl;
+
+      // Standard search direction  p -> z + b p 
       b = (rtzp)/rtz;
+      std::cout << GridLogMessage << " pCG b "<< b<<std::endl;
 
       int northog;
       //    northog     = (peri_kp==0)?1:peri_kp; // This is the fCG(mmax) algorithm
@@ -174,6 +287,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
 	RealD beta = -pbApk/pAp[peri_back];
 	axpy(p[peri_kp],beta,p[peri_back],p[peri_kp]);
       }
+      //      std::cout << GridLogMessage << " pCG p[peri_kp] orthog "<< norm2(p[peri_kp])<<std::endl;
 
       RealD rrn=sqrt(rn/ssq);
       std::cout<<GridLogMessage<<"TwoLevelfPcg: k= "<<k<<" residual = "<<rrn<<std::endl;
@@ -181,7 +295,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
       // Stopping condition
       if ( rn <= rsq ) { 
 
-	HermOp(x,mmp); // Shouldn't this be something else?
+	_FineLinop.HermOp(x,mmp[0]); // Shouldn't this be something else?
 	axpy(tmp,-1.0,src,mmp[0]);
 	
 	RealD psinorm = sqrt(norm2(x));
@@ -190,7 +304,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
 	RealD true_residual = tmpnorm/srcnorm;
 	std::cout<<GridLogMessage<<"TwoLevelfPcg:   true residual is "<<true_residual<<std::endl;
 	std::cout<<GridLogMessage<<"TwoLevelfPcg: target residual was"<<Tolerance<<std::endl;
-	return k;
+	return;
       }
     }
     // Non-convergence
@@ -199,52 +313,42 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
 
  public:
 
-  virtual void M(Field & in,Field & out,Field & tmp) {
+  virtual void PcgM1(Field & in, Field & out)
-
+  {
-  }
-
-  virtual void M1(Field & in, Field & out) {// the smoother
-
     // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
+
     Field tmp(grid);
     Field Min(grid);
+    CoarseField PleftProj(coarsegrid);
+    CoarseField PleftMss_proj(coarsegrid);
 
-    PcgM(in,Min); // Smoother call
+    _Smoother(in,Min);
 
-    HermOp(Min,out);
+    _FineLinop.HermOp(Min,out);
     axpy(tmp,-1.0,out,in);          // tmp  = in - A Min
 
-    ProjectToSubspace(tmp,PleftProj);     
+    _Aggregates.ProjectToSubspace(PleftProj,tmp);     
-    ApplyInverse(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
+    _CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
-    PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]  
+    _Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]  
+
     axpy(out,1.0,Min,tmp); // Min+tmp
+
   }
 
-  virtual void M2(const Field & in, Field & out) {
+  virtual void PcgM2(const Field & in, Field & out) {
     out=in;
-    // Must override for Def2 only
-    //  case PcgDef2:
-    //    Pright(in,out);
-    //    break;
   }
 
-  virtual RealD M3(const Field & p, Field & mmp){
+  virtual RealD PcgM3(const Field & p, Field & mmp){
-    double d,dd;
+    RealD dd;
-    HermOpAndNorm(p,mmp,d,dd);
+    _FineLinop.HermOp(p,mmp);
+    ComplexD dot = innerProduct(p,mmp);
+    dd=real(dot);
     return dd;
-    // Must override for Def1 only
-    //  case PcgDef1:
-    //    d=linop_d->Mprec(p,mmp,tmp,0,1);// Dag no
-    //      linop_d->Mprec(mmp,mp,tmp,1);// Dag yes
-    //    Pleft(mp,mmp);
-    //    d=real(linop_d->inner(p,mmp));
   }
 
-  virtual void VstartDef2(Field & xconst Field & src){
+  virtual void Vstart(Field & x,const Field & src)
-    //case PcgDef2:
+  {
-    //case PcgAdef2: 
-    //case PcgAdef2f:
-    //case PcgV11f:
     ///////////////////////////////////
     // Choose x_0 such that 
     // x_0 = guess +  (A_ss^inv) r_s = guess + Ass_inv [src -Aguess]
@@ -258,140 +362,22 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     ///////////////////////////////////
     Field r(grid);
     Field mmp(grid);
+    CoarseField PleftProj(coarsegrid);
+    CoarseField PleftMss_proj(coarsegrid);
 
-    HermOp(x,mmp);
+    _Aggregates.ProjectToSubspace(PleftProj,src);     
-    axpy (r, -1.0, mmp, src);        // r_{-1} = src - A x
+    _CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s
-    ProjectToSubspace(r,PleftProj);     
+    _Aggregates.PromoteFromSubspace(PleftMss_proj,x);  
-    ApplyInverseCG(PleftProj,PleftMss_proj); // Ass^{-1} r_s
-    PromoteFromSubspace(PleftMss_proj,mmp);  
-    x=x+mmp;
 
   }
 
-  virtual void Vstart(Field & x,const Field & src){
-    return;
-  }
-
   /////////////////////////////////////////////////////////////////////
-  // Only Def1 has non-trivial Vout. Override in Def1
+  // Only Def1 has non-trivial Vout.
   /////////////////////////////////////////////////////////////////////
   virtual void   Vout  (Field & in, Field & out,Field & src){
     out = in;
-    //case PcgDef1:
-    //    //Qb + PT x
-    //    ProjectToSubspace(src,PleftProj);     
-    //    ApplyInverse(PleftProj,PleftMss_proj); // Ass^{-1} r_s
-    //    PromoteFromSubspace(PleftMss_proj,tmp);  
-    //    
-    //    Pright(in,out);
-    //    
-    //    linop_d->axpy(out,tmp,out,1.0);
-    //    break;
   }
+};
 
-  ////////////////////////////////////////////////////////////////////////////////////////////////
+NAMESPACE_END(Grid);
-  // Pright and Pleft are common to all implementations
-  ////////////////////////////////////////////////////////////////////////////////////////////////
-  virtual void Pright(Field & in,Field & out){
-    // P_R  = [ 1              0 ] 
-    //        [ -Mss^-1 Msb    0 ] 
-    Field in_sbar(grid);
-
-    ProjectToSubspace(in,PleftProj);     
-    PromoteFromSubspace(PleftProj,out);  
-    axpy(in_sbar,-1.0,out,in);       // in_sbar = in - in_s 
-
-    HermOp(in_sbar,out);
-    ProjectToSubspace(out,PleftProj);           // Mssbar in_sbar  (project)
-
-    ApplyInverse     (PleftProj,PleftMss_proj); // Mss^{-1} Mssbar 
-    PromoteFromSubspace(PleftMss_proj,out);     // 
-
-    axpy(out,-1.0,out,in_sbar);     // in_sbar - Mss^{-1} Mssbar in_sbar
-  }
-  virtual void Pleft (Field & in,Field & out){
-    // P_L  = [ 1  -Mbs Mss^-1] 
-    //        [ 0   0         ] 
-    Field in_sbar(grid);
-    Field    tmp2(grid);
-    Field    Mtmp(grid);
-
-    ProjectToSubspace(in,PleftProj);     
-    PromoteFromSubspace(PleftProj,out);  
-    axpy(in_sbar,-1.0,out,in);      // in_sbar = in - in_s
-
-    ApplyInverse(PleftProj,PleftMss_proj); // Mss^{-1} in_s
-    PromoteFromSubspace(PleftMss_proj,out);
-
-    HermOp(out,Mtmp);
-
-    ProjectToSubspace(Mtmp,PleftProj);      // Msbar s Mss^{-1}
-    PromoteFromSubspace(PleftProj,tmp2);
-
-    axpy(out,-1.0,tmp2,Mtmp);
-    axpy(out,-1.0,out,in_sbar);     // in_sbar - Msbars Mss^{-1} in_s
-  }
-}
-
-template<class Field>
-class TwoLevelFlexiblePcgADef2 : public TwoLevelFlexiblePcg<Field> {
- public:
-  virtual void M(Field & in,Field & out,Field & tmp){
-
-  } 
-  virtual void M1(Field & in, Field & out,Field & tmp,Field & mp){
-
-  }
-  virtual void M2(Field & in, Field & out){
-
-  }
-  virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp){
-
-  }
-  virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp){
-
-  }
-}
-/*
-template<class Field>
-class TwoLevelFlexiblePcgAD : public TwoLevelFlexiblePcg<Field> {
- public:
-  virtual void M(Field & in,Field & out,Field & tmp); 
-  virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
-  virtual void M2(Field & in, Field & out);
-  virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
-  virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
-}
-
-template<class Field>
-class TwoLevelFlexiblePcgDef1 : public TwoLevelFlexiblePcg<Field> {
- public:
-  virtual void M(Field & in,Field & out,Field & tmp); 
-  virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
-  virtual void M2(Field & in, Field & out);
-  virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
-  virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
-  virtual void   Vout  (Field & in, Field & out,Field & src,Field & tmp);
-}
-
-template<class Field>
-class TwoLevelFlexiblePcgDef2 : public TwoLevelFlexiblePcg<Field> {
- public:
-  virtual void M(Field & in,Field & out,Field & tmp); 
-  virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
-  virtual void M2(Field & in, Field & out);
-  virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
-  virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
-}
-
-template<class Field>
-class TwoLevelFlexiblePcgV11: public TwoLevelFlexiblePcg<Field> {
- public:
-  virtual void M(Field & in,Field & out,Field & tmp); 
-  virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
-  virtual void M2(Field & in, Field & out);
-  virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
-  virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
-}
-*/
 #endif

Grid/algorithms/iterative/ConjugateGradient.h

@@ -183,13 +183,13 @@ public:
 		  << "\tTrue residual " << true_residual
 		  << "\tTarget " << Tolerance << std::endl;
 
-        std::cout << GridLogMessage << "Time breakdown "<<std::endl;
 	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
-	std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
 	std::cout << GridLogDebug << "\tMobius flop rate " << DwfFlops/ usecs<< " Gflops " <<std::endl;
 

Grid/algorithms/iterative/NormalEquations.h

@@ -33,7 +33,7 @@ NAMESPACE_BEGIN(Grid);
 ///////////////////////////////////////////////////////////////////////////////////////////////////////
 // Take a matrix and form an NE solver calling a Herm solver
 ///////////////////////////////////////////////////////////////////////////////////////////////////////
-template<class Field> class NormalEquations {
+template<class Field> class NormalEquations : public LinearFunction<Field>{
 private:
   SparseMatrixBase<Field> & _Matrix;
   OperatorFunction<Field> & _HermitianSolver;
@@ -60,7 +60,7 @@ public:
   }     
 };
 
-template<class Field> class HPDSolver {
+template<class Field> class HPDSolver : public LinearFunction<Field> {
 private:
   LinearOperatorBase<Field> & _Matrix;
   OperatorFunction<Field> & _HermitianSolver;
@@ -84,7 +84,7 @@ public:
 };
 
 
-template<class Field> class MdagMSolver {
+template<class Field> class MdagMSolver : public LinearFunction<Field> {
 private:
   SparseMatrixBase<Field> & _Matrix;
   OperatorFunction<Field> & _HermitianSolver;

Grid/algorithms/iterative/PowerMethod.h

@@ -20,7 +20,7 @@ template<class Field> class PowerMethod
     RealD evalMaxApprox = 0.0; 
     auto src_n = src; 
     auto tmp = src; 
-    const int _MAX_ITER_EST_ = 50; 
+    const int _MAX_ITER_EST_ = 100; 
 
     for (int i=0;i<_MAX_ITER_EST_;i++) { 
       

Grid/lattice/Lattice_rng.h

@@ -361,9 +361,14 @@ public:
     _bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
     _uid.resize(_vol,std::uniform_int_distribution<uint32_t>() );
   }
-
+  template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist)
-  template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist){
+  {
-
+    if ( l.Grid()->_isCheckerBoarded ) {
+      Lattice<vobj> tmp(_grid);
+      fill(tmp,dist);
+      pickCheckerboard(l.Checkerboard(),l,tmp);
+      return;
+    }
     typedef typename vobj::scalar_object scalar_object;
     typedef typename vobj::scalar_type scalar_type;
     typedef typename vobj::vector_type vector_type;

Grid/util/Coordinate.h

@@ -94,6 +94,13 @@ static constexpr int MaxDims = GRID_MAX_LATTICE_DIMENSION;
 
 typedef AcceleratorVector<int,MaxDims> Coordinate;
 
+template<class T,int _ndim>
+inline bool operator==(const AcceleratorVector<T,_ndim> &v,const AcceleratorVector<T,_ndim> &w)
+{
+  if (v.size()!=w.size()) return false;
+  for(int i=0;i<v.size();i++) if ( v[i]!=w[i] ) return false;
+  return true;
+}
 template<class T,int _ndim>
 inline std::ostream & operator<<(std::ostream &os, const AcceleratorVector<T,_ndim> &v)
 {

tests/debug/Test_general_coarse.cc

@@ -37,6 +37,9 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;
 
+///////////////////////
+// Tells little dirac op to use MdagM as the .Op()
+///////////////////////
 template<class Field>
 class HermOpAdaptor : public LinearOperatorBase<Field>
 {
@@ -56,84 +59,6 @@ public:
   void HermOp(const Field &in, Field &out){
     wrapped.HermOp(in,out);
   }
-  
-};
-
-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(tmp,out);
-    _Mat.Mdag(in,tmp);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp"<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    _PV.M(out,tmp);
-    _Mat.Mdag(tmp,out);
-    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-    
-  }
-};
-
-template<class Field> class DumbOperator  : public LinearOperatorBase<Field> {
-public:
-  LatticeComplex scale;
-  DumbOperator(GridBase *grid) : scale(grid)
-  {
-    scale = 0.0;
-    LatticeComplex scalesft(grid);
-    LatticeComplex scaletmp(grid);
-    for(int d=0;d<4;d++){
-      Lattice<iScalar<vInteger> > x(grid); LatticeCoordinate(x,d+1);
-      LatticeCoordinate(scaletmp,d+1);
-      scalesft = Cshift(scaletmp,d+1,1);
-      scale = 100.0*scale + where( mod(x    ,2)==(Integer)0, scalesft,scaletmp);
-    }
-    std::cout << " scale\n" << scale << std::endl;
-  }
-  // Support for coarsening to a multigrid
-  void OpDiag (const Field &in, Field &out) {};
-  void OpDir  (const Field &in, Field &out,int dir,int disp){};
-  void OpDirAll  (const Field &in, std::vector<Field> &out) {};
-
-  void Op     (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void AdjOp  (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void HermOp(const Field &in, Field &out){
-    double n1, n2;
-    HermOpAndNorm(in,out,n1,n2);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
-    ComplexD dot;
-
-    out = scale * in;
-
-    dot= innerProduct(in,out);
-    n1=real(dot);
-
-    dot = innerProduct(out,out);
-    n2=real(dot);
-  }
 };
 
 
@@ -143,7 +68,9 @@ int main (int argc, char ** argv)
 
   const int Ls=4;
 
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
 
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
@@ -154,7 +81,10 @@ int main (int argc, char ** argv)
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/2;
   }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
+
+  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt,
+							    GridDefaultSimd(Nd,vComplex::Nsimd()),
+							    GridDefaultMpi());;
   GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
 
   std::vector<int> seeds4({1,2,3,4});
@@ -170,19 +100,15 @@ int main (int argc, char ** argv)
   LatticeFermion    tmp(FGrid);
   LatticeFermion    err(FGrid);
   LatticeGaugeField Umu(UGrid);
-  //SU<Nc>::HotConfiguration(RNG4,Umu);
+  SU<Nc>::HotConfiguration(RNG4,Umu);
-  SU<Nc>::ColdConfiguration(Umu);
+  //  Umu=Zero();
-  //  auto U = peekLorentz(Umu,0);
-  //  Umu=Zero(); // Make operator local for now
-  //  pokeLorentz(Umu,U,0);
   
-  RealD mass=0.5;
+  RealD mass=0.1;
   RealD M5=1.8;
 
   DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
 
-  const int nbasis = 20;
+  const int nbasis = 4;
   const int cb = 0 ;
   LatticeFermion prom(FGrid);
 
@@ -190,40 +116,52 @@ int main (int argc, char ** argv)
 
   std::cout<<GridLogMessage<<"Calling Aggregation class" <<std::endl;
 
-  // Possible tactics -- with zero gauge field, verify block locality of dirac op
+  ///////////////////////////////////////////////////////////
-  // Possible tactics -- with zero gauge field, take inner products
+  // Squared operator is in HermOp
-  
+  ///////////////////////////////////////////////////////////
-  // Squared operator
   MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
-  DumbOperator<LatticeFermion> Diagonal(FGrid);
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
 
+  ///////////////////////////////////////////////////
+  // Random aggregation space
+  ///////////////////////////////////////////////////
+  std::cout<<GridLogMessage << "Building random aggregation class"<< std::endl;
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
   Subspace Aggregates(Coarse5d,FGrid,cb);
   Aggregates.CreateSubspaceRandom(RNG5);
 
-  std::cout<<GridLogMessage << "Called aggregation class"<< std::endl;
+  ///////////////////////////////////////////////////
+  // Build little dirac op
+  ///////////////////////////////////////////////////
+  std::cout<<GridLogMessage << "Building little Dirac operator"<< std::endl;
 
   typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
-  NextToNearestStencilGeometry5D geom;
+  NextToNearestStencilGeometry5D geom(Coarse5d);
   LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
-  LittleDiracOp.CoarsenOperator(HermDefOp,Aggregates);
+  LittleDiracOperator LittleDiracOpCol(geom,FGrid,Coarse5d);
-  //  LittleDiracOp.CoarsenOperator(Diagonal,Aggregates);
 
-  std::cout<<GridLogMessage<<"Coarsened operator "<<std::endl;
+  HermOpAdaptor<LatticeFermionD> HOA(HermDefOp);
 
+  int pp=16;
+  //  LittleDiracOpCol.CoarsenOperator(HOA,Aggregates);
+  //  std::cout << "LittleDiracOp old " << LittleDiracOpCol._A[pp]<<std::endl;
+  LittleDiracOp.CoarsenOperatorColoured(HOA,Aggregates);
+  //  std::cout << "LittleDiracOp new " << LittleDiracOp._A[pp]<<std::endl;
+  
+  ///////////////////////////////////////////////////
+  // Test the operator
+  ///////////////////////////////////////////////////
   CoarseVector c_src (Coarse5d);
   CoarseVector c_res (Coarse5d);
+  CoarseVector c_res_dag(Coarse5d);
   CoarseVector c_proj(Coarse5d);
 
   subspace=Aggregates.subspace;
 
-  Complex one(1.0);
+  //  random(CRNG,c_src);
-  c_src = one;  // 1 in every element for vector 1.
+  c_src = 1.0;
-  Coordinate coor(5,0);
 
-  std::cout << "c_src"<< c_src<< std::endl;
   blockPromote(c_src,err,subspace);
 
   prom=Zero();
@@ -235,25 +173,23 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
   std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
 
-  //  blockPick(Coarse5d,c_src,c_src,coor);
-  //  blockPromote(c_src,prom,subspace);
-
-  //  Diagonal.HermOp(prom,tmp);
   HermDefOp.HermOp(prom,tmp);
-  //  HermDefOp.Op(prom,tmp);
+
   blockProject(c_proj,tmp,subspace);
   std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
 
   LittleDiracOp.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
+  LittleDiracOp.Mdag(c_src,c_res_dag);
 
   std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
+  std::cout<<GridLogMessage<<"Little dop dag : "<<norm2(c_res_dag)<<std::endl;
   std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
+
   c_proj = c_proj - c_res;
   std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
+
+  c_res_dag = c_res_dag - c_res;
+  std::cout<<GridLogMessage<<"Little dopDag - dop: "<<norm2(c_res_dag)<<std::endl;
 
   std::cout<<GridLogMessage << "Testing Hermiticity stochastically "<< std::endl;
   CoarseVector phi(Coarse5d);
@@ -267,8 +203,10 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage<<"Made randoms "<<norm2(phi)<<" " << norm2(chi)<<std::endl;
 
   LittleDiracOp.M(phi,Aphi);
+
   LittleDiracOp.Mdag(chi,Achi);
-  std::cout<<GridLogMessage<<"Aphi "<<norm2(Aphi)<<" Adag chi" << norm2(Achi)<<std::endl;
+
+  std::cout<<GridLogMessage<<"Aphi "<<norm2(Aphi)<<" A chi" << norm2(Achi)<<std::endl;
 
   ComplexD pAc = innerProduct(chi,Aphi);
   ComplexD cAp = innerProduct(phi,Achi);
@@ -276,7 +214,6 @@ int main (int argc, char ** argv)
   ComplexD pAp = innerProduct(phi,Aphi);
 
   std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
-
   std::cout<<GridLogMessage<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
 
   std::cout<<GridLogMessage<<"Testing linearity"<<std::endl;
@@ -295,34 +232,6 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  
-  PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM(Ddwf,Dpv);
-  HermOpAdaptor<LatticeFermionD> HOA(PVdagM);
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;   PM(HOA,src);
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  Subspace AggregatesPD(Coarse5d,FGrid,cb);
-  AggregatesPD.CreateSubspaceChebyshev(RNG5,
-				     HOA,
-				     nbasis,
-				     5000.0,
-				     0.02,
-				     100,
-				     50,
-				     50,
-				     0.0);
-  
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
   
   Grid_finalize();
   return 0;

tests/debug/Test_general_coarse_hdcg.cc

@@ -0,0 +1,272 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_general_coarse_hdcg.cc
+
+    Copyright (C) 2023
+
+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 */
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+#include <Grid/algorithms/GeneralCoarsenedMatrix.h>
+#include <Grid/algorithms/iterative/AdefGeneric.h>
+
+using namespace std;
+using namespace Grid;
+
+template<class Field> class TestSolver : public LinearFunction<Field> {
+public:
+  TestSolver() {};
+  void operator() (const Field &in, Field &out){    out = Zero();  }     
+};
+
+
+RealD InverseApproximation(RealD x){
+  return 1.0/x;
+}
+
+// Want Op in CoarsenOp to call MatPcDagMatPc
+template<class Field>
+class HermOpAdaptor : public LinearOperatorBase<Field>
+{
+  LinearOperatorBase<Field> & wrapped;
+public:
+  HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme)  {};
+  void Op     (const Field &in, Field &out)   { wrapped.HermOp(in,out);  }
+  void HermOp(const Field &in, Field &out)    { wrapped.HermOp(in,out); }
+  void AdjOp     (const Field &in, Field &out){ wrapped.HermOp(in,out);  }
+  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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
+};
+template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
+{
+public:
+  using LinearFunction<Field>::operator();
+  typedef LinearOperatorBase<Field> FineOperator;
+  FineOperator   & _SmootherOperator;
+  Chebyshev<Field> Cheby;
+  ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
+    _SmootherOperator(SmootherOperator),
+    Cheby(_lo,_hi,_ord,InverseApproximation)
+  {
+    std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
+  };
+  void operator() (const Field &in, Field &out) 
+  {
+    Field tmp(in.Grid());
+    tmp = in;
+    Cheby(_SmootherOperator,tmp,out);
+  }
+};
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls=16;
+
+  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);
+
+  // Construct a coarsened grid with 4^4 cell
+  Coordinate clatt = GridDefaultLatt();
+  for(int d=0;d<clatt.size();d++){
+    clatt[d] = clatt[d]/4;
+  }
+  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt,
+							    GridDefaultSimd(Nd,vComplex::Nsimd()),
+							    GridDefaultMpi());;
+  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+
+  ///////////////////////// RNGs /////////////////////////////////
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  std::vector<int> cseeds({5,6,7,8});
+
+  GridParallelRNG          RNG5(FGrid);   RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
+
+  ///////////////////////// Configuration /////////////////////////////////
+  LatticeGaugeField Umu(UGrid);
+
+  FieldMetaData header;
+  std::string file("ckpoint_lat.4000");
+  NerscIO::readConfiguration(Umu,header,file);
+
+  //////////////////////// Fermion action //////////////////////////////////
+  RealD mass=0.01;
+  RealD M5=1.8;
+  RealD b=1.5;
+  RealD c=0.5;
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+
+  SchurDiagMooeeOperator<MobiusFermionD, LatticeFermion> HermOpEO(Ddwf);
+
+  typedef HermOpAdaptor<LatticeFermionD> HermFineMatrix;
+  HermFineMatrix FineHermOp(HermOpEO);
+  
+  LatticeFermion result(FrbGrid); result=Zero();
+
+  LatticeFermion    src(FrbGrid); random(RNG5,src);
+
+  // Run power method on FineHermOp
+  PowerMethod<LatticeFermion>       PM;   PM(HermOpEO,src);
+
+ 
+  ////////////////////////////////////////////////////////////
+  ///////////// Coarse basis and Little Dirac Operator ///////
+  ////////////////////////////////////////////////////////////
+  const int nbasis = 40;
+  const int cb = 0 ;
+  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
+  typedef LittleDiracOperator::CoarseVector CoarseVector;
+
+  NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
+  
+  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FrbGrid,cb);
+  Aggregates.CreateSubspaceChebyshev(RNG5,
+				     HermOpEO,
+				     nbasis,
+				     //				     100.0,
+				     //	0.1, // Low pass is pretty high still -- 311 iters
+				     //				     250.0,
+				     //				     0.01, // subspace too low filter power wrong
+				     //				     250.0,
+				     //				     0.2, // slower
+				     95.0,
+				     //				     0.05, // nbasis 12 - 311 -- wrong coarse inv
+				     //				     0.05, // nbasis 12 - 154 -- right filt
+				     //				     0.1, // nbasis 12 - 169 oops
+				     //				     0.05, // nbasis 16 -- 127 iters
+				     //				     0.03, // nbasis 16 -- 13-
+				     //				     0.1,  // nbasis 16 -- 142; sloppy solve
+				     0.1,  // nbasis 24 
+				     300);
+  ////////////////////////////////////////////////////////////
+  // Need to check about red-black grid coarsening
+  ////////////////////////////////////////////////////////////
+  LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d);
+  LittleDiracOp.CoarsenOperatorColoured(FineHermOp,Aggregates);
+
+  typedef HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermMatrix;
+  HermMatrix CoarseOp (LittleDiracOp);
+
+  //////////////////////////////////////////
+  // Build a coarse lanczos
+  //////////////////////////////////////////
+  Chebyshev<CoarseVector>      IRLCheby(0.02,50.0,71);  // 1 iter
+  FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,CoarseOp);
+  PlainHermOp<CoarseVector>    IRLOp    (CoarseOp);
+  int Nk=64;
+  int Nm=128;
+  int Nstop=Nk;
+  ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-5,20);
+
+  int Nconv;
+  std::vector<RealD>            eval(Nm);
+  std::vector<CoarseVector>     evec(Nm,Coarse5d);
+  CoarseVector c_src(Coarse5d); c_src=1.0;
+  IRL.calc(eval,evec,c_src,Nconv);
+  DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
+  
+  //////////////////////////////////////////
+  // Build a coarse space solver
+  //////////////////////////////////////////
+  int maxit=20000;
+  ConjugateGradient<CoarseVector>  CG(1.0e-8,maxit,false);
+  ConjugateGradient<LatticeFermionD>  CGfine(1.0e-8,10000,false);
+  ZeroGuesser<CoarseVector> CoarseZeroGuesser;
+
+  //  HPDSolver<CoarseVector> HPDSolve(CoarseOp,CG,CoarseZeroGuesser);
+  HPDSolver<CoarseVector> HPDSolve(CoarseOp,CG,DeflCoarseGuesser);
+
+  //////////////////////////////////////////
+  // Build a smoother
+  //////////////////////////////////////////
+  //  ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(10.0,100.0,10,FineHermOp); //499
+  //  ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(3.0,100.0,10,FineHermOp);  //383
+  //  ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(1.0,100.0,10,FineHermOp);  //328
+  //  std::vector<RealD> los({0.5,1.0,3.0}); // 147/142/146 nbasis 1
+  //  std::vector<RealD> los({1.0,2.0}); // Nbasis 24: 88,86 iterations
+  //  std::vector<RealD> los({2.0,4.0}); // Nbasis 32 == 52, iters
+  //  std::vector<RealD> los({2.0,4.0}); // Nbasis 40 == 36,36 iters
+
+  //
+  // Turns approx 2700 iterations into 340 fine multiplies with Nbasis 40
+  // Need to measure cost of coarse space.
+  //
+  // -- i) Reduce coarse residual   -- 0.04
+  // -- ii) Lanczos on coarse space -- done
+  // -- iii) Possible 1 hop project and/or preconditioning it - easy - PrecCG it and
+  //         use a limited stencil. Reread BFM code to check on evecs / deflation strategy with prec
+  //
+  std::vector<RealD> los({3.0}); // Nbasis 40 == 36,36 iters
+  std::vector<int> ords({7,8,10}); // Nbasis 40 == 40,38,36 iters (320,342,396 mults)
+
+  // Standard CG
+  //      result=Zero();
+  //      CGfine(HermOpEO, src, result);
+
+  for(int l=0;l<los.size();l++){
+
+    RealD lo = los[l];
+
+    for(int o=0;o<ords.size();o++){
+
+      ConjugateGradient<CoarseVector>  CGsloppy(4.0e-2,maxit,false);
+      HPDSolver<CoarseVector> HPDSolveSloppy(CoarseOp,CGsloppy,DeflCoarseGuesser);
+      
+      //    ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(lo,92,10,FineHermOp); // 36 best case
+      ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(lo,92,ords[o],FineHermOp);  // 311
+
+      //////////////////////////////////////////
+      // Build a HDCG solver
+      //////////////////////////////////////////
+      TwoLevelFlexiblePcg<LatticeFermion,CoarseVector,Subspace>
+	HDCG(1.0e-8, 3000,
+	     FineHermOp,
+	     Smoother,
+	     HPDSolveSloppy,
+	     HPDSolve,
+	     Aggregates);
+
+      //    result=Zero();
+      //    HDCG(src,result);
+    
+      result=Zero();
+      HDCG.Inflexible(src,result);
+    }
+  }
+  
+  Grid_finalize();
+  return 0;
+}

tests/debug/Test_general_coarse_pvdagm.cc

@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -141,7 +141,7 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  const int Ls=16;
+  const int Ls=2;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
@@ -152,7 +152,7 @@ int main (int argc, char ** argv)
   // Construct a coarsened grid
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
+    clatt[d] = clatt[d]/4;
   }
   GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
   GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
@@ -174,6 +174,7 @@ int main (int argc, char ** argv)
   FieldMetaData header;
   std::string file("ckpoint_lat.4000");
   NerscIO::readConfiguration(Umu,header,file);
+  //Umu = 1.0;
   
   RealD mass=0.5;
   RealD M5=1.8;
@@ -181,7 +182,7 @@ int main (int argc, char ** argv)
   DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
 
-  const int nbasis = 20;
+  const int nbasis = 1;
   const int cb = 0 ;
   LatticeFermion prom(FGrid);
 
@@ -216,6 +217,47 @@ int main (int argc, char ** argv)
   LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
   LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
 
+  std::cout<<GridLogMessage<<std::endl;
+  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
+  std::cout<<GridLogMessage<<std::endl;
+  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
+  
+  CoarseVector c_src (Coarse5d);
+  CoarseVector c_res (Coarse5d);
+  CoarseVector c_proj(Coarse5d);
+
+  std::vector<LatticeFermion> subspace(nbasis,FGrid);
+  subspace=AggregatesPD.subspace;
+
+  Complex one(1.0);
+  c_src = one;  // 1 in every element for vector 1.
+  blockPromote(c_src,err,subspace);
+
+  prom=Zero();
+  for(int b=0;b<nbasis;b++){
+    prom=prom+subspace[b];
+  }
+  err=err-prom; 
+  std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
+  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
+  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
+
+  PVdagM.Op(prom,tmp);
+  blockProject(c_proj,tmp,subspace);
+  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
+
+  LittleDiracOpPV.M(c_src,c_res);
+  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
+
+  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
+  //  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
+
+  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
+  //  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
+  c_proj = c_proj - c_res;
+  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
+  //  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
+
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
   std::cout<<GridLogMessage<<std::endl;