Working getting closer to HDCG but some low level engineering work still needed

+ MUCH work on optimisation

Grid/algorithms/GeneralCoarsenedMatrix.h

@@ -62,7 +62,7 @@ public:
 // Need to worry about red-black now
 class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
 public:
-  NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(2)
+  NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(1)
   {
       this->BuildShifts();
   };
@@ -95,7 +95,7 @@ public:
 // Need to worry about red-black now
 class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
 public:
-  NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(4)
+  NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(1)
   {
     this->BuildShifts();
   };
@@ -142,7 +142,7 @@ public:
 };
 class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
 public:
-  NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(2)
+  NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(1)
   {
       this->BuildShifts();
   };
@@ -175,7 +175,7 @@ public:
 // Need to worry about red-black now
 class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
 public:
-  NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(4)
+  NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(1)
   {
     this->BuildShifts();
   };
@@ -239,7 +239,7 @@ public:
   // Data members
   ////////////////////
   int hermitian;
-  GridCartesian *       _FineGrid; 
+  GridBase      *       _FineGrid; 
   GridCartesian *       _CoarseGrid; 
   NonLocalStencilGeometry &geom;
   PaddedCell Cell;
@@ -251,11 +251,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 +263,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 +294,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);
@@ -318,17 +332,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++) {
@@ -347,7 +359,7 @@ public:
 	    Aggregation<Fobj,CComplex,nbasis> & Subspace)
   {
     // Create a random
-    GridCartesian *grid = FineGrid();
+    GridBase *grid = FineGrid();
     FineField MbV(grid);
     FineField tmp(grid);
     FineField f_src(grid);
@@ -370,8 +382,12 @@ public:
   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++){
@@ -393,19 +409,28 @@ 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);
 	  for(int bb=0;bb<nbasis;bb++){
@@ -415,14 +440,77 @@ public:
 	  pokeSite(Ab,_A[p],scoor);
 	  pokeSite(Adagb,_Adag[p],bcoor);
 	}
+	tpeek+=usecond();
       }
     }
-    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]);
+    }
+    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;
+  }
+
+#if 0
+  void CoarsenOperatorColoured(LinearOperatorBase<Lattice<Fobj> > &linop,
+		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
+  {
+    std::cout << GridLogMessage<< "CoarsenMatrixColoured "<< std::endl;
+    GridBase *grid = FineGrid();
+
+    /////////////////////////////////////////////////////////////
+    // Orthogonalise the subblocks over the basis
+    /////////////////////////////////////////////////////////////
+    CoarseScalar InnerProd(CoarseGrid()); 
+    blockOrthogonalise(InnerProd,Subspace.subspace);
+
+    const int npoint = geom.npoint
+      
+    // Now compute the matrix elements of linop between this orthonormal
+    // set of vectors.
+    FineField bV(grid);
+    FineField MbV(grid);
+    FineField tmp(grid);
+    CoarseVector coarseInner(CoarseGrid());
+    Coordinate clatt = CoarseGrid()->GlobalDimensions();
+    Coordinate steclatt = CoarseGrid()->GlobalDimensions();
+    for(int v=0;v<nbasis;v++){
+      for(int p=0;p<npoint;p++){ // Loop over momenta
+	/////////////////////////////////////////////////////
+	// Stick a different phase on every block
+	/////////////////////////////////////////////////////
+	CoarseVector pha(CoarseGrid());
+	std::vector<double> dmom = 
+	for(int mu=0;mu<5;mu++){
+	  dmom[mu] = imom[mu]*2*M_PI/global_nb[mu];
+	}
+	/*
+	 *     Here, k,l index which possible momentum/shift within the N-points connected by MdagM.
+	 *     Matrix index i is mapped to this shift via 
+	 *               GetDelta(imat_ball_idx[i]).
+	 *
+	 *     conj(phases[block]) proj[k][ block*Nvec+j ] =  \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}
+	 */
+      }
+    }
+    
     for(int p=0;p<geom.npoint;p++){
       _A[p] = Cell.Exchange(_A[p]);
       _Adag[p] = Cell.Exchange(_Adag[p]);
     }
   }
+#endif
+
   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);};