Patches for beginnings of an overlap multigrid

lib/Algorithms.h

@@ -3,7 +3,6 @@
 
 #include <algorithms/SparseMatrix.h>
 #include <algorithms/LinearOperator.h>
-#include <algorithms/CoarsenedMatrix.h>
 
 #include <algorithms/approx/Zolotarev.h>
 #include <algorithms/approx/Chebyshev.h>
@@ -17,6 +16,9 @@
 
 #include <algorithms/iterative/ConjugateGradientMultiShift.h>
 
+
+#include <algorithms/CoarsenedMatrix.h>
+
 // Eigen/lanczos
 // EigCg
 // MCR

lib/algorithms/CoarsenedMatrix.h

@@ -12,9 +12,6 @@ namespace Grid {
     std::vector<int> directions   ;
     std::vector<int> displacements;
 
-    // FIXME -- don't like xposing the operator directions
-    // as different to the geometrical dirs
-    // Also don't like special casing five dim.. should pass an object in template
   Geometry(int _d)  {
   
       int base = (_d==5) ? 1:0;
@@ -64,6 +61,99 @@ namespace Grid {
 
   };
   
+  template<class Fobj,class CComplex,int nbasis>
+  class Aggregation   {
+  public:
+    typedef iVector<CComplex,nbasis >             siteVector;
+    typedef Lattice<siteVector>                 CoarseVector;
+    typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+
+    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
+    typedef Lattice<Fobj >        FineField;
+
+    GridBase *CoarseGrid;
+    GridBase *FineGrid;
+    std::vector<Lattice<Fobj> > subspace;
+
+    Aggregation(GridBase *_CoarseGrid,GridBase *_FineGrid) : 
+      CoarseGrid(_CoarseGrid),
+      FineGrid(_FineGrid),
+      subspace(nbasis,_FineGrid)
+	{
+	};
+  
+    void Orthogonalise(void){
+      CoarseScalar InnerProd(CoarseGrid); 
+      blockOrthogonalise(InnerProd,subspace);
+#if 1
+      //      CheckOrthogonal();
+#endif
+
+    } 
+    void CheckOrthogonal(void){
+      CoarseVector iProj(CoarseGrid); 
+      CoarseVector eProj(CoarseGrid); 
+      Lattice<CComplex> pokey(CoarseGrid);
+
+      
+      for(int i=0;i<nbasis;i++){
+	blockProject(iProj,subspace[i],subspace);
+
+	eProj=zero; 
+	for(int ss=0;ss<CoarseGrid->oSites();ss++){
+	  eProj._odata[ss](i)=CComplex(1.0);
+	}
+	eProj=eProj - iProj;
+	std::cout<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
+      }
+      std::cout <<"CheckOrthog done"<<std::endl;
+    }
+    void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
+      blockProject(CoarseVec,FineVec,subspace);
+    }
+    void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
+      blockPromote(CoarseVec,FineVec,subspace);
+    }
+    void CreateSubspaceRandom(GridParallelRNG &RNG){
+      for(int i=0;i<nbasis;i++){
+	random(RNG,subspace[i]);
+	std::cout<<" norm subspace["<<i<<"] "<<norm2(subspace[i])<<std::endl;
+      }
+      Orthogonalise();
+    }
+    virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop) {
+
+      RealD scale;
+
+      ConjugateGradient<FineField> CG(1.0e-4,10000);
+      FineField noise(FineGrid);
+      FineField Mn(FineGrid);
+
+      for(int b=0;b<nbasis;b++){
+	
+	gaussian(RNG,noise);
+	scale = std::pow(norm2(noise),-0.5); 
+	noise=noise*scale;
+
+	hermop.Op(noise,Mn); std::cout << "Noise    "<<b<<" nMdagMn "<<norm2(Mn)<< " "<< norm2(noise)<<std::endl;
+	for(int i=0;i<2;i++){
+
+	  CG(hermop,noise,subspace[b]);
+
+	  noise = subspace[b];
+	  scale = std::pow(norm2(noise),-0.5); 
+	  noise=noise*scale;
+
+	}
+
+	hermop.Op(noise,Mn); std::cout << "filtered    "<<b<<" <s|MdagM|s> "<<norm2(Mn)<< " "<< norm2(noise)<<std::endl;
+	subspace[b] = noise;
+      }
+
+      Orthogonalise();
+
+    }
+  };
   // Fine Object == (per site) type of fine field
   // nbasis      == number of deflation vectors
   template<class Fobj,class CComplex,int nbasis>
@@ -145,7 +235,8 @@ namespace Grid {
       comm_buf.resize(Stencil._unified_buffer_size);
     };
 
-    void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,std::vector<Lattice<Fobj> > & subspace){
+    void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
+			 Aggregation<Fobj,CComplex,nbasis> & Subspace){
 
       FineField iblock(FineGrid); // contributions from within this block
       FineField oblock(FineGrid); // contributions from outwith this block
@@ -162,8 +253,11 @@ namespace Grid {
       CoarseScalar InnerProd(Grid()); 
 
       // Orthogonalise the subblocks over the basis
-      blockOrthogonalise(InnerProd,subspace);
+      blockOrthogonalise(InnerProd,Subspace.subspace);
-      blockProject(iProj,subspace[0],subspace);
+      //Subspace.Orthogonalise();
+      //      Subspace.CheckOrthogonal();
+      //Subspace.Orthogonalise();
+      //      Subspace.CheckOrthogonal();
 
       // Compute the matrix elements of linop between this orthonormal
       // set of vectors.
@@ -177,7 +271,7 @@ namespace Grid {
       assert(self_stencil!=-1);
 
       for(int i=0;i<nbasis;i++){
-	phi=subspace[i];
+	phi=Subspace.subspace[i];
 	for(int p=0;p<geom.npoint;p++){ 
 
 	  int dir   = geom.directions[p];
@@ -210,8 +304,10 @@ namespace Grid {
 	    assert(0);
 	  }
 
-	  blockProject(iProj,iblock,subspace);
+	  Subspace.ProjectToSubspace(iProj,iblock);
-	  blockProject(oProj,oblock,subspace);
+	  Subspace.ProjectToSubspace(oProj,oblock);
+	  //	  blockProject(iProj,iblock,Subspace.subspace);
+	  //	  blockProject(oProj,oblock,Subspace.subspace);
 	  for(int ss=0;ss<Grid()->oSites();ss++){
 	    for(int j=0;j<nbasis;j++){
 	      if( disp!= 0 ) {
@@ -234,19 +330,19 @@ namespace Grid {
       }
       std::cout<< " picking by block0 "<< self_stencil <<std::endl;
 
-      phi=subspace[0];
+      phi=Subspace.subspace[0];
       std::vector<int> bc(FineGrid->_ndimension,0);
 
       blockPick(Grid(),phi,tmp,bc);      // Pick out a block
       linop.Op(tmp,Mphi);                // Apply big dop
-      blockProject(iProj,Mphi,subspace); // project it and print it
+      blockProject(iProj,Mphi,Subspace.subspace); // project it and print it
       std::cout<< " Computed matrix elements from block zero only "<<std::endl;
       std::cout<< iProj <<std::endl;
       std::cout<<"Computed Coarse Operator"<<std::endl;
 #endif
-      //      AssertHermitian();
       //      ForceHermitian();
-      //      ForceDiagonal();
+      AssertHermitian();
+      // ForceDiagonal();
     }
     void ForceDiagonal(void) {
 
@@ -263,7 +359,7 @@ namespace Grid {
 
       Complex one(1.0);
 
-      iMatrix<Complex,nbasis> ident;  ident=one;
+      iMatrix<CComplex,nbasis> ident;  ident=one;
 
       val = val*adj(val);
       val = val + 1.0;
@@ -279,7 +375,7 @@ namespace Grid {
 	int dd=d+1;
 	A[2*d] = adj(Cshift(A[2*d+1],dd,1));
       }
-      A[8] = 0.5*(A[8] + adj(A[8]));
+      //      A[8] = 0.5*(A[8] + adj(A[8]));
     }
     void AssertHermitian(void) {
       CoarseMatrix AA    (Grid());

lib/algorithms/iterative/ConjugateGradient.h

@@ -15,7 +15,7 @@ public:
     Integer MaxIterations;
     int verbose;
     ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { 
-      verbose=1;
+      verbose=0;
     };
 
 
@@ -58,7 +58,7 @@ public:
 	return;
       }
       
-      std::cout << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
+      if(verbose) std::cout << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
       
       int k;
       for (k=1;k<=MaxIterations;k++){
@@ -69,23 +69,29 @@ public:
 
 	RealD    qqck = norm2(mmp);
 	ComplexD dck  = innerProduct(p,mmp);
-	//	if (verbose) std::cout <<std::setprecision(4)<< "ConjugateGradient:  d,qq "<<d<< " "<<qq <<" qqcheck "<< qqck<< " dck "<< dck<<std::endl;
+	if (verbose) std::cout <<std::setprecision(4)<< "ConjugateGradient:  d,qq "<<d<< " "<<qq <<" qqcheck "<< qqck<< " dck "<< dck<<std::endl;
       
 	a      = c/d;
 	b_pred = a*(a*qq-d)/c;
 
 
-	//	if (verbose) std::cout <<std::setprecision(4)<< "ConjugateGradient:  a,bp "<<a<< " "<<b_pred <<std::endl;
+	if (verbose) std::cout <<std::setprecision(4)<< "ConjugateGradient:  a,bp "<<a<< " "<<b_pred <<std::endl;
 	cp = axpy_norm(r,-a,mmp,r);
 	b = cp/c;
-	//	std::cout <<std::setprecision(4)<< "ConjugateGradient:  cp,b "<<cp<< " "<<b <<std::endl;
+	if (verbose) std::cout <<std::setprecision(4)<< "ConjugateGradient:  cp,b "<<cp<< " "<<b <<std::endl;
 	
 	// Fuse these loops ; should be really easy
 	psi= a*p+psi;
 	p  = p*b+r;
 	  
 	if (verbose) std::cout<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
-
+	
+	if (0) { 
+	  Field tmp(src._grid);
+	  Linop.HermOpAndNorm(psi,tmp,qqck,qqck);
+	  tmp=tmp-src;
+	  std::cout<<"ConjugateGradient: true   residual  is "<< norm2(tmp);
+	}
 	// Stopping condition
 	if ( cp <= rsq ) { 
 	  
@@ -98,9 +104,10 @@ public:
 	  RealD resnorm = sqrt(norm2(p));
 	  RealD true_residual = resnorm/srcnorm;
 
-	  std::cout<<"ConjugateGradient: Converged on iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
+	  std::cout<<"ConjugateGradient: Converged on iteration " <<k
-	  std::cout<<"ConjugateGradient: true   residual  is "<<true_residual<<" sol "<<psinorm<<" src "<<srcnorm<<std::endl;
+		   <<" computed residual "<<sqrt(cp/ssq)
-	  std::cout<<"ConjugateGradient: target residual was "<<Tolerance<<std::endl;
+		   <<" true residual     "<<true_residual
+		   <<" target "<<Tolerance<<std::endl;
 	  return;
 	}
       }

lib/algorithms/iterative/ConjugateResidual.h

@@ -16,7 +16,7 @@ namespace Grid {
     int verbose;
 
     ConjugateResidual(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { 
-      verbose=1;
+      verbose=0;
     };
 
     void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi){
@@ -37,8 +37,8 @@ namespace Grid {
       Linop.HermOpAndNorm(p,Ap,pAp,pAAp);
       Linop.HermOpAndNorm(r,Ar,rAr,rAAr);
 
-      std::cout << "pAp, pAAp"<< pAp<<" "<<pAAp<<std::endl;
+      if(verbose) std::cout << "pAp, pAAp"<< pAp<<" "<<pAAp<<std::endl;
-      std::cout << "rAr, rAAr"<< rAr<<" "<<rAAr<<std::endl;
+      if(verbose) std::cout << "rAr, rAAr"<< rAr<<" "<<rAAr<<std::endl;
 
       cp =norm2(r);
       ssq=norm2(src);
@@ -63,15 +63,16 @@ namespace Grid {
 	axpy(p,b,p,r);
 	pAAp=axpy_norm(Ap,b,Ap,Ar);
 
-	std::cout<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
+	if(verbose) std::cout<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
 
 	if(cp<rsq) {
 	  Linop.HermOp(psi,Ap);
 	  axpy(r,-1.0,src,Ap);
 	  RealD true_resid = norm2(r);
-	  std::cout<<"ConjugateResidual: Converged on iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
+	  std::cout<<"ConjugateResidual: Converged on iteration " <<k
-	  std::cout<<"ConjugateResidual: true   residual  is "<<true_resid<<std::endl;
+		   << " computed residual "<<sqrt(cp/ssq)
-	  std::cout<<"ConjugateResidual: target residual was "<<Tolerance <<std::endl;
+	           << " true residual "<<true_resid
+	           << " target "       <<Tolerance <<std::endl;
 	  return;
 	}
 

lib/lattice/Lattice_ET.h

@@ -296,7 +296,7 @@ GRID_DEF_UNOP(toReal,UnaryToReal);
 GRID_DEF_UNOP(toComplex,UnaryToComplex);
 GRID_DEF_UNOP(abs  ,UnaryAbs); //abs overloaded in cmath C++98; DON'T do the abs-fabs-dabs-labs thing
 GRID_DEF_UNOP(sqrt ,UnarySqrt);
-GRID_DEF_UNOP(rsqrt,UnarySqrt);
+GRID_DEF_UNOP(rsqrt,UnaryRsqrt);
 GRID_DEF_UNOP(sin  ,UnarySin);
 GRID_DEF_UNOP(cos  ,UnaryCos);
 GRID_DEF_UNOP(log  ,UnaryLog);

lib/lattice/Lattice_transfer.h

@@ -168,7 +168,7 @@ inline void blockNormalise(Lattice<CComplex> &ip,Lattice<vobj> &fineX)
   GridBase *coarse = ip._grid;
   Lattice<vobj> zz(fineX._grid); zz=zero;
   blockInnerProduct(ip,fineX,fineX);
-  ip = rsqrt(ip);
+  ip = pow(ip,-0.5);
   blockZAXPY(fineX,ip,fineX,zz);
 }
 // useful in multigrid project;

lib/qcd/utils/SUn.h

@@ -522,22 +522,66 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
   }
 
   // reunitarise??
+  static void LieRandomize(GridParallelRNG     &pRNG,LatticeMatrix &out,double scale=1.0){
+    GridBase *grid = out._grid;
+    LatticeComplex ca (grid);
+    LatticeMatrix  lie(grid);
+    LatticeMatrix  la (grid);
+    Complex ci(0.0,scale);
+    Matrix ta;
+
+    lie=zero;
+    for(int a=0;a<generators();a++){
+
+      random(pRNG,ca); ca=real(ca)-0.5;
+      generator(a,ta);
+
+      la=ci*ca*ta;
+
+      lie = lie+la; // e^{i la ta}
+    }
+    taExp(lie,out);
+  }
+
+  static void HotConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
+    LatticeMatrix Umu(out._grid);
+    for(int mu=0;mu<Nd;mu++){
+      LieRandomize(pRNG,Umu,1.0);
+      pokeLorentz(out,Umu,mu);
+    }
+  }
+  static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
+    LatticeMatrix Umu(out._grid);
+    for(int mu=0;mu<Nd;mu++){
+      LieRandomize(pRNG,Umu,0.01);
+      pokeLorentz(out,Umu,mu);
+    }
+  }
+  static void ColdConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
+    LatticeMatrix Umu(out._grid);
+    Umu=1.0;
+    for(int mu=0;mu<Nd;mu++){
+      pokeLorentz(out,Umu,mu);
+    }
+  }
 
   static void taProj( const LatticeMatrix &in,  LatticeMatrix &out){
     out = Ta(in);
   }
   static void taExp( const LatticeMatrix &x,  LatticeMatrix &ex){ 
-    LatticeMatrix xn   = x;
 
+    LatticeMatrix xn(x._grid);
     RealD nfac = 1.0;
-    ex = 1+x; // 1+x
+
+    xn = x;
+    ex =xn+Complex(1.0); // 1+x
 
     // Do a 12th order exponentiation
-    for(int i= 2; i <= 12; ++i)
+    for(int i=2; i <= 12; ++i)
     {
-      nfac = nfac/i; 
+      nfac = nfac/RealD(i); //1/2, 1/2.3 ...
       xn   = xn * x; // x2, x3,x4....
-      ex  += xn*nfac;// x2/2!, x3/3!....
+      ex   = ex+ xn*nfac;// x2/2!, x3/3!....
     }
   }
 

tests/Test_cayley_coarsen_support.cc

@@ -83,20 +83,27 @@ int main (int argc, char ** argv)
   std::cout<<"Error "<<norm2(err)<<std::endl;
 
   const int nbasis = 2;
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
   LatticeFermion prom(FGrid);
 
-  for(int b=0;b<nbasis;b++){
+  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-    random(RNG5,subspace[b]);
+
-  }
+  std::cout<<"Calling Aggregation class" <<std::endl;
-  std::cout << "Computed randoms"<< std::endl;
+
+  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FGrid);
+  Aggregates.CreateSubspaceRandom(RNG5);
+
+  subspace=Aggregates.subspace;
+
+  std::cout << "Called aggregation class"<< std::endl;
 
   typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
   LittleDiracOperator LittleDiracOp(*Coarse5d);
 
-  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
+  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
   
   CoarseVector c_src (Coarse5d);
   CoarseVector c_res (Coarse5d);

tests/Test_cayley_ldop_cg.cc

@@ -1,4 +1,6 @@
 #include <Grid.h>
+#include <qcd/utils/WilsonLoops.h>
+#include <qcd/utils/SUn.h>
 
 using namespace std;
 using namespace Grid;
@@ -59,67 +61,84 @@ int main (int argc, char ** argv)
   LatticeFermion    ref(FGrid); ref=zero;
   LatticeFermion    tmp(FGrid);
   LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu);
+  LatticeGaugeField Umu(UGrid); 
 
 
+
+  //gaussian(RNG4,Umu);
   //random(RNG4,Umu);
+
   NerscField header;
   std::string file("./ckpoint_lat.400");
   readNerscConfiguration(Umu,header,file);
+  //  SU3::ColdConfiguration(RNG4,Umu);
+  //  SU3::TepidConfiguration(RNG4,Umu);
+  //  SU3::HotConfiguration(RNG4,Umu);
+  //  Umu=zero;
 
 #if 0
   LatticeColourMatrix U(UGrid);
-  U=zero;
   for(int nn=0;nn<Nd;nn++){
-    if (nn>2) {
+    U=peekIndex<LorentzIndex>(Umu,nn);
-      pokeIndex<LorentzIndex>(Umu,U,nn);
+    U=U*adj(U)-1.0;
-    }
+    std::cout<<"SU3 test "<<norm2(U)<<std::endl;
   }
 #endif  
   RealD mass=0.1;
-  RealD M5=1.0;
+  RealD M5=1.5;
 
   DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
 
   const int nbasis = 8;
+
+#if 0
   std::vector<LatticeFermion> subspace(nbasis,FGrid);
   LatticeFermion noise(FGrid);
   LatticeFermion ms(FGrid);
   for(int b=0;b<nbasis;b++){
-    Gamma g5(Gamma::Gamma5);
+
     gaussian(RNG5,noise);
     RealD scale = pow(norm2(noise),-0.5); 
     noise=noise*scale;
 
-    HermIndefOp.HermOp(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
+    HermIndefOp.Op(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-
 
     //    filter(HermIndefOp,noise,subspace[b]);
     // inverse iteration
     MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
-    ConjugateGradient<LatticeFermion> CG(1.0e-5,10000);
+    ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
 
-    for(int i=0;i<4;i++){
+    for(int i=0;i<1;i++){
 
       CG(HermDefOp,noise,subspace[b]);
       noise = subspace[b];
 
       scale = pow(norm2(noise),-0.5); 
       noise=noise*scale;
-      HermDefOp.HermOp(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
+      HermDefOp.Op(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
     }
 
     subspace[b] = noise;
-    HermIndefOp.HermOp(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
+    HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
+
   }
   std::cout << "Computed randoms"<< std::endl;
+#else
+  std::cout<<"Calling Aggregation class" <<std::endl;
+  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FGrid);
+  Aggregates.CreateSubspace(RNG5,HermDefOp);
+  std::cout << "Called aggregation class"<< std::endl;
+#endif
+
 
   typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
   LittleDiracOperator LittleDiracOp(*Coarse5d);
-  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
+  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
   
   CoarseVector c_src (Coarse5d);
   CoarseVector c_res (Coarse5d);
@@ -128,7 +147,7 @@ int main (int argc, char ** argv)
 
   std::cout << "Solving CG on coarse space "<< std::endl;
   MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
-  ConjugateGradient<CoarseVector> CG(1.0e-8,10000);
+  ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
   CG(PosdefLdop,c_src,c_res);
 
   std::cout << "Done "<< std::endl;

tests/Test_cayley_ldop_cr.cc

@@ -1,4 +1,6 @@
 #include <Grid.h>
+#include <qcd/utils/WilsonLoops.h>
+#include <qcd/utils/SUn.h>
 
 using namespace std;
 using namespace Grid;
@@ -59,67 +61,100 @@ int main (int argc, char ** argv)
   LatticeFermion    ref(FGrid); ref=zero;
   LatticeFermion    tmp(FGrid);
   LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu);
+  LatticeGaugeField Umu(UGrid); 
 
 
+
+  //gaussian(RNG4,Umu);
   //random(RNG4,Umu);
-  NerscField header;
-  std::string file("./ckpoint_lat.4000");
-  readNerscConfiguration(Umu,header,file);
 
-  RealD mass=0.01;
+  NerscField header;
-  RealD M5=1.8;
+  std::string file("./ckpoint_lat.400");
+  readNerscConfiguration(Umu,header,file);
+  //  SU3::ColdConfiguration(RNG4,Umu);
+  //  SU3::TepidConfiguration(RNG4,Umu);
+  //  SU3::HotConfiguration(RNG4,Umu);
+  //  Umu=zero;
+
+#if 0
+  LatticeColourMatrix U(UGrid);
+  for(int nn=0;nn<Nd;nn++){
+    U=peekIndex<LorentzIndex>(Umu,nn);
+    U=U*adj(U)-1.0;
+    std::cout<<"SU3 test "<<norm2(U)<<std::endl;
+  }
+#endif  
+  RealD mass=0.1;
+  RealD M5=1.5;
 
   DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
 
   const int nbasis = 8;
+
+#if 0
   std::vector<LatticeFermion> subspace(nbasis,FGrid);
   LatticeFermion noise(FGrid);
   LatticeFermion ms(FGrid);
   for(int b=0;b<nbasis;b++){
-    Gamma g5(Gamma::Gamma5);
+
     gaussian(RNG5,noise);
     RealD scale = pow(norm2(noise),-0.5); 
     noise=noise*scale;
 
-    HermIndefOp.HermOp(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
+    HermIndefOp.Op(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-
 
     //    filter(HermIndefOp,noise,subspace[b]);
     // inverse iteration
     MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
     ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
 
-    for(int i=0;i<4;i++){
+    for(int i=0;i<1;i++){
 
       CG(HermDefOp,noise,subspace[b]);
       noise = subspace[b];
 
       scale = pow(norm2(noise),-0.5); 
       noise=noise*scale;
-      HermDefOp.HermOp(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
+      HermDefOp.Op(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
     }
 
     subspace[b] = noise;
-    HermIndefOp.HermOp(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
+    HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
+
   }
   std::cout << "Computed randoms"<< std::endl;
+#else
+  std::cout<<"Calling Aggregation class" <<std::endl;
+  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FGrid);
+  Aggregates.CreateSubspace(RNG5,HermDefOp);
+  std::cout << "Called aggregation class"<< std::endl;
+#endif
+
 
   typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
   LittleDiracOperator LittleDiracOp(*Coarse5d);
-  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
+  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
   
   CoarseVector c_src (Coarse5d);
   CoarseVector c_res (Coarse5d);
   gaussian(CRNG,c_src);
   c_res=zero;
 
+
+  std::cout << "Solving CG on coarse space "<< std::endl;
+  MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
+  ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
+  CG(PosdefLdop,c_src,c_res);
+
+
   std::cout << "Solving MCR on coarse space "<< std::endl;
   HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
-  ConjugateResidual<CoarseVector> MCR(1.0e-8,10000);
+  ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
   MCR(HermIndefLdop,c_src,c_res);
 
   std::cout << "Done "<< std::endl;

tests/Test_cayley_ldop_cr_chebyshev.cc

@@ -1,124 +0,0 @@
-#include <Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-
-template<class d>
-struct scal {
-  d internal;
-};
-
-Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
-};
-
-double lowpass(double x)
-{
-  return pow(x*x,-2);
-}
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  Chebyshev<LatticeFermion> filter(-120.0, 120.0,256, lowpass);
-  ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
-  filter.csv(csv);
-  csv.close();
-
-  const int Ls=8;
-
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  // Construct a coarsened grid
-  std::vector<int> clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
-  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
-
-  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);
-
-  LatticeFermion    src(FGrid); gaussian(RNG5,src);
-  LatticeFermion result(FGrid); result=zero;
-  LatticeFermion    ref(FGrid); ref=zero;
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu);
-
-
-  //random(RNG4,Umu);
-  NerscField header;
-  std::string file("./ckpoint_lat.400");
-  readNerscConfiguration(Umu,header,file);
-
-#if 0
-  LatticeColourMatrix U(UGrid);
-  Complex cone(1.0,0.0);
-  for(int nn=0;nn<Nd;nn++){
-    if (nn==3) { 
-      U=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl;
-    //    else       { U[nn]= cone;std::cout << "unit gauge field in dir "<<nn<<std::endl; }
-      pokeIndex<LorentzIndex>(Umu,U,nn);
-    }
-  }
-#endif  
-  RealD mass=0.5;
-  RealD M5=1.8;
-
-  DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
-
-  const int nbasis = 8;
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-  LatticeFermion noise(FGrid);
-  LatticeFermion ms(FGrid);
-  for(int b=0;b<nbasis;b++){
-    Gamma g5(Gamma::Gamma5);
-    gaussian(RNG5,noise);
-    RealD scale = pow(norm2(noise),-0.5); 
-    noise=noise*scale;
-
-    HermIndefOp.HermOp(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-
-    filter(HermIndefOp,noise,subspace[b]);
-    scale = pow(norm2(subspace[b]),-0.5); 
-    subspace[b]=subspace[b]*scale;
-    HermIndefOp.HermOp(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
-  }
-  std::cout << "Computed randoms"<< std::endl;
-
-  typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  LittleDiracOperator LittleDiracOp(*Coarse5d);
-  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  gaussian(CRNG,c_src);
-  c_res=zero;
-
-  std::cout << "Solving MCR on coarse space "<< std::endl;
-  HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
-  ConjugateResidual<CoarseVector> MCR(1.0e-8,10000);
-  MCR(HermIndefLdop,c_src,c_res);
-
-  std::cout << "Done "<< std::endl;
-  Grid_finalize();
-}