ADEF1 and 1 hop projection

Grid/algorithms/GeneralCoarsenedMatrix.h

@@ -191,7 +191,8 @@ public:
 template<class Fobj,class CComplex,int nbasis>
 class GeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
 public:
-    
+
+  typedef GeneralCoarsenedMatrix<Fobj,CComplex,nbasis> GeneralCoarseOp;
   typedef iVector<CComplex,nbasis >           siteVector;
   typedef iMatrix<CComplex,nbasis >           siteMatrix;
   typedef Lattice<iScalar<CComplex> >         CoarseComplexField;
@@ -222,27 +223,22 @@ public:
   GridCartesian * CoarseGrid(void)     { return _CoarseGrid; };   // this is all the linalg routines need to know
 
 
-  void ProjectNearestNeighbour(RealD shift)
+  void ProjectNearestNeighbour(RealD shift, GeneralCoarseOp &CopyMe)
   {
-    int Nd = geom.grid->Nd();
+    int nfound=0;
-    int point;
-    std::cout << "ProjectNearestNeighbour "<<std::endl;
     for(int p=0;p<geom.npoint;p++){
-      int nhops = 0;
+      for(int pp=0;pp<CopyMe.geom.npoint;pp++){
-      for(int s=0;s<Nd;s++){
+ 	// Search for the same relative shift
-	nhops+=abs(geom.shifts[p][s]);
+	// Avoids brutal handling of Grid pointers
-      }
+	if ( CopyMe.geom.shifts[pp]==geom.shifts[p] ) {
-      if(nhops>1) {
+	  _A[p] = CopyMe.Cell.Extract(CopyMe._A[pp]);
-	std::cout << "setting geom "<<p<<" shift "<<geom.shifts[p]<<" to zero "<<std::endl;
+	  _Adag[p] = CopyMe.Cell.Extract(CopyMe._Adag[pp]);
-	_A[p]=Zero();
+	  ExchangeCoarseLinks();
-	_Adag[p]=Zero();
+	  nfound++;
-      }
+	}
-      if(nhops==0) {
-	std::cout << " Adding IR shift "<<shift<<" to "<<geom.shifts[p]<<std::endl;
-	_A[p]=_A[p]+shift;
-	_Adag[p]=_Adag[p]+shift;
       }
     }
+    assert(nfound==geom.npoint);
   }
   
   GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid)
@@ -317,27 +313,13 @@ public:
       autoView( out_v , pout, AcceleratorWrite);
       autoView( Stencil_v  , Stencil, AcceleratorRead);
       tviews+=usecond();
-
-      std::cout << "Calling accelerator for loop " <<std::endl;
       
       for(int point=0;point<npoint;point++){
 	tviews-=usecond();
 	autoView( A_v, A[point],AcceleratorRead);
 	tviews+=usecond();
 	tmult-=usecond();
-#if 0
+	accelerator_for(sss, osites*nbasis, Nsimd, {
-	prof_accelerator_for(ss, osites, Nsimd, {
-	  // Junk load is annoying -- need to sort out the types better.
-	  //////////////////////////////
-	  // GPU chokes on gpermute - want coalescedReadPermute()
-	  //	gpermute(nbr,SE->_permute);
-	  //////////////////////////////
-	  auto SE = Stencil_v.GetEntry(point,ss);
-	  int o = SE->_offset;
-	  coalescedWrite(out_v[ss],out_v(ss) + A_v(ss)*in_v(o));
-	});
-#else
-	prof_accelerator_for(sss, osites*nbasis, Nsimd, {
 
 	    typedef decltype(coalescedRead(in_v[0]))    calcVector;
 
@@ -352,10 +334,9 @@ public:
 	    }
 	    coalescedWrite(out_v[ss](b),res);
 	});
-#endif
+
 	tmult+=usecond();
       }
-      std::cout << "Called accelerator for loop " <<std::endl;
     }
     text-=usecond();
     out = Cell.Extract(pout);
@@ -392,84 +373,6 @@ 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;
-    GridBase *grid = FineGrid();
-
-    ////////////////////////////////////////////////
-    // Orthogonalise the subblocks over the basis
-    ////////////////////////////////////////////////
-    CoarseScalar InnerProd(CoarseGrid()); 
-    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);
-    CoarseVector coarseInner(CoarseGrid());
-    
-    // Very inefficient loop of order coarse volume.
-    // First pass hack
-    // Could replace with a coloring scheme our phase scheme
-    // as in BFM
-    for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
-      Coordinate bcoor;
-      CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
-
-      for(int b=0;b<nbasis;b++){
-	tpick-=usecond();
-	blockPick(CoarseGrid(),Subspace.subspace[b],bV,bcoor);
-	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);
-	  auto Ab    = peekSite(_A[p],scoor);
-	  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[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);
-    }
-    ExchangeCoarseLinks();
-    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

Grid/algorithms/iterative/AdefGeneric.h

@@ -43,7 +43,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 NAMESPACE_BEGIN(Grid);
 
 template<class Field, class CoarseField, class Aggregation>
-class TwoLevelFlexiblePcg : public LinearFunction<Field>
+class TwoLevelADEF2 : public LinearFunction<Field>
 {
  public:
   RealD   Tolerance;
@@ -64,14 +64,14 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
   Aggregation &_Aggregates;                    
   
   // more most opertor functions
-  TwoLevelFlexiblePcg(RealD tol,
+  TwoLevelADEF2(RealD tol,
-		      Integer maxit,
+		   Integer maxit,
-		      LinearOperatorBase<Field>   &FineLinop,
+		   LinearOperatorBase<Field>   &FineLinop,
-		      LinearFunction<Field>   &Smoother,
+		   LinearFunction<Field>   &Smoother,
-		      LinearFunction<CoarseField>  &CoarseSolver,
+		   LinearFunction<CoarseField>  &CoarseSolver,
-		      LinearFunction<CoarseField>  &CoarseSolverPrecise,
+		   LinearFunction<CoarseField>  &CoarseSolverPrecise,
-		      Aggregation &Aggregates
+		   Aggregation &Aggregates
-		     ) : 
+		   ) : 
       Tolerance(tol), 
       MaxIterations(maxit),
       _FineLinop(FineLinop),
@@ -84,7 +84,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     grid       = Aggregates.FineGrid;
   };
   
-  void Inflexible(const Field &src,Field &psi)
+  virtual void operator() (const Field &src, Field &psi)
   {
     Field resid(grid);
     RealD f;
@@ -192,10 +192,6 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     return ;
   }
   
-  virtual void operator() (const Field &in, Field &out)
-  {
-    this->Inflexible(in,out);
-  }
 
  public:
 
@@ -285,5 +281,49 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
   }
 };
 
+template<class Field, class CoarseField, class Aggregation>
+class TwoLevelADEF1ev : public TwoLevelADEF2<Field,CoarseField,Aggregation>
+{
+  TwoLevelADEF1ev(RealD tol,
+		   Integer maxit,
+		   LinearOperatorBase<Field>   &FineLinop,
+		   LinearFunction<Field>   &Smoother,
+		   LinearFunction<CoarseField>  &CoarseSolver,
+		   LinearFunction<CoarseField>  &CoarseSolverPrecise,
+		   Aggregation &Aggregates ) : 
+    TwoLevelADEF2<Field,CoarseField,Aggregation>(tol,maxit,FineLinop,Smoother,CoarseSolver,CoarseSolverPrecise,Aggregates)
+  {};
+
+  // Can just inherit existing Vstart
+  // Can just inherit existing Vout
+  // Can just inherit existing M2
+  // Can just inherit existing M3
+  // Override PcgM1
+  virtual void PcgM1(Field & in, Field & out)
+  {
+    CoarseField PleftProj(this->coarsegrid);
+    CoarseField PleftMss_proj(this->coarsegrid);
+    Field tmp(this->grid);
+    Field Pin(this->grid);
+
+    //MP  + Q = M(1-AQ) + Q = M
+    // // If we are eigenvector deflating in coarse space
+    // // Q   = Sum_i |phi_i> 1/lambda_i <phi_i|
+    // // A Q = Sum_i |phi_i> <phi_i|
+    // // M(1-AQ) = M(1-proj) + Q
+    this->_Aggregates.ProjectToSubspace(PleftProj,in);     
+    this->_Aggregates.PromoteFromSubspace(PleftProj,tmp);// tmp = Qin
+
+    Pin = in - tmp;
+    this->_Smoother(Pin,out);
+
+    this->_CoarseSolver(PleftProj,PleftMss_proj); 
+    this->_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Qin
+
+    out = out + tmp;
+  }
+};
+
 NAMESPACE_END(Grid);
+
 #endif

tests/debug/Test_general_coarse_hdcg.cc

@@ -149,6 +149,7 @@ int main (int argc, char ** argv)
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
   NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
+  NearestStencilGeometry5D geom_nn(Coarse5d);
   
   // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
   typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
@@ -178,8 +179,8 @@ int main (int argc, char ** argv)
   LittleDiracOp.CoarsenOperatorColoured(FineHermOp,Aggregates);
 
   // Try projecting to one hop only
-  LittleDiracOperator LittleDiracOpProj(LittleDiracOp);
+  LittleDiracOperator LittleDiracOpProj(geom_nn,FrbGrid,Coarse5d);
-  LittleDiracOpProj.ProjectNearestNeighbour(0.5);
+  LittleDiracOpProj.ProjectNearestNeighbour(0.5,LittleDiracOp);
 
   typedef HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermMatrix;
   HermMatrix CoarseOp (LittleDiracOp);
@@ -260,7 +261,7 @@ int main (int argc, char ** argv)
       //////////////////////////////////////////
       // Build a HDCG solver
       //////////////////////////////////////////
-      TwoLevelFlexiblePcg<LatticeFermion,CoarseVector,Subspace>
+      TwoLevelADEF2<LatticeFermion,CoarseVector,Subspace>
 	HDCG(1.0e-8, 3000,
 	     FineHermOp,
 	     Smoother,
@@ -268,11 +269,8 @@ int main (int argc, char ** argv)
 	     HPDSolve,
 	     Aggregates);
 
-      //    result=Zero();
-      //    HDCG(src,result);
-    
       result=Zero();
-      HDCG.Inflexible(src,result);
+      HDCG(src,result);
     }
   }